// Gen_MultiCombo.cc
// Contact person: Phil Jones <p.g.jones@qmul.ac.uk>
// See Gen_MultiCombo.hh for more details
//———————————————————————//

#include <RAT/Gen_MultiCombo.hh>

#include <Randomize.hh>
#include <G4ThreeVector.hh>
#include <RAT/StringUtil.hh>
#include <RAT/string_utilities.hpp>
#include <RAT/GLG4TimeGen.hh>
#include <RAT/GLG4VertexGen.hh>
#include <RAT/GLG4PosGen.hh>
#include <RAT/Factory.hh>
#include <RAT/Log.hh>
#include <sstream>

using std::string;
using std::vector;

namespace RAT {
Gen_MultiCombo::Gen_MultiCombo() : GLG4Gen("MultiCombo"),
    fGenPairs(vector<GenPair * >()), fStateStr(""), fPairStateStr(""),
    fVertexGenType(""), fPosGenType(""), fWeightSum(0)
{
    // Default time generator, in case none is specified during SetState
    fTimeGen = new GLG4TimeGen_Poisson();
}

Gen_MultiCombo::~Gen_MultiCombo()
{
    //Only the generators are newed, delete all of them
   // delete fTimeGen;
    for (size_t i = 0; i < fGenPairs.size(); ++i) {
      delete fGenPairs.at(i);
    }
}

//Constructor for a new GenPair, the struct packaging a Vertex Generator,
//Position Generator, and weight
Gen_MultiCombo::GenPair::GenPair(G4String vtxType, G4String vtxState,
        G4String posType, G4String posState, double weight) :
    fWeight(weight)
{
    try {
        fVtxGen =
            GlobalFactory<GLG4VertexGen>::New(vtxType);
        fVtxGen->SetState(vtxState);

        fPosGen =
            GlobalFactory<GLG4PosGen>::New(posType);
        fPosGen->SetState(posState);

    } catch (FactoryUnknownID &unknown) {
        Log::Die("Gen_MultiCombo: Unknown generator \"" + unknown.id +
                "\", terminating.");
    }
}

Gen_MultiCombo::GenPair::~GenPair()
{
    delete fVtxGen;
    delete fPosGen;
}

void Gen_MultiCombo::BeginOfRun() {
  // Call the base BeginofRun to deal with local instances
  GLG4Gen::BeginOfRun();
  for (size_t i = 0; i < fGenPairs.size(); ++i) {
    fGenPairs.at(i)->fVtxGen->BeginOfRun();
    fGenPairs.at(i)->fPosGen->BeginOfRun();
  }

}

void Gen_MultiCombo::EndOfRun() {
  GLG4Gen::EndOfRun();
  for (size_t i = 0; i < fGenPairs.size(); ++i) {
    fGenPairs.at(i)->fVtxGen->EndOfRun();
    fGenPairs.at(i)->fPosGen->EndOfRun();
  }
}

void Gen_MultiCombo::SetState(G4String state)
{
    state = trim(state);
    vector<string> parts = split(state, ":");

    try {
        switch (parts.size()) {
            case 3:
                // last parameter is optional time generator
                // Clear fTimeGen before resetting, in case of exception
                delete fTimeGen; fTimeGen = NULL;
                fTimeGen = RAT::GlobalFactory<GLG4TimeGen>::New(parts[2]);
            case 2:
                fVertexGenType = parts[0];
                fPosGenType = parts[1];
                break;
            default:
                Log::Die("Combo generator syntax error: "+state);
                break;
        }
        // Save for later call to GetState()
        fStateStr = state;

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

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

// wrapper for the Time Generator's own SetState method
void Gen_MultiCombo::SetTimeState(G4String state)
{
    if (fTimeGen)
        fTimeGen->SetState(state);
    else
        warn << "Gen_MultiCombo error: " <<
            "Cannot set time state, no time generator selected"
            << newline;
}

// wrapper for the Time Generator's own GetState method
G4String Gen_MultiCombo::GetTimeState() const
{
    if (fTimeGen)
        return fTimeGen->GetState();
    else
        return
            G4String("Gen_MultiCombo error: no time generator selected");
}

// SetVertexState and SetPosState have identical functionality:
//
// If the string passed is empty (or consists entirely of whitespace),
// display the current state on info,
// and the format of a state string on detail.
//
// Otherwise, use AddGenPair to try and interpret the string as a list
// of triples with which to construct a list of GenPairs, and add those
// GenPairs to fGenPairs
void Gen_MultiCombo::SetVertexState(G4String state)
{
    if (state.length() == 0)
        DisplayFormat();
    else
        AddGenPair(state);
}

// GetVertexState and GetPosState have identical functionality:
// return the StateStr as reconstructed by AddGenPair,
// which should be in a format that can be passed to Set(Vertex|Pos)State
G4String Gen_MultiCombo::GetVertexState() const
{
    return fPairStateStr;
}

// See comment for SetVertexState
void Gen_MultiCombo::SetPosState(G4String state)
{
    if (state.length() == 0)
        DisplayFormat();
    else
        AddGenPair(state);
}

// See comment for GetVertexState
G4String Gen_MultiCombo::GetPosState() const
{
    return fPairStateStr;
}

// GenerateEvents:
//
// If no GenPairs have been added yet, nothing can be done, so warn and exit
//
// Otherwise, choose a random GenPair,
// weighted by the weight member of each, using a flat distribution.
// Then use that GenPair to generate an event.
void Gen_MultiCombo::GenerateEvent(G4Event *event)
{
    if(fGenPairs.empty()){
        warn << "Multicombo generator has no vertex/position pairs "
             << "- unable to generate event!" << newline;
    } else {
        // Create a random number in (0,fWeightSum)
        double rand = CLHEP::RandFlat::shoot(fWeightSum);
        // Now we want to find which weight interval the random number falls
        // into: (0,fWeightSum) = (0,w1)U(w1,w1+w2)U(w1+w2,w1+w2+w3)U...

        vector<GenPair * >::iterator it;
        for(it=fGenPairs.begin(); it<fGenPairs.end(); it++){
            if(rand < (*it)->fWeight){
                // If the random number is less than current GenPair's
                // weight, we've found the right interval, exit the loop
                break;
            } else {
                // Not in this interval, so subtract current weight from
                // random number, so that we can just compare to weight
                // above, rather than an increasing sum.
                rand -= (*it)->fWeight;
            }
        }
        // 'it' now points to the right generator.
        G4ThreeVector pos;
        (*it)->fPosGen->GeneratePosition(pos);
        (*it)->fVtxGen->GeneratePrimaryVertex(event, pos, NextTime());
    }
}

void Gen_MultiCombo::ResetTime(double offset)
{
    fNextTime = fTimeGen->GenerateEventTime() + offset;
}

// Here we try and interpret a string as a list of fields
// with which to construct a list of GenPairs.
void Gen_MultiCombo::AddGenPair(G4String newValues)
{
    newValues = trim(newValues);
    vector<string> values = split(newValues, "|");
    Log::Assert((values.size() >= 3),
       "Call to Multicombo Set(Vtx/Pos)State must have at least 3 fields.");
    double w;
    try{
        w = to_double(values[2]);
    } catch(std::invalid_argument error) {
        Log::Die("Gen_MultiCombo: "+values[2]+
                " could not be parsed as a number");
    }
    Log::Assert((w >= 0),
        "Weight cannot be a negative number - skipping this origin set.");
    fGenPairs.push_back(
            new GenPair(fVertexGenType, values[0],
                fPosGenType, values[1], w));
    fWeightSum += w;
    fPairStateStr += newline + newValues;
    //TODO: need to check for success in creating the GenPair before
    //      adding it to the vector, weight sum and state string.
}

void Gen_MultiCombo::DisplayFormat()
{
    info
      << "Current state of this Gen_MultiCombo:" << newline
      << " \"" << GetPosState() << "\""             << newline;
    detail
      << "Format of call to Gen_MultiCombo::SetPosState"
      << " or Gen_MultiCombo::SetVtxState:"             << newline
      << "\"/generator/vtx/set vtxState|posState|weight\"" << newline
      << "where:"                                          << newline
      << "\t vtxState is the state to pass to VertexGen"   << newline
      << "\t posState is the state to pass to PosGen"      << newline
      << "\t and weight is the weight to give this pair"   << newline;
}
} //namespace RAT