// VertexGen_Laserball.cc
// Contact person: Jose Maneira <maneira@lip.pt>
// See VertexGen_Laserball.hh for more details

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

#include <TMath.h>

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

#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>
using CLHEP::halfpi;

using namespace CLHEP;
using namespace std;

namespace RAT {

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

//*****************************************************************************************
//*****************************************************************************************

VertexGen_Laserball::~VertexGen_Laserball()
{
}

//*****************************************************************************************
//*****************************************************************************************


void VertexGen_Laserball::Init()
{
    // This function should be run once on the first event (once DB all set up)
    warn << "VertexGen_Laserball::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_Laserball::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 lLASERBALL = DB::Get()->GetLink("LASERBALL_SIMULATION");

    try { fSimulationMode = lLASERBALL->GetD("simulation_mode"); }
    catch( DBNotFoundError& e ) { Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL simulation_mode info"); };

    if (fSimulationMode == 3)  Log::Die("VertexGen_LaserBall::Init: Loading laserball parameters from Camera tables not implemented yet.");

    // Load blocks of parameters for the LB
    //  1) Wavelength
    //  2) Position and Orientation (see also PosGen_Laserball)
    //  3) Angular distribution
    //  4) Intensity
    //  5) Timing

    // Load parameters from Manip if simulation_mode = 2
    if (fSimulationMode == 2){
        DBLinkPtr lCAL_SOURCE = DB::Get()->GetLink("CALIB_COMMON_RUN_LEVEL","MANIP");
        DBLinkPtr lN2DYELASER = DB::Get()->GetLink("N2DYELASER_RUN_LEVEL");

        // 1) Wavelength
        try{
            fDyeCell = lN2DYELASER->GetI("dyelaser_cell_nbr");
            info << "VertexGen_Laserball::Init: Setting dye cell nbr to " << fDyeCell<< "\n";
        } catch (DBNotFoundError &e) {
            Log::Die("VertexGen_Laserball::Init: can't find ratdb N2DYELASER_RUN_LEVEL dye_cell_nbr info");
        }
        if (fDyeCell < 0 || fDyeCell > 5)
            Log::Die(dformat("VertexGen_Laserball::Init: N2DYELASER_RUN_LEVEL dye_cell_nbr \"%d\" out of range.",fDyeCell));

        try{
            fDyeCellSet = lN2DYELASER->GetS("dyelaser_cell_set");
            info << "VertexGen_Laserball::Init: Setting dye cell set to " + fDyeCellSet + "\n";
        } catch (DBNotFoundError &e) {
            Log::Die("VertexGen_Laserball::Init: can't find ratdb N2DYELASER_RUN_LEVEL dye_cell_set info");
        }

        DBLinkPtr lDYECELL;
        try{
            lDYECELL = DB::Get()->GetLink("DYECELL",fDyeCellSet);
            info <<  "VertexGen_Laserball::Init: Getting dye cell set info for " << fDyeCellSet << "\n";
        }catch(DBNotFoundError &e) {
            Log::Die("VertexGen_Laserball::Init: can't find ratdb DYECELL info for " + fDyeCellSet);
        }

        try{
            fDyeCellOrder = lDYECELL->GetSArray("dye_order");
            info << "VertexGen_Laserball::Init: Setting dye cell order to "<< fDyeCellOrder << "\n";
        } catch (DBNotFoundError &e) {
            Log::Die("VertexGen_Laserball::Init: can't find ratdb DYECELL dye_order info");
        }
        fDyeName = fDyeCellOrder[fDyeCell];

        // 2) Position and Orientation
        // x,y,z is loaded in PosGen_Laserball
        try{
            fOrientation = lCAL_SOURCE->GetD("orientation");
            info <<  "VertexGen_Laserball::Init: Orientation: "<< fOrientation << "\n";
        } catch (DBNotFoundError &e) {
            Log::Die("VertexGen_Laserball::Init: can't find ratdb CALIB_COMMON_RUN_LEVEL orientation info");
        };
        if (fOrientation < 0 || fOrientation > 4) {
            Log::Die(dformat("VertexGen_Laserball::Init: CALIB_COMMON_RUN_LEVEL orientation  \"%d\" not recognized.",fOrientation));
        }

    }
    // If simulation_mode = 0 load from LASERBALL_SIMULATION
    else if (fSimulationMode == 0){

        // 1) Wavelength
        // Which dye to simulate - identify by name:
        try { fDyeName = lLASERBALL->GetS("dye_name"); }
        catch( DBNotFoundError& e ) { Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL dye_name info"); };

        // Now get the wl info specific for this dye
        DBLinkPtr lDYECELL;
        // wavelength info
        try{
            fWl_Mode = lLASERBALL->GetS("wl_mode");
            info <<  "VertexGen_Laserball::Init: Wavelength mode is "<< fWl_Mode <<"\n";
        } catch (DBNotFoundError &e) {
            Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL wl_mode info.");
        };

        if(fWl_Mode=="MONO"){
            try { fWl_Mono = lLASERBALL->GetD("wl_mono"); }
            catch( DBNotFoundError& e ) { Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL wl_mono info"); };
            fwaveAvg = fWl_Mono;
            info <<  "VertexGen_Laserball::Init: Setting mono wavelength to "<< fWl_Mono <<" nm\n";

        } else if(fWl_Mode=="DYENAME"){
            try{
                fDyeName = lLASERBALL->GetS("dye_name");
                info << "VertexGen_Laserball::Init: Setting wavelength distribution to " + fDyeName + "\n";
            } catch (DBNotFoundError &e) {
                Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL dye_name info.");
            }

        } else if(fWl_Mode=="DYECELL"){
            try{
                fDyeCell = lLASERBALL->GetI("dye_cell_nbr");
                info << "VertexGen_Laserball::Init: Setting dye cell nbr to " << fDyeCell<< "\n";
            } catch (DBNotFoundError &e) {
                Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL dye_cell_nbr info");
            }
            if (fDyeCell < 0 || fDyeCell > 5)
                Log::Die(dformat("VertexGen_Laserball::Init: LASERBALL dye_cell_nbr \"%d\" out of range.",fDyeCell));

            try{
                fDyeCellSet = lLASERBALL->GetS("dye_cell_set");
                info << "VertexGen_Laserball::Init: Setting dye cell set to " + fDyeCellSet + "\n";
            } catch (DBNotFoundError &e) {
                Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL dye_cell_set info");
            }
            try{
                lDYECELL = DB::Get()->GetLink("DYECELL", fDyeCellSet);
                info <<  "VertexGen_Laserball::Init: Getting dye cell set info for " << fDyeCellSet << "\n";
            }catch(DBNotFoundError &e) {
                Log::Die("VertexGen_Laserball::Init: can't find ratdb DYECELL info for " + fDyeCellSet);
            }
            try{
                fDyeCellOrder = lDYECELL->GetSArray("dye_order");
                info << "VertexGen_Laserball::Init: Setting dye cell order to "<< fDyeCellOrder << "\n";
            } catch (DBNotFoundError &e) {
                Log::Die("VertexGen_Laserball::Init: can't find ratdb DYECELL dye_order info");
            }
            fDyeName = fDyeCellOrder[fDyeCell];

        } else {
            Log::Die("VertexGen_Laserball::Init: LASERBALL wl_mode "+ fWl_Mode + " not recognized");
        }

        // 2) Position and Orientation
        // x,y,z is loaded in PosGen_Laserball
        try{
            fOrientation = lLASERBALL->GetD("orientation");
            info <<  "VertexGen_Laserball::Init: Orientation: "<< fOrientation << "\n";
        } catch (DBNotFoundError &e) {
            Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL_SIMULATION orientation info");
        };
        if (fOrientation < 0.0 || fOrientation > 4.0) {
            Log::Die(dformat("VertexGen_Laserball::Init: LASERBALL orientation  \"%d\" not recognized.",fOrientation));
        }

    } else {
        Log::Die(dformat("PosGen_Laserball: Unexpected simulation mode. Received %d",fSimulationMode));
    }

    // Dye wavelength distribution
    DBLinkPtr lDYEWAVE;
    if(fWl_Mode!="MONO"){

        try{
            lDYEWAVE = DB::Get()->GetLink("LBDYEWAVE", fDyeName);
            info <<  "VertexGen_Laserball::Init: Getting dye wl info for " << fDyeName << "\n";
        }catch(DBNotFoundError &e) {
            Log::Die("VertexGen_Laserball::Init: can't find ratdb LBDYEWAVE info for " + fDyeName);
        }

        // firstly get the wavelengths and store endpoints
        try { fwaveAvg = lDYEWAVE->GetD("dist_wl_avg"); }
        catch( DBNotFoundError& e ) { Log::Die("VertexGen_Laserball::Init: can't find ratdb LBDYEWAVE dist_wl_avg info"); };
        try { fwave = lDYEWAVE->GetFArrayFromD("dist_wl"); }
        catch( DBNotFoundError& e ) { Log::Die("VertexGen_Laserball::Init: can't find ratdb LBDYEWAVE dist_wl info"); };
        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_Laserball::Init: dist_wl[" << iv << "] = " << fwave[iv] << " [" << iv-1 << "] = "  << fwave[iv-1] << "!\n";
                Log::Die("VertexGen_Laserball::Init: dist_wl must increase monatonically. Please check values.");
            }
        }
        // now get the amplitudes
        try { fwaveInt = lDYEWAVE->GetFArrayFromD("dist_wl_intensity"); }
        catch( DBNotFoundError& e ) { Log::Die("VertexGen_Laserball::Init: can't find ratdb LBDYEWAVE 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";
    }

    // 3) Angular distribution
    DBLinkPtr lDYEANGLE;
    try{
        fLBID = lLASERBALL->GetS("lb_id");
        info <<  "VertexGen_Laserball::Init: Laserball ID is "<< fLBID <<"\n";

    } catch (DBNotFoundError &e) {
        Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL lb_id info");
    };
    try{
        fAngleMode = lLASERBALL->GetS("angle_mode");
        info <<  "VertexGen_Laserball::Init: Angular mode is "<< fAngleMode <<"\n";
    } catch (DBNotFoundError &e) {
        Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL angle_mode info");
    };

    if (fAngleMode=="TABLE2D" || fAngleMode=="SINCOEF") {
        try{
            fAngleDistDye = lLASERBALL->GetS("angle_dist_dye");
        } catch (DBNotFoundError &e) {
            Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL angle_dist_dye info");
        };
        if (fAngleDistDye=="DEFAULT") {
            if (fDyeName == "")
                Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL dye_name info, needed for angle_dist_dye ANGLEDEFAULT option"    );
            else {
                info <<  "VertexGen_Laserball::Init: default: angular dist. of "<< fDyeName <<" dye same as for wl\n";
                fAngleDistDye = fDyeName;
            }
        }
        fFullAngleIndex = fAngleDistDye + ":" + fLBID;
        try{
            lDYEANGLE = DB::Get()->GetLink("LBDYEANGLE", fFullAngleIndex);
            info <<  "VertexGen_Laserball::Init: Getting dye angle info for " << fFullAngleIndex << "\n";
        }catch(DBNotFoundError &e) {
            Log::Die("VertexGen_Laserball::Init: can't find ratdb LBDYEANGLE info for " + fFullAngleIndex);
        }
        try{
            fMaskCoef = lDYEANGLE->GetFArrayFromD("lb_mask_coef");
        } catch (DBNotFoundError &e) {
            Log::Die("VertexGen_Laserball::Init: can't find ratdb  LBDYEANGLE lb_mask_coef info");
        };
        if (fAngleMode=="TABLE2D") {
            try{
                fAngIntensity = lDYEANGLE->GetFArrayFromD("lb_ang_intensity");
            } catch (DBNotFoundError &e) {
                Log::Die("VertexGen_Laserball::Init: can't find ratdb  LBDYEANGLE lb_ang_intensity info");
            };
            info <<  "VertexGen_Laserball::Init: Setting angle_mode to TABLE2D (default)\n";
        }
        if (fAngleMode=="SINCOEF") {
            try{
                fAngSinCoef = lDYEANGLE->GetFArrayFromD("lb_ang_sincoef");
            } catch (DBNotFoundError &e) {
                Log::Die("VertexGen_Laserball::Init: can't find ratdb  LBDYEANGLE lb_ang_sincoef info");
            };
            info <<  "VertexGen_Laserball::Init: Setting angle_mode to SINCOEF\n";
        }

        fLBD = new TF2("myfun",this,&VertexGen_Laserball::funlbdist,-1,1,0,twopi,1,"VertexGen_Laserball","funlbdist");
        //fLBD->SetNpx(60);
        info <<  "VertexGen_Laserball::Init: TF2 Npx= "<<fLBD->GetNpx()<< " Npy= "<< fLBD->GetNpy()<<"\n";
        if (fAngleMode=="SINCOEF") fLBD->SetParameter(0,1);
        else fLBD->SetParameter(0,0);
        fLBDMaximum = fLBD->GetMaximum();
        info <<  "VertexGen_Laserball::Init: TF2 Maximum= "<<fLBDMaximum<<"\n";
        if ( fLBDMaximum > 2.0 ){
          info <<  "Above maximum could be errorneous, setting to new maximum value of 1.0\n";
          fLBDMaximum = 1.0;
        }
    } else if (fAngleMode=="FLAT") {
        info <<  "VertexGen_Laserball::Init: Setting angle_mode to FLAT\n";
        fFullAngleIndex = fAngleMode;
    } else {
        Log::Die("VertexGen_Laserball::Init: LASERBALL angle_mode "+ fAngleMode + " not recognized");
    }

    // 4) Intensity
    try {
        fPhotonsPerEvent = int(double(lLASERBALL->GetI("intensity"))/photon_thin);
        debug << "VertexGen_Laserball::Init: Intensity = " << fPhotonsPerEvent << "\n";
    } catch (DBNotFoundError &e) { Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL intensity info"); };
    // Optional setting
    try {
        string temp = lLASERBALL->GetS("pulse_mode");
        if(temp == "POISSON"||temp == "poisson"){
            fPoisson = true;
            info <<  "VertexGen_Laserball::Init: Setting pulse mode to poisson\n";
        }else if(temp == "fixed" || temp == "FIXED"){
            fPoisson = false;
            info <<  "VertexGen_Laserball::Init: Setting pulse mode to fixed\n";
        }else{
            Log::Die(dformat("VertexGen_Laserball::Init: LASERBALL.pulse_mode = \"%s\" is invalid.",
            temp.c_str()));
        }
    } catch (DBNotFoundError &e) { Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL pulse_mode info"); };

    // 5) Timing distribution
    try{
        fTimeMode = lLASERBALL->GetS("time_mode");
        info <<  "VertexGen_Laserball::Init: Time mode is "<< fTimeMode <<"\n";
    } catch (DBNotFoundError &e) {
        Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL time_mode info");
    };
    try{
        fTimeParams = lLASERBALL->GetFArrayFromD("time_param");
    } catch (DBNotFoundError &e) {
        Log::Die("VertexGen_Laserball::Init: can't find ratdb LASERBALL time_param info");
    };

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

}



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



void VertexGen_Laserball::GeneratePrimaryVertex(G4Event *event,
                                                G4ThreeVector &fPos0,
                                                G4double fTime0)
{
    // How many photons?
    int nphot = 0;
    if(fPoisson){
        nphot = (G4int)( CLHEP::RandPoisson::shoot(fPhotonsPerEvent) );
    }else{
        nphot = fPhotonsPerEvent;
    }

    // Record in extra G4Event information
    EventInfo *exinfo = dynamic_cast<EventInfo*>(event->GetUserInformation());
    //DS::Calib *calib = Extensible<DS::Calib>::New(); // new for dbreview branch
    DS::Calib *calib = new DS::Calib();
    calib->SetSourceName(fFullAngleIndex);
    calib->SetMode((int) fwaveAvg);
    calib->SetIntensity(nphot);
    calib->SetPos( TVector3(fPos0.x(), fPos0.y(),fPos0.z()) );

    // Set the LB orientation ID and vector
    float or_x = TMath::Cos( fOrientation * halfpi );
    float or_y = TMath::Sin( fOrientation * halfpi );

    // Check if Simulation Mode is 10, which is the flag for laserball runs without side ropes (external or internal)
    // If that is the case, set the ID in Calib as 10 - default flag for this type of runs
    if(fSimulationMode == 10)
        calib->SetID( (int) 10.0 );
    else{
        calib->SetID( (int)fOrientation );
    }
    calib->SetDir( TVector3(or_x, or_y, 0.0) );

    calib->SetTime(fTimeParams[3]);
    exinfo->SetCalib(calib);

    // Add vertices for each photon
    for (int i=0; i < nphot; i++) {

        // Set time to a GasDev distribution (default) or a single value;
        double time = fTimeParams[3];
        if(!(fTimeMode=="MONO")){
            time = fTimeParams[1] + 1.0;
            while(time < fTimeParams[0] || time > fTimeParams[1]){
                time = GasDev() *   fTimeParams[2] + fTimeParams[3];
            }
        }

        // Set wavelength to input distribution (default) or a single value;
        float wavelength;
        if(fWl_Mode=="MONO"){
            wavelength = fWl_Mono; // nm
            debug << "VertexGen_Laserball::GeneratePrimaryVertex:  LASERBALL wavelength fixed = "
            << wavelength <<"\n";
        }else{
            debug << "VertexGen_Laserball::GeneratePrimaryVertex:  LASERBALL wavelength randomly chosen = ";
            wavelength = SampleRandom(fwave, fwaveInt, fwaveCumMag);
            //wavelength = fRandWl->shoot() * (fMaxWl - fMinWl) + fMinWl; // randomly sample for recording in MC file
            debug  << wavelength <<"\n";
        }
        float energy = hbarc * twopi / (wavelength * nm);
        float momentum = energy; // GEANT uses momentum in same units as energy

        G4ThreeVector dir(0.,0.,-1.0);
        // Angular distribution
        if (fAngleMode=="FLAT") {
            dir = SphericallyUniform3Vector();
        } else {
            //dir = FunctionRandom3Vector();  // Don't use for now: ultra-slow
            dir = BoxMethodRandom3Vector();
        }

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


        G4PrimaryVertex* vertex= new G4PrimaryVertex(fPos0, 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);

    } // loop over photons
}

G4ThreeVector VertexGen_Laserball::SphericallyUniform3Vector()
{

    G4ThreeVector result(0.,0.,0.);
    double aTheta = acos( 1.0 - 2.0 * G4UniformRand() );
    double aSTheta = sin(aTheta);
    double aPhi = twopi * G4UniformRand();

    result.setX( aSTheta * cos(aPhi) );
    result.setY( aSTheta * sin(aPhi) );
    result.setZ( cos(aTheta) );

    return result;
}

G4ThreeVector VertexGen_Laserball::RotatePhiAntiClock(double aPhi, G4ThreeVector aVector3) const
{
  //  Allow for a phi rotation of the laserball
  //  anti-clockwise
  G4ThreeVector result(0.,0.,0.);

  result.setX( aVector3.x() *  cos( aPhi ) -  aVector3.y() *  sin( aPhi ) );
  result.setY( aVector3.x() *  sin( aPhi ) +  aVector3.y() *  cos( aPhi ) );
  result.setZ( aVector3.z() );

  return result;
}

G4ThreeVector VertexGen_Laserball::RotatePhiClock(double aPhi, G4ThreeVector aVector3) const
{
   //  Allow for a phi rotation of the laserball
   //  clockwise
   return RotatePhiAntiClock(-aPhi, aVector3);
}

G4ThreeVector VertexGen_Laserball::FunctionRandom3Vector()
{

    // This is ultra-slow (100 s per each event of 10k photons)
    // There's room to improve on TF2::GetRandom2, maybe we can make our own version
    G4ThreeVector result(0.,0.,0.);
    double aCTheta;
    double aPhi;
    fLBD->GetRandom2(aCTheta,aPhi);
    double aSTheta = sqrt(1 - aCTheta*aCTheta);

    result.setX( aSTheta * cos(aPhi) );
    result.setY( aSTheta * sin(aPhi) );
    result.setZ( aCTheta );

    return result;
}

G4ThreeVector VertexGen_Laserball::BoxMethodRandom3Vector()
{

    G4ThreeVector result(0.,0.,0.);
    G4ThreeVector result_lb_rotated(0.,0.,0.);
    double aCTheta;
    double aPhi;
    double value = 0;
    double rotPhi;
    rotPhi = (double) fOrientation* pi/2.;

    bool accept = false;
    while (!accept){
        result = SphericallyUniform3Vector();
        // Need to go from detector to LB coordinates.
        // If the laserball is in the North position (+90 deg w/r to East) and I'm generating
        // a photon at phi=90, I need to get the intensity value for phi_LB = 0
        // so we need to rotate by -90, i.e. clockwise
        result_lb_rotated = RotatePhiClock(rotPhi,result);
        aCTheta = cos(result_lb_rotated.theta());
        aPhi    = result_lb_rotated.phi();
        // Get the intensity from the ratdb distributions, that are in LB coordinates
        value  = fLBD->Eval(aCTheta,aPhi)/(1.01*fLBDMaximum);
        if (G4UniformRand() < value) accept = true;
    }

    return result;
}


double VertexGen_Laserball::funlbdist(double *x, double *par)
{

    double cTheta = x[0];
    double phi = x[1];
    double fOTheta = 1.0 + cTheta;
    int nMaskPars = (int)fMaskCoef.size();
    double aMask = fMaskCoef[0];
    // Number of mask parameters is (degree+1), so we only
    // want to sum with powers between 1 and (nMaskPars-1).
    for ( int iPar = 1; iPar < nMaskPars; iPar++ ){
      aMask += fMaskCoef[iPar] * pow(fOTheta,iPar);
    }

    int nBinsCT     = 12;
    int nBinsPhi    = 36;
    double binWidthCT = 2.0 /(double)(nBinsCT);
    double binWidthPhi = twopi /(double)(nBinsPhi);
    double aIntensity,aAmplitude_ac,aAmplitude_dc,aFrequency,aPhase;

    int binct0 = (int) ((cTheta+1.0)/binWidthCT);
    if (par[0] > 0.5){
        // sinusoidal model can use 24 or 12 theta bins, should be defined by half the size of the parameters list
        // range from -1.0 to 1.0 (=2.0).
        const double cosThBins = fAngSinCoef.size()/2.0;
        const double cosThRange = 2.0;
        binct0 = (int) ((cTheta+1.0)/(cosThRange/cosThBins));
        // sin wave parametrization
        aAmplitude_ac = fAngSinCoef[binct0*2];
        aAmplitude_dc = 1.0;
        aFrequency    = 1.0;
        aPhase        = fAngSinCoef[binct0*2 +1];
        aIntensity    = aAmplitude_dc + aAmplitude_ac*sin(aFrequency*phi + aPhase);
    } else {
        // 2d table parametrization. interpolate in both theta and phi (to do)
        if (phi<= -pi)     phi += twopi;
        if (phi>   pi)     phi -= twopi;

        // Binning in ratdb follows -pi, to pi convention. i.e, bin 0 is -pi, bin 35 is pi
        // So we need to add pi before converting to bin number
        int binphi0 = (int) ((phi+pi)/binWidthPhi);
        aIntensity = fAngIntensity[binct0*nBinsPhi + binphi0];
    }

    aIntensity *= aMask;
    return aIntensity;

}

double VertexGen_Laserball::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;
}

float VertexGen_Laserball::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_Laserball::SetState(G4String newValues)
{

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

}


G4String VertexGen_Laserball::GetState()
{
    stringstream result;
    result << fDyeName << ":" << fLBID;
    return result.str();
}




} // namespace RAT