| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // class G4GeomTestVolume implementation <<
27 // 26 //
28 // Author: G.Cosmo, CERN << 27 // $Id$
>> 28 //
>> 29 // --------------------------------------------------------------------
>> 30 // GEANT 4 class source file
>> 31 //
>> 32 // G4GeomTestVolume
>> 33 //
>> 34 // Author: D.C.Williams, UCSC (davidw@scipp.ucsc.edu)
29 // ------------------------------------------- 35 // --------------------------------------------------------------------
30 36
31 #include <queue> << 37 #include <vector>
32 #include <set> 38 #include <set>
>> 39 #include <algorithm>
>> 40 #include <iomanip>
33 41
34 #include "G4GeomTestVolume.hh" 42 #include "G4GeomTestVolume.hh"
>> 43
35 #include "G4PhysicalConstants.hh" 44 #include "G4PhysicalConstants.hh"
>> 45 #include "G4GeomTestLogger.hh"
>> 46 #include "G4GeomTestVolPoint.hh"
>> 47 #include "G4GeomTestSegment.hh"
>> 48
36 #include "G4VPhysicalVolume.hh" 49 #include "G4VPhysicalVolume.hh"
37 #include "G4LogicalVolume.hh" 50 #include "G4LogicalVolume.hh"
38 #include "G4VSolid.hh" 51 #include "G4VSolid.hh"
39 52
40 // 53 //
41 // Constructor 54 // Constructor
42 // 55 //
43 G4GeomTestVolume::G4GeomTestVolume( G4VPhysica << 56 G4GeomTestVolume::G4GeomTestVolume( const G4VPhysicalVolume *theTarget,
44 G4double t << 57 G4GeomTestLogger *theLogger,
45 G4int numb << 58 G4double theTolerance )
46 G4bool the << 59 : target(theTarget),
47 : target(theTarget), tolerance(theTolerance) << 60 logger(theLogger),
48 resolution(numberOfPoints), verbosity(theV << 61 tolerance(theTolerance),
>> 62 extent(theTarget->GetLogicalVolume()->GetSolid()->GetExtent())
49 {;} 63 {;}
50 64
51 // 65 //
52 // Destructor 66 // Destructor
53 // 67 //
54 G4GeomTestVolume::~G4GeomTestVolume() {;} 68 G4GeomTestVolume::~G4GeomTestVolume() {;}
55 69
56 // 70 //
57 // Get error tolerance 71 // Get error tolerance
58 // 72 //
59 G4double G4GeomTestVolume::GetTolerance() cons 73 G4double G4GeomTestVolume::GetTolerance() const
60 { 74 {
61 return tolerance; 75 return tolerance;
62 } 76 }
63 77
64 // 78 //
65 // Set error tolerance 79 // Set error tolerance
66 // 80 //
67 void G4GeomTestVolume::SetTolerance(G4double t << 81 void G4GeomTestVolume::SetTolerance(G4double val)
68 { 82 {
69 tolerance = tol; << 83 tolerance = val;
70 } 84 }
71 85
72 // 86 //
73 // Get number of points to check (resolution) << 87 // TestCartGridXYZ
74 // 88 //
75 G4int G4GeomTestVolume::GetResolution() const << 89 void G4GeomTestVolume::TestCartGridXYZ( G4int nx, G4int ny, G4int nz )
76 { 90 {
77 return resolution; << 91 TestCartGridX( ny, nz );
>> 92 TestCartGridY( nz, nx );
>> 93 TestCartGridZ( nx, ny );
78 } 94 }
79 95
80 // 96 //
81 // Set number of points to check (resolution) << 97 // TestCartGridX
82 // 98 //
83 void G4GeomTestVolume::SetResolution(G4int np) << 99 void G4GeomTestVolume::TestCartGridX( G4int ny, G4int nz )
84 { 100 {
85 resolution = np; << 101 TestCartGrid( G4ThreeVector(0,1,0), G4ThreeVector(0,0,1),
>> 102 G4ThreeVector(1,0,0), ny, nz );
86 } 103 }
87 104
88 // 105 //
89 // Get verbosity << 106 // TestCartGridY
90 // 107 //
91 G4bool G4GeomTestVolume::GetVerbosity() const << 108 void G4GeomTestVolume::TestCartGridY( G4int nz, G4int nx )
92 { 109 {
93 return verbosity; << 110 TestCartGrid( G4ThreeVector(0,0,1), G4ThreeVector(1,0,0),
>> 111 G4ThreeVector(0,1,0), nz, nx );
94 } 112 }
95 113
96 // 114 //
97 // Set verbosity << 115 // TestCartGridZ
98 // 116 //
99 void G4GeomTestVolume::SetVerbosity(G4bool ver << 117 void G4GeomTestVolume::TestCartGridZ( G4int nx, G4int ny )
100 { 118 {
101 verbosity = verb; << 119 TestCartGrid( G4ThreeVector(1,0,0), G4ThreeVector(0,1,0),
>> 120 G4ThreeVector(0,0,1), nx, ny );
102 } 121 }
103 122
104 // 123 //
105 // Get errors reporting threshold << 124 // TestRecursiveCartGrid
106 // 125 //
107 G4int G4GeomTestVolume::GetErrorsThreshold() c << 126 void G4GeomTestVolume::TestRecursiveCartGrid( G4int nx, G4int ny, G4int nz,
>> 127 G4int slevel, G4int depth )
108 { 128 {
109 return maxErr; << 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 }
110 } 180 }
111 181
112 // 182 //
113 // Set maximum number of errors to report << 183 // TestRecursiveCylinder
114 // 184 //
115 void G4GeomTestVolume::SetErrorsThreshold(G4in << 185 void
>> 186 G4GeomTestVolume::TestRecursiveCylinder( G4int nPhi, G4int nZ, G4int nRho,
>> 187 G4double fracZ, G4double fracRho,
>> 188 G4bool usePhi,
>> 189 G4int slevel, G4int depth )
116 { 190 {
117 maxErr = max; << 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 }
118 } 244 }
119 245
120 // 246 //
121 // Test overlap in tree << 247 // TestCylinder
122 // 248 //
123 void G4GeomTestVolume::TestOverlapInTree() con << 249 void G4GeomTestVolume::TestCylinder( G4int nPhi, G4int nZ, G4int nRho,
>> 250 G4double fracZ, G4double fracRho,
>> 251 G4bool usePhi )
124 { 252 {
125 std::queue<G4VPhysicalVolume*> volumes; << 253 //
126 std::set<G4LogicalVolume*> checked; << 254 // Get size of our volume
127 << 255 //
128 volumes.push(target); << 256 G4double xMax = std::max(extent.GetXmax(),-extent.GetXmin());
129 while (!volumes.empty()) << 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
130 { 274 {
131 G4VPhysicalVolume* current = volumes.front << 275 G4double cosPhi = std::cos(phi);
132 volumes.pop(); << 276 G4double sinPhi = std::sin(phi);
>> 277 G4ThreeVector vPhi1(sinPhi,-cosPhi,0);
133 278
134 // check overlaps for daughters << 279 //
135 G4LogicalVolume* logical = current->GetLog << 280 // Loop over rho
136 std::size_t ndaughters = logical->GetNoDau << 281 //
137 for (std::size_t i=0; i<ndaughters; ++i) << 282 G4double rho = rhoMax;
138 { << 283 G4int iRho = nRho;
139 G4VPhysicalVolume* daughter = logical->G << 284 do
140 daughter->CheckOverlaps(resolution, tole <<
141 } <<
142 <<
143 // append the queue of volumes <<
144 G4LogicalVolume* previousLogical = nullptr <<
145 for (std::size_t i=0; i<ndaughters; ++i) <<
146 { 285 {
147 G4VPhysicalVolume* daughter = logical->G << 286 G4ThreeVector p(rho*cosPhi,rho*sinPhi,0);
148 G4LogicalVolume* daughterLogical = daugh << 287 static G4ThreeVector v(0,0,1);
149 if (daughterLogical->GetNoDaughters() == << 288
150 G4bool found = (daughterLogical == previ << 289 TestOneLine( p, v );
151 if (!found) { found = (checked.find(daug << 290
152 if (!found) << 291 if (usePhi)
153 { <<
154 checked.emplace(daughterLogical); <<
155 previousLogical = daughterLogical; <<
156 volumes.push(daughter); <<
157 } <<
158 else <<
159 { 292 {
160 if (verbosity) << 293 //
>> 294 // Loop over z
>> 295 //
>> 296 G4double zScale = 1.0;
>> 297 G4int iZ=nZ;
>> 298 do
161 { 299 {
162 G4cout << "Checking overlaps in tree << 300 p.setZ( z0 + zScale*zHalfLength );
163 << " (" << daughterLogical->G << 301 TestOneLine(p,vPhi1);
164 << " is omitted, to avoid dup << 302 p.setZ( z0 - zScale*zHalfLength );
165 } << 303 TestOneLine(p,vPhi1);
>> 304 } while( zScale *= fracZ, --iZ );
166 } 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 //
>> 333 void 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()", "GeomNav0002",
>> 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 );
167 } 371 }
168 } 372 }
169 } << 373 }
170 374
171 // 375 //
172 // TestRecursiveOverlap << 376 // TestRecursiveLine
173 // 377 //
174 void G4GeomTestVolume::TestRecursiveOverlap( G << 378 void
>> 379 G4GeomTestVolume::TestRecursiveLine( const G4ThreeVector& p,
>> 380 const G4ThreeVector& v,
>> 381 G4int slevel, G4int depth )
175 { 382 {
176 // If reached requested level of depth (i.e. 383 // If reached requested level of depth (i.e. set to 0), exit.
177 // If not depth specified (i.e. set to -1), 384 // If not depth specified (i.e. set to -1), visit the whole tree.
178 // If requested initial level of depth is no 385 // If requested initial level of depth is not zero, visit from beginning
179 // 386 //
180 if (depth == 0) { return; } << 387 if (depth == 0) return;
181 if (depth != -1) { depth--; } << 388 if (depth != -1) depth--;
182 if (slevel != 0) { slevel--; } << 389 if (slevel != 0) slevel--;
183 390
184 // 391 //
185 // As long as we reached the requested << 392 // As long as we aren't a replica and we reached the requested
186 // initial level of depth, test ourselves 393 // initial level of depth, test ourselves
187 // 394 //
188 if ( slevel==0 ) << 395 if ( (!target->IsReplicated()) && (slevel==0) )
189 { 396 {
190 target->CheckOverlaps(resolution, toleranc << 397 TestOneLine( p, v );
>> 398 ReportErrors();
191 } 399 }
192 400
193 // 401 //
194 // Loop over unique daughters 402 // Loop over unique daughters
195 // 403 //
196 std::set<const G4LogicalVolume *> tested; 404 std::set<const G4LogicalVolume *> tested;
197 405
198 const G4LogicalVolume *logical = target->Get 406 const G4LogicalVolume *logical = target->GetLogicalVolume();
199 auto nDaughter = (G4int)logical->GetNoDaugh << 407 G4int nDaughter = logical->GetNoDaughters();
200 for( auto iDaughter=0; iDaughter<nDaughter; << 408 G4int iDaughter;
>> 409 for( iDaughter=0; iDaughter<nDaughter; ++iDaughter )
201 { 410 {
202 G4VPhysicalVolume *daughter = logical->Get << 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;
203 420
>> 421 //
204 // Tested already? 422 // Tested already?
205 // 423 //
206 // const G4LogicalVolume *daughterLogical << 424 std::pair<std::set<const G4LogicalVolume *>::iterator,G4bool>
207 // daughter->GetLogicalVolume(); << 425 there = tested.insert(daughterLogical);
208 // std::pair<std::set<const G4LogicalVolum << 426 if (!there.second) continue;
209 // there = tested.insert(daughterLog <<
210 // if (!there.second) continue; <<
211 427
212 // 428 //
213 // Recurse 429 // Recurse
214 // 430 //
215 G4GeomTestVolume vTest( daughter, toleranc << 431 G4GeomTestVolume vTest( daughter, logger, tolerance );
216 vTest.SetErrorsThreshold(maxErr); << 432 vTest.TestRecursiveLine( p, v, slevel, depth );
217 vTest.TestRecursiveOverlap( slevel,depth ) << 433 }
>> 434 }
>> 435
>> 436 //
>> 437 // TestOneLine
>> 438 //
>> 439 // Test geometry consistency along a single line
>> 440 //
>> 441 void 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 overshoot = result.first;
>> 601 }
>> 602
>> 603 if (kEntering)
>> 604 (*overshoot).second.AddError( points[i].GetPosition(),
>> 605 points[k].GetPosition() );
>> 606 else
>> 607 (*overshoot).second.AddError( points[k].GetPosition(),
>> 608 points[j].GetPosition() );
>> 609 }
>> 610 else
>> 611 {
>> 612 //
>> 613 // Ignore small overlaps if they are within tolerance
>> 614 //
>> 615 if (kS-iS < tolerance && (!kEntering)) continue;
>> 616 if (jS-kS < tolerance && kEntering ) continue;
>> 617
>> 618 //
>> 619 // We appear to overlap with another daughter
>> 620 //
>> 621 G4long key = iDaug < kDaug ?
>> 622 (iDaug*nDaughter + kDaug) : (kDaug*nDaughter + iDaug);
>> 623
>> 624 std::map<G4long,G4GeomTestOverlapList>::iterator overlap =
>> 625 overlaps.find(key);
>> 626 if (overlap == overlaps.end())
>> 627 {
>> 628 std::pair<std::map<G4long,G4GeomTestOverlapList>::iterator,G4bool>
>> 629 result =
>> 630 overlaps.insert( std::pair<const G4long,G4GeomTestOverlapList>
>> 631 (key,G4GeomTestOverlapList(target,iDaug,kDaug)) );
>> 632 overlap = result.first;
>> 633 }
>> 634
>> 635 if (points[k].Entering())
>> 636 (*overlap).second.AddError( points[k].GetPosition(),
>> 637 points[j].GetPosition() );
>> 638 else
>> 639 (*overlap).second.AddError( points[i].GetPosition(),
>> 640 points[k].GetPosition() );
>> 641 }
>> 642 }
218 } 643 }
>> 644 }
>> 645
>> 646 //
>> 647 // ReportErrors
>> 648 //
>> 649 void G4GeomTestVolume::ReportErrors()
>> 650 {
>> 651 //
>> 652 // Report overshoots
>> 653 //
>> 654 if (overshoots.empty())
>> 655 logger->NoProblem("GeomTest: no daughter volume extending outside mother detected.");
>> 656 else
>> 657 {
>> 658 std::map<G4long,G4GeomTestOvershootList>::iterator overshoot =
>> 659 overshoots.begin();
>> 660 while( overshoot != overshoots.end() )
>> 661 {
>> 662 logger->OvershootingDaughter( &(*overshoot).second );
>> 663 ++overshoot;
>> 664 }
>> 665 }
>> 666
>> 667 //
>> 668 // Report overlaps
>> 669 //
>> 670 if (overlaps.empty())
>> 671 logger->NoProblem("GeomTest: no overlapping daughters detected.");
>> 672 else
>> 673 {
>> 674 std::map<G4long,G4GeomTestOverlapList>::iterator overlap =
>> 675 overlaps.begin();
>> 676 while( overlap != overlaps.end() )
>> 677 {
>> 678 logger->OverlappingDaughters( &(*overlap).second );
>> 679 ++overlap;
>> 680 }
>> 681 }
>> 682 }
>> 683
>> 684 //
>> 685 // ClearErrors
>> 686 //
>> 687 void G4GeomTestVolume::ClearErrors()
>> 688 {
>> 689 overshoots.clear();
>> 690 overlaps.clear();
219 } 691 }
220 692