source: trunk/source/geometry/solids/test/fred/src/FredTest3.cc@ 1350

Last change on this file since 1350 was 1316, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 19.3 KB
RevLine 
[1316]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// FredTest3
28//
29// Implementation of fred's third geometry tester
30//
31#include "FredTest3.hh"
32
33#include "Randomize.hh"
34#include "G4VSolid.hh"
35#include "G4UImanager.hh"
36#include "G4GeometryTolerance.hh"
37
38#include <time.h>
39#include <iomanip>
40#include <sstream>
41
42//
43// Constructor
44//
45FredTest3::FredTest3()
46{
47 //
48 // Defaults
49 //
50 SetDefaults();
51
52 //
53 // Zero error list
54 //
55 errorList = 0;
56}
57
58
59//
60// Destructor
61//
62FredTest3::~FredTest3()
63{
64 ClearErrors();
65}
66
67
68//
69// SetDefaults
70//
71// Set defaults values
72//
73void FredTest3::SetDefaults()
74{
75 target = G4ThreeVector( 0, 0, 0 );
76 widths = G4ThreeVector( 1*m, 1*m, 1*m );
77 grids = G4ThreeVector( 0, 0, 0 );
78
79 maxPoints = 10000;
80 maxErrors = 100;
81}
82
83
84//
85// GetRandomPoint
86//
87// Return a random point in three dimensions using the current
88// specs
89//
90G4ThreeVector FredTest3::GetRandomPoint() const {
91 G4double dx = widths.x()*GaussianRandom(10*m/widths.x()),
92 dy = widths.y()*GaussianRandom(10*m/widths.y()),
93 dz = widths.z()*GaussianRandom(10*m/widths.z());
94
95
96 if (grids.x() > 0) dx = grids.x()*G4int( dx/grids.x()+1 );
97 if (grids.y() > 0) dy = grids.y()*G4int( dx/grids.y()+1 );
98 if (grids.z() > 0) dz = grids.z()*G4int( dx/grids.z()+1 );
99
100 G4ThreeVector randvec( dx, dy, dz );
101
102 return target + randvec;
103}
104
105
106//
107// GaussianRandom
108//
109// Return Gaussian random number of unit width
110//
111// A classic, slow, but remarkably effective algorithm. Certainly good
112// enough for our purposes.
113//
114G4double FredTest3::GaussianRandom(const G4double cutoff) const {
115 if (cutoff <= 0) G4Exception( "Illegal cutoff" );
116
117 G4double answer;
118 do {
119 answer = -3.0;
120 for( G4int j = 0; j < 6; j++ ) answer += G4UniformRand();
121 answer *= std::sqrt(2.0);
122 } while( std::fabs(answer) > cutoff );
123
124 return(answer);
125}
126
127
128//
129// RunTest
130//
131// Do your stuff!
132//
133void FredTest3::RunTest( const G4VSolid *testVolume, std::ostream &logger )
134{
135 //
136 // Clear error list
137 //
138 ClearErrors();
139
140 //
141 // Output test parameters
142 //
143 time_t now;
144 time(&now);
145 G4String dateTime(ctime(&now));
146
147 logger << "% Fred test3 logged output " << dateTime;
148 logger << "% volume = " << testVolume->GetName() << G4endl;
149 logger << "% target = " << target << G4endl;
150 logger << "% widths = " << widths << G4endl;
151 logger << "% grids = " << grids << G4endl;
152 logger << "% maxPoints = " << maxPoints << G4endl;
153 logger << "% maxErrors = " << maxErrors << G4endl;
154
155 //
156 // Setup lists
157 //
158 FredTest3PointList inside(100);
159 FredTest3PointList outside(100);
160 FredTest3PointList surface(100);
161
162 //
163 // Set iostream precision to 14 digits
164 //
165 logger << std::setprecision(14);
166
167 //
168 // Set clock
169 //
170 // clock_t start = clock();
171
172 //
173 // Loop over points
174 //
175 G4int nIn = 0, nOut = 0, nSurf = 0;
176
177 G4int nPoint = 0;
178 G4int nError = 0;
179
180 for(;;) {
181
182 //
183 // Generate a random point
184 //
185 G4ThreeVector point = GetRandomPoint();
186
187 //
188 // Catagorize, test, and store
189 //
190 EInside catagory = testVolume->Inside( point );
191 switch( catagory ) {
192 case kOutside:
193 nOut++;
194 TestOutsidePoint( testVolume, &nError, &inside, point, logger );
195 outside.AddPoint( point );
196 break;
197
198 case kInside:
199 nIn++;
200 TestInsidePoint( testVolume, &nError, &outside, point, logger );
201 inside.AddPoint( point );
202 break;
203
204 case kSurface:
205 nSurf++;
206// TestSurfacePoint( testVolume, &nError, point, logger );
207 surface.AddPoint( point );
208 break;
209 }
210
211 //
212 // Return early if we have accumulated enough errors
213 //
214 if (nError >= maxErrors) {
215 logger << "% End of test (maximum number errors) ";
216 break;
217 }
218 if (++nPoint >= maxPoints) {
219 logger << "% End of test (maximum number points) ";
220 break;
221 }
222
223 }
224
225 time(&now);
226 G4String dateTime2(ctime(&now));
227 logger << dateTime2;
228
229 logger << "% Statistics: points=" << nPoint << " errors reported=" << nError << G4endl;
230
231 logger << "% inside=" << nIn << " outside=" << nOut
232 << " surface=" << nSurf << G4endl;
233 logger << "% cpu time=" << clock()/CLOCKS_PER_SEC << G4endl;
234
235 logger << "%(End of file)" << G4endl;
236}
237
238
239//
240// DebugError
241//
242// Recover previously logged error and setup particle gun appropriately
243//
244void FredTest3::RunDebug( const G4VSolid* solid, std::istream &logger )
245{
246
247 G4UImanager *UI = G4UImanager::GetUIpointer();
248 G4int errorIndex = 0;
249 while ( ! DebugError(solid, logger, ++errorIndex) ) {
250
251 UI->ApplyCommand( "/run/beamOn 1" );
252 }
253}
254
255//
256// DebugError
257//
258// Recover previously logged error and setup particle gun appropriately
259//
260G4int FredTest3::DebugError( const G4VSolid *, std::istream &logger,
261 const G4int errorIndex ) const
262{
263 G4ThreeVector p, v;
264
265 //
266 // Recover information from log file
267 //
268 G4int error = GetLoggedPV( logger, errorIndex, p, v );
269 if (error) return error;
270
271 G4cout << "DebugError " << errorIndex << ": p=" << p << ", v=" << v << G4endl;
272
273 //
274 // Setup fred to simulate this event
275 //
276 // We do this using the command line, a cheesy short-cut but
277 // rather useful in hacks like this.
278 //
279 // If you are writing your own serious GEANT4 application,
280 // please do something better.
281 //
282 G4UImanager *UI = G4UImanager::GetUIpointer();
283
284 UI->ApplyCommand( "/tracking/verbose 1" );
285
286 UI->ApplyCommand( "/fred/gun G4" );
287
288 UI->ApplyCommand( "/gun/particle geantino" );
289
290 std::ostringstream formatter1;
291 formatter1 << std::setprecision(14);
292 formatter1 << "/gun/position " << p.x() << " " << p.y() << " " << p.z() << " mm" << G4endl;
293 G4String commandBuffer=formatter1.str();
294 UI->ApplyCommand( commandBuffer );
295
296 std::ostringstream formatter2;
297 formatter2 << std::setprecision(14);
298 formatter2 << "/gun/direction " << v.x() << " " << v.y() << " " << v.z() << " mm" << G4endl;
299 commandBuffer=formatter2.str();
300 UI->ApplyCommand( commandBuffer );
301 return 0;
302}
303
304
305//
306// DebugInside
307//
308// Recover previously logged error and invoke G4VSolid::Inside
309//
310G4int FredTest3::DebugInside( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const
311{
312 G4ThreeVector p, v;
313
314 //
315 // Recover information from log file
316 //
317 G4int error = GetLoggedPV( logger, errorIndex, p, v );
318 if (error) return error;
319
320 //
321 // Call
322 //
323 G4cout << "testVolume->Inside(p): " << testVolume->Inside( p ) << G4endl;
324 return 0;
325}
326
327
328//
329// DebugToInP
330//
331// Recover previously logged error and invoke G4VSolid::DistanceToIn(p)
332//
333G4int FredTest3::DebugToInP( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const
334{
335 G4ThreeVector p, v;
336
337 //
338 // Recover information from log file
339 //
340 G4int error = GetLoggedPV( logger, errorIndex, p, v );
341 if (error) return error;
342
343 //
344 // Call
345 //
346 G4cout << "testVolume->DistanceToIn(p): " << testVolume->DistanceToIn( p ) << G4endl;
347 return 0;
348}
349
350
351//
352// DebugToInPV
353//
354// Recover previously logged error and invoke G4VSolid::DistanceToIn(p,v)
355//
356G4int FredTest3::DebugToInPV( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const
357{
358 G4ThreeVector p, v;
359
360 //
361 // Recover information from log file
362 //
363 G4int error = GetLoggedPV( logger, errorIndex, p, v );
364 if (error) return error;
365
366 //
367 // Call
368 //
369 G4double answer = testVolume->DistanceToIn( p, v );
370 G4cout << "testVolume->DistanceToIn(p,v): " << answer << G4endl;
371
372 p += answer*v;
373
374 G4cout << "testVolume->Inside(p+=answer*v):" << p << " " << testVolume->Inside(p) << G4endl;
375 return 0;
376}
377
378
379//
380// DebugToOutP
381//
382// Recover previously logged error and invoke G4VSolid::DistanceToOut(p)
383//
384G4int FredTest3::DebugToOutP( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const
385{
386 G4ThreeVector p, v;
387
388 //
389 // Recover information from log file
390 //
391 G4int error = GetLoggedPV( logger, errorIndex, p, v );
392 if (error) return error;
393
394 //
395 // Call
396 //
397 G4cout << "testVolume->DistanceToOut(p): " << testVolume->DistanceToOut( p ) << G4endl;
398
399 return 0;
400}
401
402
403//
404// DebugToOutPV
405//
406// Recover previously logged error and invoke G4VSolid::DistanceToOut(p,v)
407//
408G4int FredTest3::DebugToOutPV( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const
409{
410 G4ThreeVector p, v;
411
412 //
413 // Recover information from log file
414 //
415 G4int error = GetLoggedPV( logger, errorIndex, p, v );
416 if (error) return error;
417
418 //
419 // Call
420 //
421 G4bool validNorm;
422 G4ThreeVector norm;
423 G4double answer = testVolume->DistanceToOut( p, v, true, &validNorm, &norm);
424 G4cout << "testVolume->DistanceToOut( p, v ): " << answer << " validNorm: " << validNorm << G4endl;
425
426 p += answer*v;
427 G4cout << "testVolume->Inside(p += answer*v): " << testVolume->Inside(p) << G4endl;
428
429 return 0;
430}
431
432
433//
434// --------------------------------------
435// Tests
436//
437
438//
439// TestOutsidePoint
440//
441void FredTest3::TestOutsidePoint( const G4VSolid *testVolume, G4int *nError,
442 const FredTest3PointList *inside, const G4ThreeVector point, std::ostream &logger )
443{
444 // Check if point is really outside
445 G4double safeDistance = testVolume->DistanceToIn( point );
446 if (safeDistance <= 0.0) {
447 ReportError( nError, point, 0, "T0: DistanceToIn(p) <= 0", logger );
448 return;
449 }
450
451 // Loop over inside points
452 for( G4int i=0; i < inside->NumPoints(); i++ ) {
453
454 // Direction vector from tested point to the inside point
455 G4ThreeVector vr = (*inside)[i] - point;
456 G4ThreeVector v = vr.unit();
457
458 // Distance from the outside point to solid in direction
459 // to the tested inside point
460 G4double dist = testVolume->DistanceToIn( point, v );
461
462 if (dist <= 0) {
463 ReportError( nError, point, v, "T0: DistanceToIn(p,v) <= 0", logger );
464 continue;
465 }
466 if (dist == kInfinity) {
467 ReportError( nError, point, v, "T0: DistanceToIn(p,v) == kInfinity", logger );
468 continue;
469 }
470 if (dist < safeDistance+1E-10) {
471 ReportError( nError, point, v, "T0: DistanceToIn(p,v) < DistanceToIn(p)", logger );
472 continue;
473 }
474
475 // Tested point moved to the solid surface
476 // (in direction to the tested inside point):
477 G4ThreeVector p = point + dist*v;
478
479 // Inside(p) should return kSurface.
480 // If kOutside is returned, report undershoots; if kInside, report overshoots
481 EInside insideOrNot = testVolume->Inside( p );
482 if (insideOrNot == kOutside) {
483 ReportError( nError, point, v, "T0: DistanceToIn(p,v) undershoots", logger );
484 continue;
485 }
486 if (insideOrNot == kInside) {
487 ReportError( nError, point, v, "TO: DistanceToIn(p,v) overshoots", logger );
488 continue;
489 }
490
491 // DistanceToIn for the computed point on the surface
492 dist = testVolume->DistanceToIn( p );
493 //if (dist != 0) {
494 if (dist > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()) {
495 ReportError( nError, p, v, "T02: DistanceToIn(p) should be zero", logger );
496 continue;
497 }
498
499 // DistanceToOut for the computed point on the surface
500 dist = testVolume->DistanceToOut( p );
501 //if (dist != 0) {
502 if (dist > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()) {
503 ReportError( nError, p, v, "T02: DistanceToOut(p) should be zero", logger );
504 continue;
505 }
506
507 // DistanceToIn with direction for the computed point on the surface
508 dist = testVolume->DistanceToIn( p, v );
509 //if (dist != 0) {
510 if (dist > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()) {
511 ReportError( nError, p, v, "T02: DistanceToIn(p,v) should be zero", logger );
512 continue;
513 }
514
515 // DistanceToOut with direction for the computed point on the surface
516 G4bool validNorm;
517 G4ThreeVector norm;
518
519 dist = testVolume->DistanceToOut( p, v, true, &validNorm, &norm );
520 // if (dist == 0) continue;
521 // if (dist < G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()) continue;
522 // Why quit here withour Error ???
523
524 if (dist == kInfinity) {
525 ReportError( nError, p, v, "T02: DistanceToOut(p,v) == kInfinity", logger );
526 continue;
527 }
528
529 if (validNorm) {
530 if (norm.dot(v) < 0) {
531 ReportError( nError, p, v, "T02: Outgoing normal incorrect", logger );
532 continue;
533 }
534 }
535
536 // The point on surface (entering point) moved to the second point
537 // on the surface (exiting point)
538 // (in direction to the tested inside point):
539 G4ThreeVector p2 = p + v*dist;
540
541 insideOrNot = testVolume->Inside(p2);
542 if (insideOrNot == kInside) {
543 ReportError( nError, p, v, "T02: DistanceToOut(p,v) undershoots", logger );
544 continue;
545 }
546 if (insideOrNot == kOutside) {
547 ReportError( nError, p, v, "TO2: DistanceToOut(p,v) overshoots", logger );
548 continue;
549 }
550
551 // DistanceToIn from the exiting point (in the exiting direction)
552 // If infinity, there is no intersectin with solid anymore, what means that the solid
553 // lied entirely between entering and existing points, and the validNorm should have been true.
554
555 dist = testVolume->DistanceToIn(p2,v);
556 if (validNorm) {
557 if (dist != kInfinity) {
558 ReportError( nError, p, v, "TO2: DistanceToOut incorrectly returns validNorm==true", logger );
559 continue;
560 }
561 }
562 else {
563 if (dist == kInfinity) {
564 ReportError( nError, p, v, "TO2: DistanceToOut incorrectly returns validNorm==false", logger );
565 continue;
566 }
567 }
568 }
569}
570
571
572
573//
574// TestInsidePoint
575//
576void FredTest3::TestInsidePoint( const G4VSolid *testVolume, G4int *nError,
577 const FredTest3PointList *outside, const G4ThreeVector point, std::ostream &logger )
578{
579 // Check if point is really inside
580 G4double safeDistance = testVolume->DistanceToOut( point );
581 if (safeDistance <= 0.0) {
582 ReportError( nError, point, 0, "TI: DistanceToOut(p) <= 0", logger );
583 return;
584 }
585
586 // Loop over outside points
587 for( G4int i=0; i < outside->NumPoints(); i++ ) {
588
589 // Direction vector from tested point to the outside point
590 G4ThreeVector vr = (*outside)[i] - point;
591 G4ThreeVector v = vr.unit();
592
593 G4bool validNorm;
594 G4ThreeVector norm;
595
596 // Distance from the inside point to solid border in direction
597 // to the tested outside point
598 G4double dist = testVolume->DistanceToOut( point, v, true, &validNorm, &norm );
599 if (dist <= 0) {
600 ReportError( nError, point, v, "TI: DistanceToOut(p,v) <= 0", logger );
601 continue;
602 }
603 if (dist == kInfinity) {
604 ReportError( nError, point, v, "TI: DistanceToOut(p,v) == kInfinity", logger );
605 continue;
606 }
607 if (dist < safeDistance+1E-10) {
608 ReportError( nError, point, v, "TI: DistanceToOut(p,v) < DistanceToOut(p)", logger );
609 continue;
610 }
611
612 if (validNorm) {
613 if (norm.dot(v) < 0) {
614 ReportError( nError, point, v, "TI: Outgoing normal incorrect", logger );
615 continue;
616 }
617 }
618
619 // Tested point moved to the solid surface
620 // (in direction to the tested outside point):
621 G4ThreeVector p = point + v*dist;
622
623 EInside insideOrNot = testVolume->Inside(p);
624 if (insideOrNot == kInside) {
625 ReportError( nError, point, v, "TI: DistanceToOut(p,v) undershoots", logger );
626 continue;
627 }
628 if (insideOrNot == kOutside) {
629 ReportError( nError, point, v, "TI: DistanceToOut(p,v) overshoots", logger );
630 continue;
631 }
632 }
633}
634
635
636//
637// ReportError
638//
639// Report the specified error message, but only if it has not been reported a zillion
640// times already.
641//
642void FredTest3::ReportError( G4int *nError, const G4ThreeVector p,
643 const G4ThreeVector v, const G4String comment, std::ostream &logger )
644{
645 //
646 // Have we encountered this error message before?
647 //
648 FredTest3ErrorList *last=0, *errors = errorList;
649 while( errors ) {
650 //
651 // Note: below is an expensive comparison. This could be replaced with something
652 // faster if we really wanted, like a hash. But is it worth the trouble? Nah.
653 //
654 if (errors->message == comment) {
655 //
656 // Yup. Increment count, and return now if the count is too high
657 //
658 if ( ++errors->nUsed > 5 ) return;
659 break;
660 }
661 last = errors;
662 errors = errors->next;
663 }
664
665 if (errors == 0) {
666 //
667 // New error: add it the end of our list
668 //
669 errors = new FredTest3ErrorList;
670 errors->message = comment;
671 errors->nUsed = 1;
672 errors->next = 0;
673 if (errorList)
674 last->next = errors;
675 else
676 errorList = errors;
677 }
678
679 //
680 // Output the message
681 //
682 logger << "% " << comment;
683 if (errors->nUsed == 5) logger << " (any further such errors suppressed)";
684 logger << G4endl;
685
686 logger << ++(*nError) << " " << p.x() << " " << p.y() << " " << p.z()
687 << " " << v.x() << " " << v.y() << " " << v.z() << G4endl;
688}
689
690
691//
692// ClearErrors
693//
694// Reset list of errors (and clear memory)
695//
696void FredTest3::ClearErrors()
697{
698 FredTest3ErrorList *here, *next;
699
700 here = errorList;
701 while( here ) {
702 next = here->next;
703 delete here;
704 here = next;
705 }
706}
707
708
709//
710// GetLoggedPV
711//
712// Get the p and v vectors stored in a test3 log file
713//
714G4int FredTest3::GetLoggedPV( std::istream &logger, const G4int errorIndex,
715 G4ThreeVector &p, G4ThreeVector &v ) const
716{
717 logger >> std::setprecision(14); // I wonder if this is necessary?
718 logger.seekg(0);
719
720 //
721 // Search for the requested error index, skipping comments along the way
722 //
723 for(;;) {
724 while( logger.peek() == '%' ) logger.ignore( 99999, '\n' );
725
726 if (logger.peek() == EOF) return 1;
727
728 G4int thisIndex;
729
730 logger >> thisIndex;
731
732 if (thisIndex == errorIndex) break;
733
734 logger.ignore( 99999, '\n' );
735 }
736
737 //
738 // Extract the vectors
739 //
740 // We should probably have some error checking here...
741 //
742 G4double x, y, z;
743
744 logger >> x >> y >> z;
745 p = G4ThreeVector(x*mm,y*mm,z*mm);
746
747 logger >> x >> y >> z;
748 v = G4ThreeVector(x*mm,y*mm,z*mm);
749
750 return 0;
751}
752
753
754//
755// --------------------------------------
756// FredTest3PointList stuff
757//
758
759FredTest3PointList::FredTest3PointList( G4int size )
760{
761 pointList = new G4ThreeVector[size];
762 maxPoints = size;
763 numPoints = 0;
764}
765
766FredTest3PointList::~FredTest3PointList()
767{
768 delete [] pointList;
769}
770
771void FredTest3PointList::AddPoint( G4ThreeVector newPoint )
772{
773 if (numPoints < maxPoints) {
774 //
775 // Since we still have space, append the point on the
776 // end of the list
777 //
778 pointList[numPoints] = newPoint;
779 numPoints++;
780 }
781 else {
782 //
783 // Our list is filled up, so replace a random
784 // entry
785 //
786 G4int irand = G4int(G4UniformRand()*( (G4double)maxPoints ));
787 pointList[irand] = newPoint;
788 }
789}
Note: See TracBrowser for help on using the repository browser.