// VertexGen_ELLIE.cc
// Contact person:Jeanne Wilson <j.r.wilson@qmul.ac.uk>
// See VertexGen_ELLIE.hh for more details
//———————————————————————//

// This file contains the bulk of code that used to be in Gen_ELLIE.cc -
// split up so can now be used to generate ELLIE events
// as a separate vertex. It is meant to be used with the PosGen_ELLIE position
// generator with the same fibre_id argument for both vertex and position
// to simulate ELLIE events. A separate vertex and position generator are
// required as input to the coincidence generator.
////////////////////////////////////////////////////////////////////////
#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>

#include <RAT/VertexGen_ELLIE.hh>
#include <RAT/EventInfo.hh>
#include <RAT/Log.hh>

#include <RAT/GLG4PosGen.hh>
#include <RAT/GLG4TimeGen.hh>
#include <RAT/Factory.hh>
#include <RAT/GLG4StringUtil.hh>

#include <Randomize.hh>
#include <G4Event.hh>
#include <G4ParticleDefinition.hh>
#include <G4ParticleTable.hh>
#include <G4PrimaryVertex.hh>
#include <G4PrimaryParticle.hh>
#include <G4ThreeVector.hh>
#include <RAT/LinearInterp.hh>
#include <RAT/Factory.hh>
#include <RAT/DB.hh>
#include <RAT/PhotonThinning.hh>

using namespace CLHEP;
using namespace std;

namespace RAT {


VertexGen_ELLIE::VertexGen_ELLIE(const char *arg_dbname) : GLG4VertexGen(arg_dbname)
{
    // Not yet initialised - call the Init function to setup on first event
    fInitDB = false;
    // Initialise the fibre to null - set either by state or as DB command
    fFibreID="";
}

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

    // Photon definition (assume /run/initialize already done)
    fPhotonDef = G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
    // Is there a photon thinning factor - this doesn't necessarily make sense for optical simulations
    double photon_thin =  RAT::PhotonThinning::GetFactor();
    if(photon_thin > 1.0){
        // OUTPUT WARNING
        warn << "VertexGen_ELLIE::Init: !!!WARNING!!! - You are using photon thinning with factor " << photon_thin
	       << "\n Are you sure you want to do this for an Optical Simulation (could result in loss of accuracy if low number of photons) \n";
    }

    // Load database values
    DBLinkPtr lELLIE = DB::Get()->GetLink("ELLIE");
    fPhotonsPerEvent = int(double(lELLIE->GetI("intensity"))/photon_thin);
    debug << "VertexGen_ELLIE::Init:  intensity = " << fPhotonsPerEvent << "\n";

    // Optional setting
    try {
        string temp = lELLIE->GetS("pulse_mode");
        if(temp == "poisson"){
            fPoisson = true;
        }else if(temp == "fixed"){
            fPoisson = false;
        }else{
            Log::Die(dformat("VertexGen_ELLIE::Init: ELLIE.pulse_mode = \"%s\" is invalid.",
			temp.c_str()));
        }
    } catch (DBNotFoundError &e) { };

    // Which fibre to simulate - identify by name:
    if(fFibreID==""){
        // It hasn't been set as the state of the generator so load from DB
        fFibreID = lELLIE->GetS("fibre_id");
    }

    // Now get the info specific for this fibre
    DBLinkPtr lfibre;
    try{
        lfibre = DB::Get()->GetLink("FIBRE", fFibreID);
        info <<  "VertexGen_ELLIE::Init: Getting fibre info for " << fFibreID << "\n";
    }catch(DBNotFoundError &e) {
        Log::Die("VertexGen_ELLIE::Init: can't find ratdb FIBRE info for " + fFibreID);
    }

    // ELLIE directions
    try {
        fDirX = lfibre->GetD("u");
        fDirY = lfibre->GetD("v");
        fDirZ = lfibre->GetD("w");
    } catch (DBNotFoundError &e) {
        Log::Die("VertexGen_ELLIE::Init: No directions set for fibre " + fFibreID);
    };

    // wavelength info
    try{
        fMono_wls = lELLIE->GetI("mono_wl");
    } catch (DBNotFoundError &e) {  };
    if(!fMono_wls){
        // chose wavelength either set by user (SMELLIE) or use default for that ELLIE fibre (AMELLIE/TELLIE)
        try{
            fELLIEWavelength = lELLIE->GetS("wavelength_dist");
            info << "VertexGen_ELLIE::Init: Setting wavelength distribution to " + fELLIEWavelength + "\n";
        } catch (DBNotFoundError &e) {
            info << "VertexGen_ELLIE::Init: Using fixed wavelengths as set in ratdb FIBRE wavelength field\n";
            fELLIEWavelength = "FIBREDEFAULT";
        }
        if(fELLIEWavelength=="FIBREDEFAULT"){
            try{
                fELLIEWavelength = lfibre->GetS("wavelength_dist");
            }catch (DBNotFoundError &e) {
                // wavelength not specified in either place
                Log::Die("VertexGen_ELLIE::Init: can't find wavelength specification for fibre " + fFibreID);
            };
        };
        // Find the right database file with this index (need exception catching here)
        debug <<  "VertexGen_ELLIE::Init: Getting Wavelength info for " << fELLIEWavelength << "\n";
        DBLinkPtr lELLIEwave;
        try{
            lELLIEwave = DB::Get()->GetLink("ELLIEWAVE", fELLIEWavelength);
        }catch(DBNotFoundError &e) {
            Log::Die("VertexGen_ELLIE::Init: can't find ratdb wavelength distribution for " + fELLIEWavelength);
        };
        // firstly get the wavelengths and store endpoints
        fwave = lELLIEwave->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 int iv=1;iv<fwave.size(); iv++){
            if(fwave[iv]<=fwave[iv-1]){
                warn << "VertexGen_ELLIE::Init: dist_wl[" << iv << "] = " << fwave[iv] << " [" << iv-1 << "] = "  << fwave[iv-1] << "!\n";
                Log::Die("VertexGen_ELLIE::Init: dist_wl must increase monotonically. Please check values.");
            }
        }
        // now get the amplitudes
        fwaveInt = lELLIEwave->GetFArrayFromD("dist_wl_intensity");
        // now create a cumulative magnitude vector
        float wavemagsum = 0;
        fwaveCumMag.clear();
        fwaveCumMag.push_back(0.0);
        detail << " 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);
            detail << fwaveCumMag[i] << " ";
        }
        detail << "\n";
    }


    // timing info
    fMono_time = 0;
    try{
        fELLIETime = lELLIE->GetS("time_dist");
        info << "VertexGen_ELLIE::Init WARNING: Setting time distribution to " + fELLIETime + "\n";
    }catch (DBNotFoundError &e) {
        info << "VertexGen_ELLIE::Init: Using fixed time distribution as set in ratdb FIBRE time field\n";
        fELLIETime = "FIBREDEFAULT";
    }
    if(fELLIETime=="mono"){
        fMono_time = 1;
        info << "VertexGen_ELLIE::Init WARNING: Setting time distribution to mono \n";
    }else if(fELLIETime=="FIBREDEFAULT"){
        try{
            fELLIETime = lfibre->GetS("time_dist");
        }catch (DBNotFoundError &e){
            // time distribution not specified - check for parameter definition
            try {
                fELLIEIPW = lfibre->GetD("ipw");
            } catch (DBNotFoundError &e) {
                Log::Die("VertexGen_ELLIE::Init: can't find time IPW for fibre " + fFibreID);
            }
            try {
                fELLIEPulserBoard = lfibre->GetS("pulser_board");
            } catch (DBNotFoundError &e) {
                Log::Die("VertexGen_ELLIE::Init: can't find time pulser-board for fibre " + fFibreID);
            }
            fELLIETime="FIBREPARAM";
            info << "VertexGen_ELLIE::Init info: Using parameterised timing distribution for " << fFibreID << " with IPW = " << fELLIEIPW << ", pulserboard = " <<fELLIEPulserBoard << "\n";
                debug << "Though actually its just a dummy distribution at the moment FIX ME!\n";
        }
    }
    if(fELLIETime=="FIBREPARAM"){
        // Extract the information to define the timing parameterisation
        DBLinkPtr lpulserboard;
        try {
            lpulserboard = DB::Get()->GetLink("ELLIEPULSERBOARD",fELLIEPulserBoard);
            debug << "VertexGen_ELLIE::Init: Getting Pulserboard info for board " << fELLIEPulserBoard << "\n";
        } catch (DBNotFoundError &e) {
            Log::Die("VertexGen_ELLIE::Init: can't find pulser board info for board " + fELLIEPulserBoard);
        }
        ////FIXME!!!!////
        // Pull the parameters out of the PulserBoard table here
        // Then do something with them (ideally make a set of vectors using SampleRandom)
        // - for now put in a dummy distribution.
        ftime.push_back(0.2);
        ftime.push_back(0.8);
        ftime.push_back(1.0);
        ftimeInt.push_back(1.0);
        ftimeInt.push_back(0.5);
        ftimeInt.push_back(0.0);
        float timemagsum = 0;
        ftimeCumMag.clear();
        ftimeCumMag.push_back(0.0);
        detail << " Cumulative time dummy dist";
        for(unsigned int i=1;i<ftime.size();++i){
            timemagsum +=(ftimeInt[i]+ftimeInt[i-1])/2.0 * (ftime[i]-ftime[i-1]);
            ftimeCumMag.push_back(timemagsum);
            detail << ftimeCumMag[i] << " ";
        }
        detail << "\n";
    }else if(!fMono_time){
        // Find the right database file with this index
        DBLinkPtr lELLIEtime;
        try{
            lELLIEtime = DB::Get()->GetLink("ELLIETIME", fELLIETime);
            debug <<  "VertexGen_ELLIE::Init: Getting Time info for " << fELLIETime << "\n";
        }catch(DBNotFoundError &e) {
            Log::Die("VertexGen_ELLIE::Init: can't find ratdb time distribution for " + fELLIETime);
        }
        // firstly get the times and store endpoints
        ftime = lELLIEtime->GetFArrayFromD("dist_time");
        fMaxTime = ftime.back();
        fMinTime = ftime.front();
        // make a check that the values in the dist_time make sense and increase monotonically
        for(unsigned int iv=1;iv<ftime.size(); iv++){
            if(ftime[iv]<=ftime[iv-1]){
                warn << "VertexGen_ELLIE::Init: dist_time[" << iv << "] = " << ftime[iv] << " [" << iv-1 << "] = "  << ftime[iv-1] << "!\n";
                Log::Die("VertexGen_ELLIE::Init: dist_time must increase monotonically. Please check values.");
            }
        }
        // now get the amplitudes
        ftimeInt = lELLIEtime->GetFArrayFromD("dist_time_intensity");
        // now create a cumulative magnitude vector
        float timemagsum = 0;
        ftimeCumMag.clear();
        ftimeCumMag.push_back(0.0);
        detail << " Cumulative time ";
        for(unsigned int i=1;i<ftime.size();++i){
            timemagsum +=(ftimeInt[i]+ftimeInt[i-1])/2.0 * (ftime[i]-ftime[i-1]);
            ftimeCumMag.push_back(timemagsum);
            detail << ftimeCumMag[i] << " ";
        }
        detail << "\n";
    }

    // angular info
    try{
        fELLIEAngle = lELLIE->GetS("angle_dist");
        info << "VertexGen_ELLIE::Init: Setting angular distribution to " + fELLIEAngle + "\n";
    }catch (DBNotFoundError &e) {
        info << "VertexGen_ELLIE::Init: Using fixed angular distribution as set in ratdb FIBRE time field\n";
        fELLIEAngle="FIBREDEFAULT";
    }
    fIso_angle = 0;
    fZero_angle = 0;
    if(fELLIEAngle=="iso"){
        fIso_angle = 1;
        info << "VertexGen_ELLIE::Init: Setting angular distribution to isotropic \n";
    }else if(fELLIEAngle=="zero"){
        fZero_angle = 1;
        info << "VertexGen_ELLIE::Init: Setting angular distribution to zero (ie no dispersion) \n";
        if(fIso_angle&&fZero_angle)Log::Die("VertexGen_ELLIE: Beam is set to be both isotropic and Zero dispersion!");
    }else{

        if(fELLIEAngle=="FIBREDEFAULT"){
            try{
                fELLIEAngle = lfibre->GetS("angle_dist");
            }catch (DBNotFoundError &e){
                // angular distribution not specified in either place
                Log::Die("VertexGen_ELLIE::Init: can't find angular information for fibre " + fFibreID);
            }
        }
         // Find the right database file with this index (need exception catching here)
        DBLinkPtr lELLIEangle;
        try{
            lELLIEangle = DB::Get()->GetLink("ELLIEANGLE", fELLIEAngle);
            debug <<  "VertexGen_ELLIE::Init: Getting Angular info for " << fELLIEAngle << "\n";
        }catch(DBNotFoundError &e) {
            Log::Die("VertexGen_ELLIE::Init: can't find ratdb angular distribution for " + fELLIEAngle);
        }
        string angle_units = lELLIEangle->GetS("dist_angle_option");
        if(angle_units == "degrees" ){
            fConvertAngle = twopi/360.;
        }else if(angle_units == "radians"){
            fConvertAngle = 1;
        }else{
            fConvertAngle = 1;
            warn << "VertexGen_ELLIE::Init: Unknown angular unit, assuming radians\n";
        }
        // firstly get the angles and store endpoints
        fang = lELLIEangle->GetFArrayFromD("dist_angle");
        fMaxAngle = fang.back();
        fMinAngle = fang.front();
        // make a check that the values in the dist_wl make sense and increase monotonically
        for(unsigned int iv=1;iv<fang.size(); iv++){
            if(fang[iv]<=fang[iv-1]){
                warn << "VertexGen_ELLIE::Init: dist_angle[" << iv << "] = " << fang[iv] << " [" << iv-1 << "] = "  << fang[iv-1] << "!\n";
                Log::Die("VertexGen_ELLIE::Init: dist_angle must increase monotonically. Please check values.");
            }
        }
        if(fMaxAngle-fMinAngle<-0 || fMaxAngle<0 || fMinAngle<0){
            // problem with angular range
            Log::Die("VertexGen_ELLIE: Nonsensical angular range for " + fELLIEAngle);
        }
        // now get the amplitudes
        fangInt = lELLIEangle->GetFArrayFromD("dist_angle_intensity");
        // now create a cumulative magnitude vector and account for solid angle weighting
        float angmagsum = 0;
        fangCumMag.clear();
        fangCumMag.push_back(0.0);
        fangInt[0] =fangInt[0]*sin(fang[0]*fConvertAngle);
        detail << " Cumulative angle ";
        for(unsigned int i=1;i<fang.size();++i){
            fangInt[i]=fangInt[i]*sin(fang[i]*fConvertAngle);
            angmagsum +=(fangInt[i]+fangInt[i-1])/2.0 * (fang[i]-fang[i-1]);
            fangCumMag.push_back(angmagsum);
            detail << fangCumMag[i] << " ";
        }
        detail << "\n";
     }

    // Should now be all initialised so set flag
    fInitDB = true;

}



VertexGen_ELLIE::~VertexGen_ELLIE()
{
}

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

void VertexGen_ELLIE::GeneratePrimaryVertex(G4Event *event,
                                              G4ThreeVector& dx,
                                              G4double dt)
 {
  if(!fInitDB)Init();

     G4ThreeVector dir(fDirX,fDirY,fDirZ);
  G4ThreeVector normal = dir.unit();
  G4ThreeVector perp = normal.orthogonal().unit();

  // How many photons?
  int nphot = 0;
  if(fPoisson){
    nphot = (G4int)( CLHEP::RandPoisson::shoot(fPhotonsPerEvent) );
      if(nphot<=0){
          detail << "VertexGen_ELLIE::GeneratePrimaryVertex:  Poisson throw of intensity results in zero photons, simulating 1 \n";
          nphot = 1;
      }
  }else{
    nphot = fPhotonsPerEvent;
  }

  // sample the wavelength up front for writing out to calib branch
  float wavelength;
  if(fMono_wls){
    wavelength = fMono_wls; // nm
    debug << "VertexGen_ELLIE::GeneratePrimaryVertex:  ELLIE wavelength fixed = " << wavelength <<"\n";
  }else{
    debug << "VertexGen_ELLIE::GeneratePrimaryVertex:  ELLIE wavelength randomly chosen = ";
    wavelength = SampleRandom(fwave, fwaveInt, fwaveCumMag);
    //wavelength = fRandWl->shoot() * (fMaxWl - fMinWl) + fMinWl; // randomly sample for recording in MC file
    debug  << wavelength <<"\n";
  }
  // and sample the time
  float time = 0.;
  if(!fMono_time){
    time = SampleRandom(ftime, ftimeInt, ftimeCumMag);
    //time = (fRandTime->shoot() * (fMaxTime - fMinTime) + fMinTime);
  }
  // Record in extra G4Event information
  EventInfo *exinfo = dynamic_cast<EventInfo*>(event->GetUserInformation());
  DS::Calib *calib = new DS::Calib();
  calib->SetSourceName(fFibreID);       // set name as fibre ID now
     calib->SetID(1);                   // This number is non zero if there is a source, but what ID to assign???
  calib->SetMode((int) wavelength);     // sampled wavelength distribution
  calib->SetIntensity(nphot);
  calib->SetDir( TVector3(normal.x(), normal.y(), normal.z()) );
  calib->SetPos( TVector3(dx.x(), dx.y(), dx.z()) );
  calib->SetTime(time);
  exinfo->SetCalib(calib);

  // Add vertices for each photon
  for (int i=0; i < nphot; i++) {
	// photon time
	double t1 = 0;  // stays as 0 for Mono_time
	if(!fMono_time){
        t1 = SampleRandom(ftime, ftimeInt, ftimeCumMag);
	}
	if(!fMono_wls){ // set above for Mono_wls
        wavelength = SampleRandom(fwave, fwaveInt, fwaveCumMag);
	}

    float energy = hbarc * twopi / (wavelength * nm);
    float momentum = energy; // GEANT uses momentum in same units as energy

	// photon direction
    double theta;
    if (fIso_angle){
        theta = acos(0.9999 * (2.0*G4UniformRand()-1.0));
    }else if(fZero_angle){
        theta = 0;
    }else{
        theta = fConvertAngle*SampleRandom(fang, fangInt, fangCumMag);
    }

    double phi = twopi * G4UniformRand();

    // Calc momentum
    G4ThreeVector mom(normal);
    mom.rotate(theta, perp); // Rotate away from ELLIE direction (polar angle)
    mom.rotate(phi, normal); // Rotate around ELLIE direction (phi)

    mom.setMag(momentum);    // Scale to right magnitude

    G4PrimaryVertex *vertex = new G4PrimaryVertex(dx, dt+t1);
    G4PrimaryParticle *particle = new G4PrimaryParticle(fPhotonDef,
							mom.x(),
							mom.y(),
							mom.z());
    // Generate random polarization
    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);
  }

}

float VertexGen_ELLIE::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;

}


void VertexGen_ELLIE::SetState( G4String newValues )
{
    newValues = trim(newValues);
    if(newValues.length() == 0)// Empty arg, print help
    {
        detail << "Format of argument to VertexGen_ELLIE::SetState:" << endl
        << " \n\"fFibreID\" which is index of a ratdb FIBRE table"
        << " default will use value from ELLIE ratdb table "<< endl;
        Log::Die("Please set correctly parameter for the VertexGen_ELLIE "+GetState());
        return;
    }
    else
    {
        fFibreID = newValues;
        debug << "Current state of this VertexGen_ELLIE:\n\"" << GetState() << " \"\n";
    }
}


G4String VertexGen_ELLIE::GetState()
{
    stringstream result;
    result << fFibreID;
    return result.str();
}

} // namespace RAT