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

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// FredTest3
//
// Implementation of fred's third geometry tester
//
#include "FredTest3.hh"

#include "Randomize.hh"
#include "G4VSolid.hh"
#include "G4UImanager.hh"
#include "G4GeometryTolerance.hh"

#include <time.h>
#include <iomanip>
#include <sstream>

//
// Constructor
//
FredTest3::FredTest3()
{
        //
        // Defaults
        //
        SetDefaults();
        
        //
        // Zero error list
        //
        errorList = 0;
}


//
// Destructor
//
FredTest3::~FredTest3()
{
        ClearErrors();
}


//
// SetDefaults
//
// Set defaults values
//
void FredTest3::SetDefaults()
{
        target = G4ThreeVector( 0, 0, 0 );
        widths = G4ThreeVector( 1*m, 1*m, 1*m );
        grids  = G4ThreeVector( 0, 0, 0 );
        
        maxPoints = 10000;
        maxErrors = 100;
}


//
// GetRandomPoint
//
// Return a random point in three dimensions using the current
// specs 
//
G4ThreeVector FredTest3::GetRandomPoint() const {
        G4double dx = widths.x()*GaussianRandom(10*m/widths.x()),
                 dy = widths.y()*GaussianRandom(10*m/widths.y()),
                 dz = widths.z()*GaussianRandom(10*m/widths.z());
                 
                 
        if (grids.x() > 0) dx = grids.x()*G4int( dx/grids.x()+1 );
        if (grids.y() > 0) dy = grids.y()*G4int( dx/grids.y()+1 );
        if (grids.z() > 0) dz = grids.z()*G4int( dx/grids.z()+1 );

        G4ThreeVector randvec( dx, dy, dz );
        
        return target + randvec;
}


//
// GaussianRandom
//
// Return Gaussian random number of unit width
//
// A classic, slow, but remarkably effective algorithm. Certainly good
// enough for our purposes.
//
G4double FredTest3::GaussianRandom(const G4double cutoff) const {
        if (cutoff <= 0) G4Exception( "Illegal cutoff" );

        G4double answer;
        do {
                answer = -3.0;
                for( G4int j = 0; j < 6; j++ ) answer += G4UniformRand();
                answer *= std::sqrt(2.0);
        } while( std::fabs(answer) > cutoff );
        
        return(answer);
}


//
// RunTest
//
// Do your stuff!
//
void FredTest3::RunTest( const G4VSolid *testVolume, std::ostream &logger )
{
        //
        // Clear error list
        //
        ClearErrors();
        
        //
        // Output test parameters
        //
        time_t now;
        time(&now);
        G4String dateTime(ctime(&now));
        
        logger << "% Fred test3 logged output " << dateTime;
        logger << "% volume =    " << testVolume->GetName() << G4endl;
        logger << "% target =    " << target << G4endl;
        logger << "% widths =    " << widths << G4endl;
        logger << "% grids  =    " << grids  << G4endl;
        logger << "% maxPoints = " << maxPoints << G4endl;
        logger << "% maxErrors = " << maxErrors << G4endl;

        //
        // Setup lists
        //
        FredTest3PointList inside(100);
        FredTest3PointList outside(100);
        FredTest3PointList surface(100);
        
        //
        // Set iostream precision to 14 digits
        //
        logger << std::setprecision(14);
        
        //
        // Set clock
        //
        // clock_t start = clock();
        
        //
        // Loop over points
        //
        G4int nIn = 0, nOut = 0, nSurf = 0;
        
        G4int nPoint = 0;
        G4int nError = 0;
        
        for(;;) {
  
                //
                // Generate a random point
                //
                G4ThreeVector   point = GetRandomPoint();
                
                //
                // Catagorize, test, and store
                //
                EInside catagory = testVolume->Inside( point );
                switch( catagory ) {
                        case kOutside:
                        nOut++;
                        TestOutsidePoint( testVolume, &nError, &inside, point, logger );
                        outside.AddPoint( point );
                        break;
                        
                        case kInside:
                        nIn++;
                        TestInsidePoint( testVolume, &nError, &outside, point, logger );
                        inside.AddPoint( point );
                        break;
                        
                        case kSurface:
                        nSurf++;
//                      TestSurfacePoint( testVolume, &nError, point, logger );
                        surface.AddPoint( point );
                        break;
                }
                
                //
                // Return early if we have accumulated enough errors
                //
                if (nError >= maxErrors) {
                        logger << "% End of test (maximum number errors) ";
                        break;
                }
                if (++nPoint >= maxPoints) {
                        logger << "% End of test (maximum number points) ";
                        break;
                }
        
        }

        time(&now);
        G4String dateTime2(ctime(&now));
        logger << dateTime2;

        logger << "% Statistics: points=" << nPoint << " errors reported=" << nError << G4endl;

        logger << "%             inside=" << nIn << " outside=" << nOut
               << " surface=" << nSurf << G4endl;
        logger << "%             cpu time=" << clock()/CLOCKS_PER_SEC << G4endl;
               
        logger << "%(End of file)" << G4endl;
}


//
// DebugError
//
// Recover previously logged error and setup particle gun appropriately
//
void FredTest3::RunDebug( const G4VSolid* solid, std::istream &logger )
{

        G4UImanager *UI = G4UImanager::GetUIpointer();
        G4int errorIndex = 0;
        while ( ! DebugError(solid, logger, ++errorIndex) ) {

                UI->ApplyCommand( "/run/beamOn 1" );
        }        
}

//
// DebugError
//
// Recover previously logged error and setup particle gun appropriately
//
G4int FredTest3::DebugError( const G4VSolid *, std::istream &logger,
                             const G4int errorIndex ) const
{
        G4ThreeVector p, v;
        
        //
        // Recover information from log file
        //
        G4int error = GetLoggedPV( logger, errorIndex, p, v );
        if (error) return error;
        
        G4cout << "DebugError  " << errorIndex << ":  p=" << p << ",  v=" << v << G4endl; 
                        
        //
        // Setup fred to simulate this event
        //
        // We do this using the command line, a cheesy short-cut but
        // rather useful in hacks like this.
        //
        // If you are writing your own serious GEANT4 application,
        // please do something better.
        //
        G4UImanager *UI = G4UImanager::GetUIpointer();

        UI->ApplyCommand( "/tracking/verbose 1" );

        UI->ApplyCommand( "/fred/gun G4" );
        
        UI->ApplyCommand( "/gun/particle geantino" );
        
        std::ostringstream formatter1;
        formatter1 << std::setprecision(14);
        formatter1 << "/gun/position "  << p.x() << " " << p.y() << " " << p.z() << " mm" << G4endl;
        G4String commandBuffer=formatter1.str();
        UI->ApplyCommand( commandBuffer );
        
        std::ostringstream formatter2;
        formatter2 << std::setprecision(14);
        formatter2 << "/gun/direction " << v.x() << " " << v.y() << " " << v.z() << " mm" << G4endl;
        commandBuffer=formatter2.str();
        UI->ApplyCommand( commandBuffer );
        return 0;
}


//
// DebugInside
//
// Recover previously logged error and invoke G4VSolid::Inside
//
G4int FredTest3::DebugInside( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const
{       
        G4ThreeVector p, v;
        
        //
        // Recover information from log file
        //
        G4int error = GetLoggedPV( logger, errorIndex, p, v );
        if (error) return error;
        
        //
        // Call
        //
        G4cout << "testVolume->Inside(p): " << testVolume->Inside( p ) << G4endl; 
        return 0;
}


//
// DebugToInP
//
// Recover previously logged error and invoke G4VSolid::DistanceToIn(p)
//
G4int FredTest3::DebugToInP( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const
{       
        G4ThreeVector p, v;
        
        //
        // Recover information from log file
        //
        G4int error = GetLoggedPV( logger, errorIndex, p, v );
        if (error) return error;
        
        //
        // Call
        //
        G4cout << "testVolume->DistanceToIn(p): " <<  testVolume->DistanceToIn( p )  << G4endl; 
        return 0;
}


//
// DebugToInPV
//
// Recover previously logged error and invoke G4VSolid::DistanceToIn(p,v)
//
G4int FredTest3::DebugToInPV( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const
{       
        G4ThreeVector p, v;
        
        //
        // Recover information from log file
        //
        G4int error = GetLoggedPV( logger, errorIndex, p, v );
        if (error) return error;
        
        //
        // Call
        //
        G4double answer = testVolume->DistanceToIn( p, v );
        G4cout << "testVolume->DistanceToIn(p,v): " << answer << G4endl; 
        
        p += answer*v;
        
        G4cout << "testVolume->Inside(p+=answer*v):" <<  p << " " << testVolume->Inside(p) << G4endl;
        return 0;
}


//
// DebugToOutP
//
// Recover previously logged error and invoke G4VSolid::DistanceToOut(p)
//
G4int FredTest3::DebugToOutP( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const
{       
        G4ThreeVector p, v;
        
        //
        // Recover information from log file
        //
        G4int error = GetLoggedPV( logger, errorIndex, p, v );
        if (error) return error;
        
        //
        // Call
        //
        G4cout << "testVolume->DistanceToOut(p): " << testVolume->DistanceToOut( p ) << G4endl; 

        return 0;
}


//
// DebugToOutPV
//
// Recover previously logged error and invoke G4VSolid::DistanceToOut(p,v)
//
G4int FredTest3::DebugToOutPV( const G4VSolid *testVolume, std::istream &logger, const G4int errorIndex ) const
{       
        G4ThreeVector p, v;
        
        //
        // Recover information from log file
        //
        G4int error = GetLoggedPV( logger, errorIndex, p, v );
        if (error) return error;

        //
        // Call
        //
        G4bool validNorm;
        G4ThreeVector norm;
        G4double answer = testVolume->DistanceToOut( p, v, true, &validNorm, &norm);
        G4cout << "testVolume->DistanceToOut( p, v ): " << answer << " validNorm: " << validNorm << G4endl; 
        
        p += answer*v;
        G4cout << "testVolume->Inside(p += answer*v): " << testVolume->Inside(p) << G4endl; 
        
        return 0;
}


//
// --------------------------------------
// Tests
//

//
// TestOutsidePoint
//
void FredTest3::TestOutsidePoint( const G4VSolid *testVolume, G4int *nError,
                                  const FredTest3PointList *inside, const G4ThreeVector point, std::ostream &logger )
{
        // Check if point is really outside
        G4double safeDistance = testVolume->DistanceToIn( point );
        if (safeDistance <= 0.0) {
                ReportError( nError, point, 0, "T0: DistanceToIn(p) <= 0", logger );
                return;
        }
        
        // Loop over inside points
        for( G4int i=0; i < inside->NumPoints(); i++ ) {
  
                // Direction vector from tested point to the inside point
                G4ThreeVector vr = (*inside)[i] - point;
                G4ThreeVector v = vr.unit();
                
                // Distance from the outside point to solid in direction
                // to the tested inside point
                G4double dist = testVolume->DistanceToIn( point, v );

                if (dist <= 0) {
                        ReportError( nError, point, v, "T0: DistanceToIn(p,v) <= 0", logger );
                        continue;
                }
                if (dist == kInfinity) {
                        ReportError( nError, point, v, "T0: DistanceToIn(p,v) == kInfinity", logger );
                        continue;
                }
                if (dist < safeDistance+1E-10) {
                        ReportError( nError, point, v, "T0: DistanceToIn(p,v) < DistanceToIn(p)", logger );
                        continue;
                }
                
                // Tested point moved to the solid surface
                // (in direction to the tested inside point): 
                G4ThreeVector p = point + dist*v;
                
                // Inside(p) should return kSurface.
                // If kOutside is returned, report undershoots; if kInside, report overshoots
                EInside insideOrNot = testVolume->Inside( p );
                if (insideOrNot == kOutside) {
                        ReportError( nError, point, v, "T0: DistanceToIn(p,v) undershoots", logger );
                        continue;
                }
                if (insideOrNot == kInside) {
                        ReportError( nError, point, v, "TO: DistanceToIn(p,v) overshoots", logger );
                        continue;
                }
                
                // DistanceToIn for the computed point on the surface
                dist = testVolume->DistanceToIn( p );
                //if (dist != 0) {
                if (dist > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()) {
                        ReportError( nError, p, v, "T02: DistanceToIn(p) should be zero", logger );
                        continue;
                }
                
                // DistanceToOut for the computed point on the surface
                dist = testVolume->DistanceToOut( p );
                //if (dist != 0) {
                if (dist > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()) {
                        ReportError( nError, p, v, "T02: DistanceToOut(p) should be zero", logger );
                        continue;
                }
                
                // DistanceToIn with direction for the computed point on the surface
                dist = testVolume->DistanceToIn( p, v );
                //if (dist != 0) {
                if (dist > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()) {
                        ReportError( nError, p, v, "T02: DistanceToIn(p,v) should be zero", logger );
                        continue;
                }       
                
                // DistanceToOut with direction for the computed point on the surface
                G4bool validNorm;
                G4ThreeVector norm;
                
                dist = testVolume->DistanceToOut( p, v, true, &validNorm, &norm );
                // if (dist == 0) continue;
                // if (dist < G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()) continue;
                // Why quit here withour Error ??? 
    
                if (dist == kInfinity) {
                        ReportError( nError, p, v, "T02: DistanceToOut(p,v) == kInfinity", logger );
                        continue;
                }
                
                if (validNorm) {
                        if (norm.dot(v) < 0) {
                                ReportError( nError, p, v, "T02: Outgoing normal incorrect", logger );
                                continue;
                        }
                }
                
                // The point on surface (entering point) moved to the second point
                // on the surface (exiting point)
                // (in direction to the tested inside point):
                G4ThreeVector p2 = p + v*dist;
                
                insideOrNot = testVolume->Inside(p2);
                if (insideOrNot == kInside) {
                        ReportError( nError, p, v, "T02: DistanceToOut(p,v) undershoots", logger );
                        continue;
                }
                if (insideOrNot == kOutside) {
                        ReportError( nError, p, v, "TO2: DistanceToOut(p,v) overshoots", logger );
                        continue;
                }
                        
                // DistanceToIn from the exiting point (in the exiting direction)
                // If infinity, there is no intersectin with solid anymore, what means that the solid
                // lied entirely between entering and existing points, and the validNorm should have been true.
                
                dist = testVolume->DistanceToIn(p2,v);
                if (validNorm) {
                        if (dist != kInfinity) {
                                ReportError( nError, p, v, "TO2: DistanceToOut incorrectly returns validNorm==true", logger );
                                continue;
                        }
                }
                else {
                        if (dist == kInfinity) {
                                ReportError( nError, p, v, "TO2: DistanceToOut incorrectly returns validNorm==false", logger );
                                continue;
                        }
                }
        }
}



//
// TestInsidePoint
//
void FredTest3::TestInsidePoint( const G4VSolid *testVolume, G4int *nError,
                                 const FredTest3PointList *outside, const G4ThreeVector point, std::ostream &logger )
{
        // Check if point is really inside
        G4double safeDistance = testVolume->DistanceToOut( point );
        if (safeDistance <= 0.0) {
                ReportError( nError, point, 0, "TI: DistanceToOut(p) <= 0", logger );
                return;
        }
        
        // Loop over outside points
        for( G4int i=0; i < outside->NumPoints(); i++ ) {

                // Direction vector from tested point to the outside point
                G4ThreeVector vr = (*outside)[i] - point;
                G4ThreeVector v = vr.unit();

                G4bool validNorm;
                G4ThreeVector norm;
                
                // Distance from the inside point to solid border in direction
                // to the tested outside point
                G4double dist = testVolume->DistanceToOut( point, v, true, &validNorm, &norm );
                if (dist <= 0) {
                        ReportError( nError, point, v, "TI: DistanceToOut(p,v) <= 0", logger );
                        continue;
                }
                if (dist == kInfinity) {
                        ReportError( nError, point, v, "TI: DistanceToOut(p,v) == kInfinity", logger );
                        continue;
                }
                if (dist < safeDistance+1E-10) {
                        ReportError( nError, point, v, "TI: DistanceToOut(p,v) < DistanceToOut(p)", logger );
                        continue;
                }

                if (validNorm) {
                        if (norm.dot(v) < 0) {
                                ReportError( nError, point, v, "TI: Outgoing normal incorrect", logger );
                                continue;
                        }
                }
                
                // Tested point moved to the solid surface
                // (in direction to the tested outside point): 
                G4ThreeVector p = point + v*dist;
                
                EInside insideOrNot = testVolume->Inside(p);
                if (insideOrNot == kInside) {
                        ReportError( nError, point, v, "TI: DistanceToOut(p,v) undershoots", logger );
                        continue;
                }
                if (insideOrNot == kOutside) {
                        ReportError( nError, point, v, "TI: DistanceToOut(p,v) overshoots", logger );
                        continue;
                }
        }
}


//
// ReportError
//
// Report the specified error message, but only if it has not been reported a zillion
// times already.
//
void FredTest3::ReportError( G4int *nError, const G4ThreeVector p, 
                             const G4ThreeVector v, const G4String comment, std::ostream &logger )
{
        //
        // Have we encountered this error message before?
        //
        FredTest3ErrorList *last=0, *errors = errorList;
        while( errors ) {
                //
                // Note: below is an expensive comparison. This could be replaced with something
                // faster if we really wanted, like a hash. But is it worth the trouble? Nah.
                //
                if (errors->message == comment) {
                        //
                        // Yup. Increment count, and return now if the count is too high
                        //
                        if ( ++errors->nUsed > 5 ) return;
                        break;
                }
                last = errors;
                errors = errors->next;
        }
        
        if (errors == 0) {
                //
                // New error: add it the end of our list
                //
                errors = new FredTest3ErrorList;
                errors->message = comment;
                errors->nUsed = 1;
                errors->next = 0;
                if (errorList) 
                        last->next = errors;
                else
                        errorList = errors;
        }

        //
        // Output the message
        //      
        logger << "% " << comment;
        if (errors->nUsed == 5) logger << " (any further such errors suppressed)";
        logger << G4endl;
        
        logger << ++(*nError) << " " << p.x() << " " << p.y() << " " << p.z() 
                              << " " << v.x() << " " << v.y() << " " << v.z() << G4endl;
}


//
// ClearErrors
//
// Reset list of errors (and clear memory)
//
void FredTest3::ClearErrors()
{
        FredTest3ErrorList *here, *next;
        
        here = errorList;
        while( here ) {
                next = here->next;
                delete here;
                here = next;
        }
}


//
// GetLoggedPV
//
// Get the p and v vectors stored in a test3 log file
//
G4int FredTest3::GetLoggedPV( std::istream &logger, const G4int errorIndex,
                              G4ThreeVector &p, G4ThreeVector &v        ) const
{
        logger >> std::setprecision(14);                // I wonder if this is necessary?
        logger.seekg(0);

        //
        // Search for the requested error index, skipping comments along the way
        //
        for(;;) {
                while( logger.peek() == '%' ) logger.ignore( 99999, '\n' );
                
                if (logger.peek() == EOF) return 1;
                
                G4int thisIndex;
                
                logger >> thisIndex;
                
                if (thisIndex == errorIndex) break;
                
                logger.ignore( 99999, '\n' );
        }
        
        //
        // Extract the vectors
        //
        // We should probably have some error checking here...
        //
        G4double x, y, z;
        
        logger >> x >> y >> z;
        p = G4ThreeVector(x*mm,y*mm,z*mm);
        
        logger >> x >> y >> z;
        v = G4ThreeVector(x*mm,y*mm,z*mm);

        return 0;
}


//
// --------------------------------------
// FredTest3PointList stuff
//

FredTest3PointList::FredTest3PointList( G4int size )
{
        pointList = new G4ThreeVector[size];
        maxPoints = size;
        numPoints = 0;
}

FredTest3PointList::~FredTest3PointList()
{
        delete [] pointList;
}

void FredTest3PointList::AddPoint( G4ThreeVector newPoint )
{
        if (numPoints < maxPoints) {
                //
                // Since we still have space, append the point on the
                // end of the list
                //
                pointList[numPoints] = newPoint;
                numPoints++;
        }
        else {
                //
                // Our list is filled up, so replace a random
                // entry
                //
                G4int irand = G4int(G4UniformRand()*( (G4double)maxPoints ));
                pointList[irand] = newPoint;
        }
}