/**VertexGen_UFO
 * contact person:  Mark Stringer <ms711@sussex.ac.uk>
 */

#include <G4Event.hh>
#include <G4ParticleTable.hh>
#include <G4PrimaryParticle.hh>
#include <G4PrimaryVertex.hh>
#include <RAT/GLG4StringUtil.hh>

#include <Randomize.hh>
#include <RAT/DB.hh>
#include <RAT/EventInfo.hh>
#include <RAT/Log.hh>
#include <RAT/VertexGen_UFO.hh>

#include <TFile.h>
#include <TH1D.h>

#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>

#include <cmath>
using namespace CLHEP;
using namespace std;

namespace RAT {

    VertexGen_UFO::VertexGen_UFO(const char *arg_dbname) : GLG4VertexGen(arg_dbname)
    {
        // Not yet initialised - call the Init function to setup on first event
        fInitDB = false;
    }

    VertexGen_UFO::~VertexGen_UFO()
    {
    }


    void VertexGen_UFO::Init()
    {
        // This function should be run once on the first event (once DB all set up)
        warn << "VertexGen_UFO::Init: Initialising simulation settings \n";

        // Photon definition (assume /run/initialize already done)
        fPhotonDef = G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");

        // Load database values getting information about geometry of ufo first
        DBLinkPtr lUFO = DB::Get()->GetLink("UFOGEO");
        fNleds = lUFO->GetI("n_leds");
        fUFORadius = lUFO->GetD("radius");

        //Use these values to generate the LED Positionsj
        double angularSpacing = twopi/fNleds;
        //FIRST LED FIRING WILL BE POINTING EAST
        for(int ledIter=0; ledIter<fNleds; ledIter++){
            fUFO_XPos.push_back(fUFORadius*cos(ledIter*angularSpacing));
            fUFO_YPos.push_back(fUFORadius*sin(ledIter*angularSpacing));
            fUFO_ZPos.push_back(0.0);
        }

        //Loading the UFO wavelength distribution
        try{
            lUFO = DB::Get()->GetLink("UFOLED");
            info << "Loaded information about the UFO LED"<<newline;
        }
        catch(DBNotFoundError &e) {
            Log::Die("VertexGen_UFO::Init: can't find ratdb UFOLED");
        }

        //COPIED FROM VertexGen_Laserball
        // firstly get the wavelengths and store endpoints
        fwave = lUFO->GetFArrayFromD("dist_wl");
        fMaxWl = fwave.back();
        fMinWl = fwave.front();
        // make a check that the values in the dist_wl make sense and increase monotonically
        for(unsigned iv=1;iv<fwave.size(); iv++){
            if(fwave[iv]<=fwave[iv-1]){
                warn << "VertexGen_UFO::Init: dist_wl[" << iv << "] = " << fwave[iv] << " [" << iv-1 << "] = "  << fwave[iv-1] << "!\n";
                Log::Die("VertexGen_UFO::Init: dist_wl must increase monatonically. Please check values.");
            }
        }
        // now get the amplitudes
        try{
            fwaveInt = lUFO->GetFArrayFromD("dist_wl_intensity");
        }
        catch(DBNotFoundError &e){
            Log::Die("VertexGen_UFO::Init: can't find ratdb UFOLED dist_wl_intensity info");
        }
        // now create a cumulative magnitude vector
        float wavemagsum = 0;
        fwaveCumMag.clear();
        fwaveCumMag.push_back(0.0);
        info << " Cumulative wavelength ";
        for(unsigned int i=1;i<fwave.size();++i){
            wavemagsum +=(fwaveInt[i]+fwaveInt[i-1])/2.0 * (fwave[i]-fwave[i-1]);
            fwaveCumMag.push_back(wavemagsum);
            info << fwaveCumMag[i] << " ";
        }
        info << "\n";
        //Getting intensity from DB
        try{
            fPhotonsPerEvent = lUFO->GetI("intensity");
        }
        catch(DBNotFoundError &e){
            Log::Die("VertexGen_UFO::Init: cant find ratdb UFOLED intensity info");
        }
        //Copied from  VertexGen_Laserball
        try {
            string temp = lUFO->GetS("pulse_mode");
            if(temp == "POISSON"||temp == "poisson"){
                fPoisson = true;
                info <<  "VertexGen_UFO::Init: Setting pulse mode to poisson\n";
            }else if(temp == "fixed" || temp == "FIXED"){
                fPoisson = false;
                info <<  "VertexGen_UFO::Init: Setting pulse mode to fixed\n";
            }else{
                Log::Die(dformat("VertexGen_UFO::Init: UFO.pulse_mode = \"%s\" is invalid.",
                            temp.c_str()));
            }
        } catch (DBNotFoundError &e) { };
        //Now loading angle vector
        try{
            fangle = lUFO->GetFArrayFromD("dist_angle");
        }
        catch(DBNotFoundError &e){
            Log::Die("VertexGen_UFO::Init: can't find ratdb  UFOLED dist_angle info");
        }
        //Now generating angular intensity array as lambertian
        fAngInt = GenerateLambertian(fangle);

        // now create a cumulative magnitude vector (Taken from VertexGen_Laserball)
        float anglemagsum = 0;
        fCumAngleDist.clear();
        fCumAngleDist.push_back(0.0);
        info << " Cumulative Angle ";
        for(unsigned int i=1;i<fangle.size();++i){
            anglemagsum +=(fAngInt[i]+fAngInt[i-1])/2.0 * (fangle[i]-fangle[i-1]);
            fCumAngleDist.push_back(anglemagsum);
            info << fCumAngleDist[i] << " ";
        }
        info << "\n";

        //now obtaining time distribution
        try{
            fTimeWidth = lUFO->GetD("time_width");
        }
        catch(DBNotFoundError &e){
            Log::Die("VertexGen_UFO::Init: can't find ratdb UFOTIME time_width info");
        }
        //Now loading the timing information of the UFO
        try{
            lUFO = DB::Get()->GetLink("UFOTIME");
            info << "Loaded information about the UFO TIMING"<<newline;
        }
        catch(DBNotFoundError &e) {
            Log::Die("VertexGen_UFO::Init: can't find ratdb UFOTIME");
        }
        try{
            fRotationTime = lUFO->GetD("rotation_time");
        }
        catch(DBNotFoundError &e){
            Log::Die("VertexGen_UFO::Init: can't find ratdb  UFOTIME rotation_time info");
        }
        //Initializing flash counter
        fLEDFlashCount = 0;

        try{
            fFlashRate = lUFO->GetD("pulse_rate");
        }
        catch(DBNotFoundError &e){
            Log::Die("VertexGen_UFO::Init: can't find ratdb  UFOTIME pulse_rate  info");
        }
        if(fFlashRate < 0){
            Log::Die("VertexGen_UFO::Init please set UFOTIME pulse_rate to the generator rate");
        }
        fFlashesPerRotation = fFlashRate*fRotationTime;
        // Should now be all initialised so set flag
        fInitDB = true;
    }

    void VertexGen_UFO::BeginOfRun() {
      Init();
    }

    std::vector<float>   VertexGen_UFO::GenerateLambertian(std::vector<float> angularDistribution){
        std::vector<float> output;
        double area = sin(angularDistribution.back())-sin(angularDistribution.front());
        if(area == 0.0 || area == 0) area = 1.0;
        info << "AREA is:  "<<area<<newline;
        for(int i=0; i<angularDistribution.size(); i++){
            //Division by area is for normalization
            output.push_back(cos(angularDistribution[i])/area);
        }
        return output;
    }

    //TAKEN FROM VertexGen_Laserball
    float VertexGen_UFO::SampleRandom(std::vector<float> mydistx, std::vector<float> mydistraw, std::vector<float> mydistcum){
        // Return value sampled from distribution with cumulative distribution mydistcum, raw distribution mydistraw for x points mydistx
        // (assume linear interpolation between specified points in spectrum)
        float start = mydistcum.front();
        float stop = mydistcum.back();
        // now throw a random number between these limits
        float random = start+G4UniformRand()*(stop-start);
        // now loop through cumulative distribution to choose energy
        int istep = 0;
        while(mydistcum[istep+1]<random){
            istep++;
        }

        // The desired value is bracketed by istep and istep+1 points.
        // Since the cumulative function is an integral of a linearly
        // interpolated spectrum, it is quadratic between data points.
        // Find x assuming quadratic shape.

        float spec_slope = (mydistraw[istep+1] - mydistraw[istep]) / (mydistx[istep+1] - mydistx[istep]);

        float x;
        float delta_y = random - mydistcum[istep];
        if (fabs(spec_slope) < 1e-6) {
            // ?? not sure I get this step?
            x = mydistx[istep] + delta_y/mydistraw[istep];

        } else {
            float a = 0.5 * spec_slope;
            float b = mydistraw[istep] - spec_slope * mydistx[istep];
            float c = -a * mydistx[istep] * mydistx[istep] - b * mydistx[istep] - delta_y;
            x = (-b + sqrt(b*b - 4 * a * c))/(2*a);
        }

        return x;
    }

    G4ThreeVector VertexGen_UFO::GeneratePhotonDirection(G4ThreeVector led_dir, float angle){
        //Getting the unit vector of the LED
        led_dir.setMag(1.0);
        //z component of vector is zero (in coordinate system local to UFO) hence do a rotation using standard 2D rotation matrix
        double photon_dirx = led_dir.x()*cos(angle)-led_dir.y()*sin(angle);
        double photon_diry = led_dir.x()*sin(angle)+led_dir.y()*cos(angle);
        G4ThreeVector photon_dir = G4ThreeVector(photon_dirx,photon_diry,0.0);
        //Generate random phi roation
        double rand_phi = G4UniformRand()*CLHEP::twopi;
        photon_dir.rotate(rand_phi,led_dir);
        return photon_dir;
    }

    /**Taken from vertexGen_Laserball*/
    double VertexGen_UFO::GasDev(){

        // This routine is from snoman where D. Wark from Oxford obtained it by
        // modifying the Numerical Recipes version.

        double v1, v2, r, fac;
        r = 1.02;
        while(r>= 1.0){
            v1 = 2.0 * G4UniformRand() - 1.0;
            v2 = 2.0 * G4UniformRand() - 1.0;
            r = v1 * v1 + v2 * v2;
        }
        fac = v2 * sqrt( -2.0 * log( r ) / r );
        return fac;
    }
    void VertexGen_UFO::GeneratePrimaryVertex(G4Event *event,
            G4ThreeVector &fPos0,
            G4double fTime0)
    {
        if (!fInitDB) Init();
        // How many photons?
        int nphot = 0;
        if(fPoisson){
            nphot = (G4int)( CLHEP::RandPoisson::shoot(fPhotonsPerEvent) );
        }else{
            nphot = fPhotonsPerEvent;
        }
        //Calculate which LED is firing

        int ledNumber = (int)((((double)(fLEDFlashCount%fFlashesPerRotation))/fFlashesPerRotation)*fNleds);
        fLEDFlashCount++;
        G4ThreeVector ledPos;
        for(int i=0; i<nphot; i++){
            //Calculating time spread on LED Flash
            double time = GasDev()*fTimeWidth;
            //Getting UFO led  position
            ledPos = G4ThreeVector(fUFO_XPos[ledNumber],fUFO_YPos[ledNumber],fUFO_ZPos[ledNumber]);
            float  wavelength = SampleRandom(fwave, fwaveInt, fwaveCumMag);
            float energy = hbarc * twopi / (wavelength * nm);
            float momentum = energy; // GEANT uses momentum in same units as energy
            //Getting Lambertian angle
            float  emissionAngle = SampleRandom(fangle,fAngInt,fCumAngleDist);
            G4ThreeVector dir = GeneratePhotonDirection(ledPos,emissionAngle);

            //Taken from vertex gen laserball
            G4ThreeVector mom(dir);
            mom.setMag(momentum);    // Scale to right magnitude

            //Adding the led direction to the centre of the UFO
            G4PrimaryVertex* vertex= new G4PrimaryVertex(fPos0+ledPos, fTime0 + time);
            G4PrimaryParticle* particle = new G4PrimaryParticle(fPhotonDef,
                    mom.x(),mom.y(),mom.z());

            // Generate random polarization
            double phi = (G4UniformRand()*2.0-1.0)*pi;
            G4ThreeVector e1 = mom.orthogonal().unit();
            G4ThreeVector e2 = mom.unit().cross(e1);
            G4ThreeVector pol = e1*cos(phi)+e2*sin(phi);
            particle->SetPolarization(pol.x(), pol.y(), pol.z());
            particle->SetMass(0.0); // Seems odd, but used in GLG4VertexGen_Gun


            vertex->SetPrimary(particle);
            event->AddPrimaryVertex(vertex);

        }


    }

    void VertexGen_UFO::SetState(G4String newValues)
    {

        Log::Die("SetState is deprecated, do not call it. Set all variables with /rat/db/set.");
        return;

    }


    G4String VertexGen_UFO::GetState()
    {
        stringstream result;
        result << "Photons per pulse: " <<fPhotonsPerEvent <<newline;
        result << "Rotation Time: "<<fRotationTime<<newline;
        return result.str();
    }

}