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

#include <RAT/Coincidence_Gen.hh>
#include <RAT/CoincidenceMessenger.hh>

#include <RAT/Log.hh>

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

#include <G4Event.hh>
#include <G4PrimaryVertex.hh>
#include <G4ThreeVector.hh>
#include <G4PhysicalConstants.hh>

#include <Randomize.hh>

#include <vector>
#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>
#include <cstring>
using std::vector;
using std::string;

#define G4std  std

//#undef DEBUG
using namespace CLHEP;

namespace RAT {
    Coincidence_Gen::Coincidence_Gen() : GLG4Gen("coincidence"),fStateStr("") {

        // Create a messenger to allow the user to change some parameters.
        fMessenger = new CoincidenceMessenger(this);

        // As in the generator, use a default time generator if the
        // user does not supply one.
        fTimeGen     = new GLG4TimeGen_Poisson();
        // also use default time window
        fTimeWindow = 400;        // 400ns
        fLoEnergyLimit     = 0;         // 0 MeV
        fHiEnergyLimit     = 99999;     // 99999 MeV
        fNExtra         = 0;        // number of additional interactions
        //Use random timing as default
        fIsExpTiming  = false;
        fIsFixedTiming = false;
        fIsGaussianTiming = false;
        fIsQuadraticTiming = false;
    }

/*--------------------------------------------------------------*/
Coincidence_Gen::~Coincidence_Gen() {
  if ( fMessenger != 0 ){
      delete fMessenger;
      fMessenger = 0;
    }
#ifdef RAT_GEN_DEBUG
  info << "Coincidence_Gen::~Coincidence_Gen : Destroying aux generators..." << newline;
  info << "                                  : Expecting " << fPosGenExtra.size() << newline;
#endif

  for (size_t i = 0; i < fPosGenExtra.size(); ++i) {
    delete fPosGenExtra.at(i);
  }
  fPosGenExtra.clear();
#ifdef RAT_GEN_DEBUG
  info << "Coincidence_Gen::~Coincidence_Gen : Destroying vertex aux generators..." << newline;
  info << "                                  : Expecting " << fVertexGenExtra.size() << newline;
#endif
  for (size_t i = 0; i < fVertexGenExtra.size(); ++i) {
    delete fVertexGenExtra.at(i);
  }
  fVertexGenExtra.clear();
#ifdef RAT_GEN_DEBUG
  info << "Coincidence_Gen::~Coincidence_Gen : Destroying state str..." << newline;
  info << "                                  : Expected " <<  fStateStrExtra .size() << newline;
#endif

  fStateStrExtra.clear();
}

/*--------------------------------------------------------------*/
    void Coincidence_Gen::GenerateEvent(G4Event* event)
    {
        G4ThreeVector     pos;
        // We will need to keep track of how much energy has been generated if the
        // total event energy is to be limited
        G4double        KEgen=0;

        // generate times for all the events
        vector<double> times;
        vector<double> MaxE;
        vector<double> MinE;
        if(fIsExpTiming || fIsFixedTiming || fIsGaussianTiming || fIsQuadraticTiming){
            //user has set a time distribution other than default
            //first vertex is at time 0
            times.push_back(0);
            double lastt = 0;           //store time of previous vertex
            for(int iEx=0; iEx<fNExtra; iEx++){
                if(fIsExpTiming){
                    lastt +=CLHEP::RandExponential::shoot(fExponent[iEx]);
                }else if(fIsGaussianTiming){
                    lastt +=CLHEP::RandGauss::shoot(fGaussMean[iEx],fGaussStdDev[iEx]);
                }else if(fIsQuadraticTiming){
                    //Transform the user-provided quadratic probability density function
                    //to the quantile function (inverse cumulative density function)
                    //Use this quantile to transform flat random numbers
                    //to random numbers weighted by the distribution

                    //calculate discriminant of the quadratic
                    double discriminant = (fQuadB[iEx]*fQuadB[iEx]);
                    discriminant -= (fQuadA[iEx]*fQuadC[iEx]*4.);
                    //lowlimit and highlimit are the 2 roots of the quadratic
                    //pdf is only defined where quadratic is positive
                    double lowlimit, highlimit, d, norm;
                    highlimit=(-fQuadB[iEx]-sqrt(discriminant))/(2.*fQuadA[iEx]);
                    lowlimit=(-fQuadB[iEx]+sqrt(discriminant))/(2.*fQuadA[iEx]);
                    d = -(fQuadA[iEx]/3.)*pow(lowlimit,3);
                    d -= (fQuadB[iEx]/2.)*lowlimit*lowlimit;
                    d -= fQuadC[iEx]*lowlimit;

                    //normalize total probability density to 1
                    //change to coefficients for the CDF (divide a by 3, b by 2)
                    norm = (fQuadA[iEx]/3.)*pow(highlimit, 3);
                    norm += (fQuadB[iEx]/2.)*(highlimit*highlimit);
                    norm += fQuadC[iEx]*highlimit + d;
                    fQuadA[iEx] = fQuadA[iEx]/(3.*norm);
                    fQuadB[iEx] = fQuadB[iEx]/(2.*norm);
                    fQuadC[iEx] = fQuadC[iEx]/norm;
                    d = d/norm;

                    //Solve CDF for quadratic-weighted time
                    //use trigonometric method for solving cubics
                    //begin by transforming to a reduced cubic
                    //p and q are coefficients for reduced cubic t^3 + pt + q = 0
                    double p, q, rand, nextt;
                    rand = CLHEP::RandFlat::shoot();
                    p = (3.*fQuadA[iEx]*fQuadC[iEx]) - (fQuadB[iEx]*fQuadB[iEx]);
                    p /= (3.*fQuadA[iEx]*fQuadA[iEx]);
                    q = (2.*pow(fQuadB[iEx],3))-(9.*fQuadA[iEx]*fQuadB[iEx]*fQuadC[iEx]);
                    q += (27.*fQuadA[iEx]*fQuadA[iEx]*(d-rand));
                    q /= (27.*pow(fQuadA[iEx],3));
                    //Use formula for second root of a cubic with 3 real roots
                    nextt = acos(q*sqrt(27.)/(2.*p*sqrt(-p))) - twopi;
                    nextt = 2.*sqrt(-p/3.)*cos((1./3.)*nextt) - fQuadB[iEx]/(3.*fQuadA[iEx]);
                    lastt += nextt;
                    //return coefficients to their previous values
                    fQuadA[iEx] = fQuadA[iEx]*(3.*norm);
                    fQuadB[iEx] = fQuadB[iEx]*(2.*norm);
                    fQuadC[iEx] = fQuadC[iEx]*norm;


                }else{
                    lastt+=fFixed[iEx];
                }
                times.push_back(lastt);
            }
        }else{
             // each event has random time in event window
            // the first one will be the start of the event so  keep track of that
            double firsttime = fTimeWindow; // set to latest possible value to begin
            for(int iEx=0; iEx<=fNExtra; ++iEx){
                double t = G4UniformRand()*fTimeWindow;
                times.push_back(t);    // for the extra interactions
                if(t<firsttime)firsttime=t;
            }
            // now subtract this first time from every time
            for(int iEx=0; iEx<=fNExtra; ++iEx){
                times[iEx] -= firsttime;
            }
        }
        // now access the maxima and minima energy each extra generator can provide
        for(int iEx=0; iEx<fNExtra; ++iEx){
            MaxE.push_back(fVertexGenExtra[iEx]->EMaximum());
            MinE.push_back(fVertexGenExtra[iEx]->EMinimum());
        }
        // But if more than one extra generator, we actually want a cumulative total of the other energy limits
        float sumLo = 0;
        float sumHi = 0;
        for(int iEx=fNExtra-1; iEx>=0; --iEx){
            sumLo += MinE[iEx];
            sumHi += MaxE[iEx];
            MinE[iEx] = sumLo;
            MaxE[iEx] = sumHi;
        }

        // A test vertex and particle (used to extract event energy)
        G4PrimaryVertex* vtemp=0;
        G4PrimaryParticle* PP;
        // A temporary event
        G4Event* etemp;
        // first position
        fPosGen->GeneratePosition(pos);
        // the vertex generators add a primary vertex to event - ready to addirst one
        // if we need to limit the upper energy for this vertex, do so
        float lo,hi;
        if(fLoEnergyLimit>0||fHiEnergyLimit<99999){
            lo = fLoEnergyLimit-sumHi;
            hi = fHiEnergyLimit-sumLo;
            if(lo<0)lo = 0;
            if(hi<0){
                warn << "Coincidence_Gen: cannot achieve selected energy range with these generators.\n";
            }
            if(fVertexGen->ELimitable()){
                fVertexGen->LimitEnergies(lo,hi);
                fVertexGen->GeneratePrimaryVertex(event, pos, times[0]);
            }else{
                // we will have to apply trial and error so generate temporary event and test it
                G4double testE = -999*MeV;
                int count=1;
                while(count<100&&!(testE/MeV>lo&&testE/MeV<hi)){
                    etemp = new G4Event();        // need to reset the temporary event
                    fVertexGen->GeneratePrimaryVertex(etemp, pos, times[0]);
                    vtemp = etemp->GetPrimaryVertex();
                    testE = 0;
                    for(int iP=0;iP<vtemp->GetNumberOfParticle(); ++iP){
                        PP = vtemp->GetPrimary(iP);
                        testE += sqrt(PP->GetMomentum().mag2()+pow(PP->GetMass(),2))-PP->GetMass();
                        // NOTE !!! problem if the particle is short lived
                        // how do I access decay particle energy???
                    }
                    count++;
                }
                if(count==100){
                    // Jump out of what could have been infinite loop
                    warn << "No luck finding correct energy ! \n";
                }
                // got a good vertex so add it to the real event
                event->AddPrimaryVertex(vtemp);
            }
            // now access event, primary vertex and extract energy
            vtemp = event->GetPrimaryVertex();
            for(int iP=0;iP<vtemp->GetNumberOfParticle(); ++iP){
                PP = vtemp->GetPrimary(iP);
                KEgen += sqrt(PP->GetMomentum().mag2()+pow(PP->GetMass(),2))-PP->GetMass();
            }
        }else{
            // No limits on the energy so just add the vertex
            fVertexGen->GeneratePrimaryVertex(event, pos, times[0]);
        }

        // Now loop through extra vertices
        for(int iEx=0; iEx<fNExtra; ++iEx){
            fPosGenExtra[iEx] -> GeneratePosition(pos);
            if(fLoEnergyLimit>0||fHiEnergyLimit<99999){
                if(iEx<fNExtra-1){
                    lo = fLoEnergyLimit-MaxE[iEx+1]-KEgen/MeV;
                    hi = fHiEnergyLimit-MinE[iEx+1]-KEgen/MeV;
                }else{
                    lo = fLoEnergyLimit-KEgen/MeV;
                    hi = fHiEnergyLimit-KEgen/MeV;
                }
                if(fVertexGenExtra[iEx]->ELimitable()){
                    // Set limits on energy and add the vertex
                    fVertexGenExtra[iEx]->LimitEnergies(lo,hi);
                    fVertexGenExtra[iEx]->GeneratePrimaryVertex(event, pos, times[iEx+1]);
                }else{
                    // we will have to apply trial and error so generate temporary event and test it
                    G4double testE = -999*MeV;
                    int count=1;
                    while(count<100&&!(testE/MeV>lo&&testE/MeV<hi)){
                        etemp = new G4Event();        // need to reset the temporary event
                        fVertexGenExtra[iEx]->GeneratePrimaryVertex(etemp, pos, times[iEx+1]);
                        vtemp = etemp->GetPrimaryVertex();
                        testE = 0;
                        for(int iP=0;iP<vtemp->GetNumberOfParticle(); ++iP){
                            PP = vtemp->GetPrimary(iP);
                            testE += sqrt(PP->GetMomentum().mag2()+pow(PP->GetMass(),2))-PP->GetMass();
                        }
                        count++;
                    }
                    if(count==100||count==1){
                        // Jump out of what could have been infinite loop
                        warn << "No luck finding correct energy ! \n";
                    }
                    // got a good vertex so add it to the real event
                    event->AddPrimaryVertex(vtemp);
                }
                // and increment our energy sum (if not last event)
                if(iEx<(fNExtra-1)){
                    vtemp = event->GetPrimaryVertex(iEx+1);
                    for(int iP=0;iP<vtemp->GetNumberOfParticle(); ++iP){
                        PP = vtemp->GetPrimary(iP);
                        KEgen += sqrt(PP->GetMomentum().mag2()+pow(PP->GetMass(),2))-PP->GetMass();
                    }
                }
            }else{
                // No limits on the energy so just add the vertex
                fVertexGenExtra[iEx]->GeneratePrimaryVertex(event, pos, times[iEx+1]);
            }
        }
    }

/*--------------------------------------------------------------*/
    void Coincidence_Gen::ResetTime(double offset)
    {
        // This applies to the overall event time
        double eventTime = fTimeGen->GenerateEventTime();
        fNextTime = eventTime + offset;
    }

/*--------------------------------------------------------------*/
  void Coincidence_Gen::SetState(G4String state)
  {
    // arguments are V1:P1:T
    // V1 = first vertex type, required
    // P1 = first vertex position, required
    // T  = overall time generator, optional, default = poisson

    state = util_strip_default(state);

    vector<string> parts = util_split(state, ":");

    try {
      switch (parts.size()) {
      case 3:
        // last is optional time generator
        delete fTimeGen; fTimeGen = 0; // In case of exception in next line
        fTimeGen = RAT::GlobalFactory<GLG4TimeGen>::New(parts[2]);
        //JW - is that a bug? fTimeGen = RAT::GlobalFactory<GLG4TimeGen>::New(parts[4]);
      case 2:
        info << "adding first interaction " << fNExtra << " " << parts[0] << " "<<
          parts[1] << "\n";

        delete fPosGen; fPosGen = 0;
        fPosGen = RAT::GlobalFactory<GLG4PosGen>::New(parts[1]);
        delete fVertexGen; fVertexGen = 0;
        fVertexGen = RAT::GlobalFactory<GLG4VertexGen>::New(parts[0]);
        break;
      default:
        G4Exception(__FILE__, "Invalid Parameter", FatalException, ("Coincidence generator syntax error: "+state).c_str());
        break;
      }
      fStateStr = state; // Save for later call to GetState()
    } catch (RAT::FactoryUnknownID &unknown) {
      warn << "Unknown generator \"" << unknown.id << "\"\n";
    }
  }

  G4String Coincidence_Gen::GetState() const
  {
    return fStateStr;
  }

/*--------------------------------------------------------------*/
    void Coincidence_Gen::SetEnergyRange(G4String newValues)
    {
        newValues = util_strip_default(newValues);
        if (newValues.length() == 0) {
            // Print help and current state
            info     << " Current limits on Coincidence_Gen generated energy range = "
                       << fLoEnergyLimit << " - " << fHiEnergyLimit << " MeV\n"
                    << " To change use syntax:  Elo  Ehi \n"
                    << " Elo = lower energy limit in MeV, Ehi = upper energy limit in MeV"
                    << " !!!! But note that this option will be very slow for vertex generators whose energy cannot be limited !!! \n";
        }
         G4std::istringstream is(newValues.c_str());
        double Elo, Ehi;
        is >> Elo >> Ehi;
        if(Elo<Ehi){
            fLoEnergyLimit = Elo;
            fHiEnergyLimit = Ehi;
            info  << "Coincidence_Gen: Restricting generated energy range to " << Elo << " - " << Ehi << " MeV\n"
                    << " !!!! But note that this option will be very slow for vertex generators whose energy cannot be limited !!!! "
                    << " !!!! And does not work when event energy comes from radioactive decay !!!! \n";
        }else{
            warn << "Coincidence_Gen error: Elo must be less than Ehi!\n";
        }
    }

/*--------------------------------------------------------------*/
    void Coincidence_Gen::SetTimeState(G4String state)
    {
        if (fTimeGen){
            fTimeGen->SetState(state);
        }else{
            warn << "Coincidence_Gen error: Cannot set time state, no time generator selected\n";
        }
    }

/*--------------------------------------------------------------*/
    G4String Coincidence_Gen::GetTimeState() const
    {
        if (fTimeGen){
            return fTimeGen->GetState();
        }else{
            return G4String("Coincidence_Gen error: no time generator selected");
        }
    }

/*--------------------------------------------------------------*/
    void Coincidence_Gen::SetPosState(G4String state)
    {
        // Set the position generator for the first interaction type
        if (fPosGen){
            fPosGen->SetState(state);
            info << "Setting first interaction " << fNExtra << ", position type to " << state << "\n";
        }else{
            warn << "Coincidence_Gen error: Cannot set position state, no position generator selected\n";
        }
    }

/*--------------------------------------------------------------*/
    G4String Coincidence_Gen::GetPosState() const
    {
        // Get the position generator for the first interaction type
        if (fPosGen){
            return fPosGen->GetState();
        }else{
            return G4String("Coincidence_Gen error: no pos generator selected");
        }
    }

/*--------------------------------------------------------------*/
    void Coincidence_Gen::SetVertexState(G4String state)
    {

      if (fNExtra>0){                                    // Set the vertex generator for extra interactions
        fVertexGenExtra[fNExtra-1]->SetState(state);
        info << "Setting extra interaction " << fNExtra << ", vertex type to " << state << "\n";
      } else if (fVertexGen){                              // Set the vertex generator for the first interaction
          fVertexGen->SetState(state);
          info << "Setting first interaction " << fNExtra << ", vertex type to " << state << "\n";
        }else{
          warn << "Coincidence_Gen error: Cannot set vertex state, no vertex generator selected\n";
        }
    }

/*--------------------------------------------------------------*/
    G4String Coincidence_Gen::GetVertexState() const
    {
        // Get the name of the vertex generator for the first interaction
        if (fVertexGen){
            return fVertexGen->GetState();
        }else{
            return G4String("Coincidence_Gen error: no vertex generator selected");
        }
    }

/*--------------------------------------------------------------*/
    void Coincidence_Gen::AddExtra(G4String state){
        // Add an extra vertex to the same event
        // arguments are V:P
        // V = vertex type, required
        // P = vertex position, required

        state = util_strip_default(state);

        vector<string> parts = util_split(state, ":");

        try {
            switch (parts.size()) {
                case 2:
                    info << "adding new interaction " << fNExtra+1 << " " << parts[0] << " "<<
                    parts[1] << "\n";

                    fPosGenExtra.push_back(RAT::GlobalFactory<GLG4PosGen>::New(parts[1]));
                    fVertexGenExtra.push_back(RAT::GlobalFactory<GLG4VertexGen>::New(parts[0]));
                    // successfully added so increment number of extra states
                    fNExtra++;
                    //Add elements to vectors with parameters for various timing distributions
                    fQuadA.push_back(0);
                    fQuadB.push_back(0);
                    fQuadC.push_back(0);
                    fGaussMean.push_back(0);
                    fGaussStdDev.push_back(0);
                    fExponent.push_back(0);
                    fFixed.push_back(0);
                    break;
                default:
                        G4Exception(__FILE__, "Invalid Parameter", FatalException, ("Coincidence generator syntax error: "+state).c_str());
                        fStateStrExtra.pop_back();
                    break;
            }
            fStateStrExtra.push_back(state); // Save for later call to GetState()

        } catch (RAT::FactoryUnknownID &unknown) {
            warn << "Unknown generator \"" << unknown.id << "\"\n";
        }
    }

/*--------------------------------------------------------------*/
    G4String Coincidence_Gen::GetExtraState(int nint) const
    {
        if (nint<=fNExtra){    // check we have added an extra interaction
            return *fStateStrExtra[nint-1];
        }else{
            return G4String("Coincidence_Gen error: that extra interaction has not been selected");
        }
    }

/*--------------------------------------------------------------*/
    void Coincidence_Gen::SetExtraPosState(G4String state)
    {
        // set the position generator name for the most recently added extra interaction
        if (fNExtra>0){    // check we have added an extra interaction
            fPosGenExtra[fNExtra-1]->SetState(state);
            info << "Setting extra interaction " << fNExtra << ", position type to " << state << "\n";
        }else{
            warn << "Coincidence_Gen error: Cannot set extra position state, no extra interaction selected\n";
        }
    }

/*--------------------------------------------------------------*/
    G4String Coincidence_Gen::GetExtraPosState(int nint) const
    {
        // return the position generator name for extra state nint
        if (nint<=fNExtra){    // check we have added an extra interaction
            return fPosGenExtra[nint-1]->GetState();
        }else{
            return G4String("Coincidence_Gen error: that extra interaction has not been selected");
        }
    }

/*--------------------------------------------------------------*/
    void Coincidence_Gen::SetExtraVertexState(G4String state)
    {
        // Set the vertex generator fo the most recently added in interaction
        if (fNExtra>0){
            fVertexGenExtra[fNExtra-1]->SetState(state);
            info << "Setting extra interaction " << fNExtra << ", vertex type to " << state << "\n";
        }else{
            warn << "Coincidence_Gen error: Cannot set extra vertex state, no extra interaction selected\n";
        }
    }

/*--------------------------------------------------------------*/
    G4String Coincidence_Gen::GetExtraVertexState(int nint) const
    {
        if (nint<=fNExtra){
            return fVertexGenExtra[nint-1]->GetState();
        }else{
            return G4String("Coincidence_Gen error: that extra interaction has not been selected");
        }
    }

/*--------------------------------------------------------------*/
    void Coincidence_Gen::SetTimeWindow(double window)
    {
        //set the time window over which 2 events are coincident
        //use random timing within the window

        fTimeWindow = window;
        fIsQuadraticTiming = false;
        fIsGaussianTiming = false;
        fIsFixedTiming = false;
        fIsExpTiming = false; // if user sets multiple time distributions, most recent distribution will be used

    }

/*--------------------------------------------------------------*/
    void Coincidence_Gen::SetExponentials(G4String newValues)
    {
        /* Set exponential time constants to separate decays.
        * Default is to give each vertex a random time in the event window but if this option is called
        * User must set a time constant for each added vertex.
        * First vertex is at t - 0, the rest have times selected from exponentials
        * each exponential applies to the time since last vertex */

        fIsQuadraticTiming = false;
        fIsGaussianTiming = false;
        fIsFixedTiming = false;
        fIsExpTiming = true; // if user sets multiple time distributions, most recent distribution will be used
        G4std::istringstream is(newValues.c_str());
        for(int i=0; i<fNExtra; ++i){
            if(is.good()){
                is >> fExponent[i];
                info << "Coincidence_Gen: Setting extra interaction " << i+1
                       << " exponential time constant to " <<  fExponent[i] << "\n";
            }else{
                warn << "Coincidence_Gen error: "
                       << "Exponential timing selected, but insufficient time constants provided \n";
            }
        }
    }


/*--------------------------------------------------------------*/
    void Coincidence_Gen::SetFixedTimes(G4String newValues)
    {
        // Set fixed timing distribution for subsequent decays within coincidence
        // Default is to give decays random time within the window, but if this option is called
        // the user must set a time between decays for each added vertex
        // First vertex is at t = 0, the rest have times defined by the user
        // Each argument gives the time to the next decay in ns

        fIsQuadraticTiming = false;
        fIsGaussianTiming = false;
        fIsFixedTiming = true;
        fIsExpTiming = false; // if user sets multiple time distributions, most recent distribution will be used
        G4std::istringstream is(newValues.c_str());
        for(int i=0; i<fNExtra; i++){
            if(is.good()){
                is >> fFixed[i];
                info << "Coincidence_Gen: Setting extra interaction " << i+1
                        << " fixed time step to " << fFixed[i] << "\n";
            }else{
                warn << "Coincidence_Gen error: "
                        << "Fixed timing selected, but insufficient time constants provided \n";
            }
        }
    }

/*--------------------------------------------------------------*/
    void Coincidence_Gen::SetGaussians(G4String newValues)
    {
        // Set Gaussian timing for separate decays
        // Default is to give decays random time within the window, but if this option is called
        // the user must set a mean time between decays for each added vertex
        // and a standard deviation from the mean
        // First vertex is at t = 0, the rest have times generated by a random normal number

        fIsQuadraticTiming = false;
        fIsGaussianTiming = true;
        fIsFixedTiming = false;
        fIsExpTiming = false; // if user sets multiple time distributions, most recent distribution will be used
        G4std::istringstream is(newValues.c_str());
        for(int i=0; i<fNExtra; i++){
            if(is.good()){
                is >> fGaussMean[i];
                info << "Coincidence_Gen: Setting extra interaction " << i+1
                        << " Gaussian mean time step to " << fGaussMean[i] << "\n";
            }else{
                warn << "Coincidence_Gen error: "
                        << "Gaussian timing selected, but insufficient time constants provided \n";
            }
            if(is.good()){
                is >> fGaussStdDev[i];
                info << "Coincidence_Gen: Setting extra interaction " << i+1
                        << " Gaussian standard deviation to " << fGaussStdDev[i] << "\n";
            }else{
                warn << "Coincidence_Gen error: "
                        << "Gaussian timing selected, but insufficient time constants provided \n";
            }
        }
    }
/*--------------------------------------------------------------*/
    void Coincidence_Gen::SetQuadratic(G4String newValues)
    {
        // Set quadratic timing distribution for separate decays
        // Default is to give decays random time within the window, but if this option is called
        // the user must give coefficients for a quadratic probability distribution function
        // Input should be of the form a1 b1 c1 a2 b2 c2 ...
        // This will form probability distribution function a1t^2 + b1t + c1 for the first decay, similar for subsequent decays
        // First vertex is at t=0

        fIsQuadraticTiming = true;
        fIsGaussianTiming = false;
        fIsFixedTiming = false;
        fIsExpTiming = false; // if user sets multiple time distributions, most recent distribution will be used
        G4std::istringstream is(newValues.c_str());
        for(int i=0; i<fNExtra; i++){
            if(is.good()){
                is >> fQuadA[i];
                info << "Coincidence_Gen: Setting extra interaction " << i+1
                        << " quadratic probability density coefficient A to " << fQuadA[i] << "\n";
            }else{
                warn << "Coincidence_Gen error: "
                        << "Quadratic timing selected, but insufficient time constants provided "
                        << "\n";
            }
            if(is.good()){
                is >> fQuadB[i];
                info << "Coincidence_Gen: Setting extra interaction " << i+1
                        << " quadratic probability density coefficient B to " << fQuadB[i] << "\n";
            }else{
                warn << "Coincidence_Gen error: "
                        << "Quadratic timing selected, but insufficient time constants provided "
                        << "\n";
            }
            if(is.good()){
                is >> fQuadC[i];
                info << "Coincidence_Gen: Setting extra interaction " << i+1
                        << " quadratic probability density coefficient C to " << fQuadC[i] << "\n";
            }else{
               warn << "Coincidence_Gen error: "
                        << "Quadratic timing selected, but insufficient time constants provided\n";
            }
            if((fQuadA[i] >= 0) || ((fQuadB[i]*fQuadB[i]) - (fQuadA[i]*fQuadC[i]*4.) < 0)){
                warn << "Coincidence_Gen error: "
                        <<"Invalid quadratic probability density function entered\n";
            }
        }
    }

    /*--------------------------------------------------------------*/
      void Coincidence_Gen::BeginOfRun() {
#ifdef RAT_GEN_DEBUG
        info << "Coincidence_Gen::BeginOfRun : Initializing aux position generators..." << newline;
#endif
        GLG4Gen::BeginOfRun(); // Initialize the base generators
        size_t i = 0;
        for (i = 0; i < fPosGenExtra.size(); ++i) {
#ifdef RAT_GEN_DEBUG
        info << "                            : " << i << newline;
#endif
          fPosGenExtra.at(i)->BeginOfRun();
        }
#ifdef RAT_GEN_DEBUG
        info << "Coincidence_Gen::BeginOfRun : Initializing aux vertex generators..." << newline;
#endif
        for (i = 0; i < fVertexGenExtra.size(); ++i) {
#ifdef RAT_GEN_DEBUG
        info << "                            : " << i << newline;
#endif
          fVertexGenExtra.at(i)->BeginOfRun();
        }
      }

      void Coincidence_Gen::EndOfRun() {
        GLG4Gen::EndOfRun(); // End the base generators
#ifdef RAT_GEN_DEBUG
        info << "Coincidence_Gen::EndOfRun : Ending aux position generators..." << newline;
#endif
        size_t i = 0;
        for (i = 0; i < fPosGenExtra.size(); ++i) {
          fPosGenExtra.at(i)->EndOfRun();
        }
#ifdef RAT_GEN_DEBUG
        info << "Coincidence_Gen::EndOfRun : Ending aux vertex generators..." << newline;
#endif
        for (i = 0; i < fVertexGenExtra.size(); ++i) {
          fVertexGenExtra.at(i)->EndOfRun();
        }
      }
}