// PosGen_Paint.cc
// Contact person:Matthew Mottram <m.mottram@qmul.ac.uk>
// See PosGen_Paint.hh for more details
//———————————————————————//
#include <G4VPhysicalVolume.hh>
#include <G4LogicalVolume.hh>
#include <G4VSolid.hh>
#include <G4TransportationManager.hh>
#include <G4Navigator.hh>
#include <G4VoxelLimits.hh>
#include <G4AffineTransform.hh>
#include <Randomize.hh>
#include <G4GeometryTolerance.hh>

#include <RAT/GLG4StringUtil.hh>

#include <sstream>

#include <RAT/PosGen_Paint.hh>
#include <RAT/Log.hh>

using std::vector;

namespace RAT {
void PosGen_Paint::SetState( G4String newValues )
{
    newValues = trim(newValues);
    if(newValues.length() == 0) // Empty arg, print help
    {
        info   << "Current state of this GLG4PosGen_Paint:\n\""
               << GetState() << "\"\n";
        detail << "Format of argument to GLG4PosGen_Paint::SetState:" << endl
               << " \n\"x y z t(+/-)\" with all inputs in [mm], a point in a volume,"
                  " or \"name t\", the name of a volume. "
                  "t is the thickness of the surface.\n"
                  "An optional + after t indicates to generate events only "
                  "on the outer surface, or an optional - indicates to "
                  "generate only on the outer surface. "
                  "With neither, generate on both." << endl;
    } else {
        vector<std::string> input = split(newValues, " ");
        Log::Assert(input.size()>=2,
                "PosGen_Paint::SetState: no thickness specified");
        fStateThicknesses.push_back(input.back());
        input.pop_back();
        fStates.push_back(join(input, " "));
        fListValid = false;
    }
}

G4String PosGen_Paint::GetState() const
{
    std::stringstream result;
    for(unsigned int i=0; i < fStates.size(); i++)
        result << fStates.at(i) << " " << fStateThicknesses.at(i) << ", ";
    return result.str();
}

void PosGen_Paint::GeneratePosition(G4ThreeVector& argResult)
{
    // Check if must regenerate list
    if(!fListValid)
    {
        fVolumes.Clear();
        fTotalSurfaceArea = 0.0;
        for(unsigned int i = 0; i < fStates.size(); i++)
        {
            unsigned int j = fVolumes.Size();
            std::string thickness = fStateThicknesses.at(i);
            // Decide which surfaces to add based on last character of
            // thickness string.
            std::string::iterator lastChar = thickness.end();
            --lastChar;
            if(*lastChar == '+'){//Outer only
                fVolumes.AddString(fStates.at(i), VolumeList::SELF);
            } else if(*lastChar == '-') {//Inner only
                fVolumes.AddString(fStates.at(i), VolumeList::DAUGHTERS);
            } else {//Both
                fVolumes.AddString(fStates.at(i), VolumeList::BOTH);
            }
            // For each added volume add the associated thickness
            for(; j < fVolumes.Size(); j++)
                fThicknesses.push_back(to_double(thickness));
        }
        for(unsigned int i = 0; i < fVolumes.Size(); i++)
            fTotalSurfaceArea += ComputeArea(fVolumes.GetVolume(i));
        fListValid = true;
    }
    // Now can generate a point, first choose the volume
    const double trialArea = G4UniformRand() * fTotalSurfaceArea;
    debug << "Trial area is " << trialArea
          << " and Total area is " << fTotalSurfaceArea << newline;
    int iVolume = -1; //iVolume incremented each time through, including first
    for(double cumulativeArea = 0; trialArea > cumulativeArea;){
        ++iVolume;
        cumulativeArea += ComputeArea(fVolumes.GetVolume(iVolume));
    }
    // Extract the solid
    G4VSolid* solid = fVolumes.GetVolume(iVolume)->GetLogicalVolume()->GetSolid();

    // Note: It is possible to ask the solid for a point on the surface,
    //       and thus save the computation below.
    //       However this seems to not weight the choices by surface area.
    //       For example a boolean solid is merely a 50/50 chance
    //       between two solids regardless of surface area.
    //       Thus the below routine is used.
    //       P Jones
    debug << "Selecting volume named " << fVolumes.GetVolume(iVolume)->GetName()
          << "with index " << iVolume << newline;
    const double thickness = fThicknesses[iVolume];
    // Now calculate the bounding box
    double x0, y0, z0, x1, y1, z1;
    {
        G4VoxelLimits voxelLimits;  // Defaults to "infinite" limits.
        G4AffineTransform affineTransform; // Defaults to no transform
        solid->CalculateExtent( kXAxis, voxelLimits, affineTransform, x0, x1 );
        solid->CalculateExtent( kYAxis, voxelLimits, affineTransform, y0, y1 );
        solid->CalculateExtent( kZAxis, voxelLimits, affineTransform, z0, z1 );
    }
    const double surfaceTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
    const double sphereRadius = 0.50001 * sqrt((x1 - x0)*(x1 - x0) + (y1 - y0)*(y1 - y0) + (z1 - z0)*(z1 - z0)) + surfaceTolerance;
    const G4ThreeVector sphereCenter(x0 + 0.5 * (x1 - x0), y0 + 0.5 * (y1 - y0), z0 + 0.5 * (z1 - z0));
    // Technique is to start outside of the geometry then trace through
    // saving each surface intercept point. One of these points will
    // then be selected at random to be the point used. This follows
    // this comment:
    // a complicated case: generate points uniformly on the
    // surface of the surrounding volume.  The mathematical technique
    // used here relies on what I think of as "Olber's Law":
    // given a surrounding surface with uniform # sources per unit area and
    // cos(theta) distribution of emitted rays (theta is angle from normal),
    // the brightness per unit area on the surface of a contained volume
    // is uniform (and the incident rays have a cos(theta) distribution,
    // although this is largely irrelevant for us).  This is true even for
    // non-convex solids if one considers all intercepts, such that the
    // interior solid does not shadow itself.

    // First find a trial point and ray direction:
    for( unsigned int iTrial = 0; iTrial < 10000; iTrial++ )
      {
        double u,v,w;
        do
        {
            u= G4UniformRand()*2.0-1.0;
            v= G4UniformRand()*2.0-1.0;
            w= 1.0- (u*u+v*v);
        } while(w < 0.0);
        w = sqrt(w);
        // (end of cute sequence.)
        // Now generate position on unit sphere:
        G4ThreeVector spos;
        double r2;
        do
        {
            spos = G4ThreeVector(G4UniformRand() * 2.0 -1.0,
                    G4UniformRand() * 2.0 - 1.0,
                    G4UniformRand() * 2.0 - 1.0);
            r2= spos.mag2();
        } while (r2 > 1.0 || r2 < 0.0625);
        spos*= sqrt(1.0/r2);
        // rotate direction to be relative to normal
        const G4ThreeVector e1 = spos.orthogonal().unit();
        const G4ThreeVector e2 = spos.cross(e1);
        const G4ThreeVector raydir = u*e1 + v*e2 - w*spos;
        // scale and offset position to be on surrounding sphere
        spos *= sphereRadius;
        spos += sphereCenter;
        // Now have a point (spos) and direction (raydir), find the intercepts
        vector<G4ThreeVector> intercepts;
        for( unsigned int iTrack = 0; iTrack < 20; iTrack++ )
          {
            if( solid->Inside( spos ) == kSurface ) // Move off the surface
              spos += 0.05 * raydir;
            if( solid->Inside( spos ) == kOutside )
              {
                const double distIn = solid->DistanceToIn( spos, raydir );
                if( distIn == kInfinity || distIn > 2.0 * sphereRadius ) // No intercepts left
                  break;
                else if( distIn > 0.0 )
                  {
                    spos += distIn * raydir;
                    intercepts.push_back( spos + solid->SurfaceNormal(spos) * G4UniformRand() * thickness );
                  }
                else // Something is wrong
                  {
                    warn << "GLG4PosGen_PointPaintFill: strange DistanceToIn "
                         << distIn << " at "
                         << spos.x() << ", " << spos.y() << ", " << spos.z()
                         << " in dir "
                         << raydir.x() << ", " << raydir.y() << ", " << raydir.z()
                         << " loop " << iTrack << endl;
                  }
              }
            else // Inside or on surface
              {
                const double distOut = solid->DistanceToOut( spos, raydir );
                if( distOut == kInfinity || distOut > 2.0 * sphereRadius ) // No intercepts left
                  break;
                else if( distOut > 0.0 )
                  {
                    spos += distOut * raydir;
                    intercepts.push_back( spos + solid->SurfaceNormal(spos) * G4UniformRand() * thickness );
                  }
                else // Something is wrong
                  {
                    warn << "GLG4PosGen_PointPaintFill: strange DistanceToOut "
                         << distOut << " " << " at "
                         << spos.x() << ", " << spos.y() << ", " << spos.z()
                         << " in dir "
                         << raydir.x() << ", " << raydir.y() << ", " << raydir.z()
                         << " loop " << iTrack << endl;
                  }
              }
          }
        // Now should have intercepts
        if( intercepts.empty() == false )
          {
            unsigned int iIntercept = static_cast<unsigned int>( G4UniformRand() * intercepts.size() );
            // Convert to global coordinates
            fVolumes.GetTransform(iVolume).ApplyPointTransform( intercepts[iIntercept] );
            // Set result and return
            argResult = intercepts[iIntercept];
            return;
          }
      }
    Log::Die( "PosGen_Fill::GeneratePosition failed in 100000 attempts to find position in volume " + fVolumes.GetVolume(iVolume)->GetName() );
}

double PosGen_Paint::ComputeArea(const G4VPhysicalVolume* volume)
{
    if(fSolidAreas.count(volume->GetName()))
        return fSolidAreas[volume->GetName()];
    else{
        double surfaceArea = volume->GetLogicalVolume()
            ->GetSolid()->GetSurfaceArea();
        fSolidAreas[volume->GetName()] = surfaceArea;
        return surfaceArea;
    }
}

std::map<std::string, double> PosGen_Paint::fSolidAreas;
} // namespace RAT