// PosGen_FillExp.cc
// Contact person: Aleksandra Bialek <abialek@ualberta.ca>
// See PosGen_FillExp.hh for more details
//———————————————————————//

#include <RAT/string_utilities.hpp>
#include <RAT/StringUtil.hh>
#include <RAT/PosGen_FillExp.hh>
#include <RAT/Log.hh>

#include <G4VPhysicalVolume.hh>
#include <G4LogicalVolume.hh>
#include <G4VSolid.hh>
#include <G4VoxelLimits.hh>
#include <G4AffineTransform.hh>
#include <G4GeometryTolerance.hh>
#include <Randomize.hh>
#include <G4Navigator.hh>
#include <G4TransportationManager.hh>
#include <G4PhysicalConstants.hh>

#include <TF1.h>
#include <TAxis.h>
#include <TMath.h>

#include <sstream>
#include <iostream>
#include <stdlib.h>

using std::string;
using std::vector;
using std::cout;
using std::endl;

namespace RAT {
void PosGen_FillExp::SetState(G4String newValues)
{
    vector<string> arg = split(newValues, " ");

    if(arg.size() == 0) // Empty arg, print help
    {
        fStates.push_back(" ");
        fStateDR.push_back(" ");
        info << "Current state of this PosGen_FillExp:\n\""
             << GetState() << " \"\n"
             << "Format of argument to PosGen_FillExp:\n"
             << " x y z dr\" a point in a volume in [mm],\n"
             << " dr - decay rate of the exponential radial distribution ('-' decay, '+' increase)\n"
             << " or\n\"name dr\", the name of a volume and the decay rate \n";
        Log::Die("Please set correctly parameters for the PosGen_FillExp "+GetState());
        return;
    }
    else if(arg.size() == 3 || arg.size()== 1) // 'dr' not given
    {
        info << "Format of arguments for PosGen_FillExp is missing 'dr':\n"
             << " dr - decay rate of the exponential radial distribution ('-' decay, '+' increase)\n"
             << "\"x y z dr\"  position within the volume and decay rate\n"
             << " or\n \"volname dr\", the name of a volume and the decay rate "
             << "\ntaking default value of 'dr'\n\n" ;

        string defDR = "0.0005";///default value
        fStateDR.push_back(defDR);
        fStates.push_back(join(arg," "));

        info << "Current state of this PosGen_FillExp:\n\""
             << GetState() << " \"\n";
    }
    else
    {
        fStateDR.push_back(arg.back());
        arg.pop_back();
        fStates.push_back(join(arg," "));
        info << "Current state of this PosGen_FillExp:\n\""
             << GetState() << " \"\n";
    }
    fListValid = false;
}

G4String PosGen_FillExp::GetState() const
{
    std::stringstream result;
    for(unsigned int it=0; it<fStates.size(); it++)
        result << fStates.at(it) << " " << fStateDR.at(it) ;
    return result.str();
}

void PosGen_FillExp::GeneratePosition(G4ThreeVector& argResult)
{
    // Check if must regenerate list
    if(!fListValid)
    {
        fVolumes.Clear();
        fTotalVolume = 0.0;
        std::string sdecayrate;
        for(unsigned int i = 0; i < fStates.size(); i++)
        {
            sdecayrate = fStateDR.at(i);
            fDecayRates.push_back(to_double(sdecayrate));
            fVolumes.AddString(fStates[i], VolumeList::SELF);
        }

        for(unsigned int i = 0; i < fVolumes.Size(); i++)
            fTotalVolume += fVolumes.ComputeVolume(fVolumes.GetVolume(i));

        fListValid = true;
    }

    // Now can generate a point, first choose the volume
    const double trialVolume = G4UniformRand() * fTotalVolume;
    int iVolume = -1; //iVolume incremented each time through, including first
    for(double cumulativeVolume = 0; trialVolume > cumulativeVolume;)
    {
        ++iVolume;
        cumulativeVolume += fVolumes.ComputeVolume(fVolumes.GetVolume(iVolume));
    }

    static G4Navigator *fNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();

    // calculate the size of the volume
    double volsize= fVolumes.ComputeVolume(fVolumes.GetVolume(iVolume));
    // Extract the solid
    G4VSolid* solid = fVolumes.GetVolume(iVolume)->GetLogicalVolume()->GetSolid();

    // Now calculate the bounding box, in local coordinates (to solid)
    double xmin, ymin, zmin, xmax, ymax, zmax;

    G4VoxelLimits voxelLimits;  // Defaults to "infinite" limits.
    G4AffineTransform affineTransform; // Defaults to no transform
    solid->CalculateExtent(kXAxis, voxelLimits, affineTransform, xmin,xmax);
    solid->CalculateExtent(kYAxis, voxelLimits, affineTransform, ymin,ymax);
    solid->CalculateExtent(kZAxis, voxelLimits, affineTransform, zmin,zmax);
    //compute the size of extent
    double extSize=(xmax-xmin)*(ymax-ymin)*(zmax-zmin);

    //calculate the ratio between sizes of logical volume and and the extent
    double sizeRatio = extSize/volsize;

    // maximum number of tries in the loop
    unsigned int maxtries = 100000 * (unsigned int)ceil(sizeRatio);

    double rmax=22600.;//~ max radius of the sphere over the cavity (circumradius)
    double r0=0., r1=rmax; //to range the distribution of the function

    if(TMath::Abs(zmin)<13000 && TMath::Abs(zmax)<13000) r1=13500; // for the smaller extent than cavity's

    const double decayrate = fDecayRates[iVolume];  // taken or from macro or default value

    TF1 funcRdecay("Rdecay","exp(x*[0])",0,rmax);
    funcRdecay.SetParameter(0,decayrate);
    funcRdecay.GetXaxis()->SetRangeUser(r0,r1);

    double r, phi, theta, rand, xx, yy, zz;
    unsigned int iTrial;
    // Now test positions until one is found in the relevant volume
    for(iTrial = 0; iTrial < maxtries; iTrial++)
    {
        // Now the trialposition is generated
        r = funcRdecay.GetRandom();
        phi = CLHEP::twopi * G4UniformRand();
        rand= G4UniformRand();
        theta = acos(2.0 * rand - 1);
        if(rand>0.5) theta = -theta;

        xx= r * sin(theta) * cos(phi);
        yy= r * sin(theta) * sin(phi);
        zz= r * cos(theta);

        G4ThreeVector trialPos( xx, yy, zz);

        // Convert this to global coordinates
        fVolumes.GetTransform(iVolume).ApplyPointTransform(trialPos);

        // NFB : The old logic was deeply flawed!!!
        //       The issue is that the solids are not aware of placement
        //       Calling  G4VSolid::Inside tests if a translation is inside the solid
        //       this means that it is the same as putting the solid centered in the origin
        //       This means that any solid that is placed off origin will be valid here
        //       the most obvious example are the NCD anchors.
        // This was the origin of issue #1392 on github
        //

        // The best way to check that the point belongs to the volume is actually using the navigator
        //
        G4VPhysicalVolume* trialVolumePhys = fNavigator->LocateGlobalPointAndSetup(trialPos,NULL,false,true);

        if (trialVolumePhys == fVolumes.GetVolume(iVolume)) {
          // This is the solid we want.
          // Set result and return
          argResult = trialPos;
          return;

        }
    }
    std::string volName = fVolumes.GetVolume(iVolume)->GetName();
    Log::Die(dformat("PosGen_Fill::GeneratePosition failed in %d attempts to find position in volume %s ",iTrial, volName.c_str() ) );
}
} // namespace RAT