// Sc48SourceGen.cc
// Contact person: Valentina Lozza <valentina.lozza@tu-dresden.de>
// See Sc48SourceGen.hh for more details
//———————————————————————//

#include <CLHEP/Units/SystemOfUnits.h>
using namespace CLHEP;

#include <RAT/Sc48SourceGen.hh>
#include <RAT/Detector.hh>
#include <RAT/Log.hh>
#include <RAT/string_utilities.hpp>


#include <G4UIcommand.hh>

#include <TMath.h>

#include <string>
#include <algorithm>

namespace RAT {

Sc48SourceGen::Sc48SourceGen() : GLG4Gen("Sc48Source") ,fStateStr(""), fActiveVolumeName("Sc48Source_sample_logical")
{
    fNextGen = NULL;

    // half lives of the nuclides.
    // only the nuclides in this map are accepted in the state string.

    fHalfLives["46Sc"] = 83.79  *  24 * 3600 * s;
    fHalfLives["47Sc"] = 3.3492 *  24 * 3600 * s;
    fHalfLives["48Sc"] =        43.67 * 3600 * s;

    if (Detector::FindPhysicalVolume(fActiveVolumeName) == NULL) {
        Log::Die("Sc48SourceGen: Did not find physical volume '" + fActiveVolumeName + "'. Please load the Sc48Source geometry.");
    }
}

Sc48SourceGen::~Sc48SourceGen()
{
    std::map<G4String, DecayChain_Gen*>::iterator itGenerators;
    for(itGenerators = fGenerators.begin(); itGenerators != fGenerators.end(); itGenerators++){
        delete itGenerators->second;
        itGenerators->second = NULL;
    }
}

void Sc48SourceGen::GenerateEvent(G4Event *event)
{
    if (fNextGen != NULL) {
        fNextGen->GenerateEvent(event);
        fNextGen = NULL;
    }
    else {
        Log::Die("Sc48SourceGen: no nuclide determined to be used next in GenerateEvent()!");
    }
}

void Sc48SourceGen::SubtractTime(double time)
{
    // Subtract time for the next generator ...
    fNextGen->SubtractTime(time);
    // ... which also affects the resulting fNextTime.
    fNextTime -= time;
}

double Sc48SourceGen::GetActivityAfterDecayTime(const G4String &nuclide, const double initialActivity, const double decayTime) const
{
    detail << "Sc48SourceGen: Activity of " << nuclide << " initially " << initialActivity << "Bq.";

    if (decayTime == 0) {
        return initialActivity;
    }
    else {
        double newActivity = initialActivity * TMath::Exp(- TMath::Log(2) * decayTime/fHalfLives.at(nuclide));
        detail << " -> Set to " << to_string(newActivity) << "Bq after " << decayTime / (3600*s) << " hours." << "\n";
        return newActivity;
    }
}

void Sc48SourceGen::SetState(G4String state)
{

    info << "Sc48SourceGen: Initializing with " << state << "\n";
    fStateStr = state;

    G4String explanation =
        "Arguments to the scandium generator should look like this: \n "
        " '/generator/add sc48source 48Sc:100 47Sc:10 46Sc:1 DecayTime:5day', \n"
        "where the numbers for the nuclides are the activity (in Bq) after activation and the DecayTime is the time passed after the activation. \n"
        "Single nuclides may be omitted, but at least one must be given.";

    std::vector<std::string> paramStrings, splitParameter;

    std::map<G4String, G4String> params;
    std::map<G4String, G4String>::iterator itParams;

    if (state != "") {
        paramStrings = split(state, " ");
        detail << "Sc48SourceGen: Found " << paramStrings.size() << " parameters in state." << "\n";

        for (unsigned int i = 0; i < paramStrings.size(); i++) {
            splitParameter = split(paramStrings.at(i), ":");
            if (params.find(splitParameter.at(0)) != params.end()) {
                Log::Die("Sc48SourceGen: more than one value given for parameter " + splitParameter.at(0) + " in state string.");
            }
            params[splitParameter.at(0)] = splitParameter.at(1);
        }

        // Check first, whether a DecayTime is given:
        double decayTime = 0;
        itParams = params.find("DecayTime");
        if (itParams != params.end()) {
            decayTime = G4UIcommand::ConvertToDimensionedDouble(itParams->second);
            if (decayTime == 0 && to_double(itParams->second) != 0) {
                Log::Die("Sc48SourceGen: No unit given for DecayTime parameter, or unit not recognized.");
            }
            info << "Sc48SourceGen: Activities will be calculated for a decay time of " << decayTime/(3600*s) << " hours." << "\n";
            params.erase(itParams);
        }
        else {
            detail << "Sc48SourceGen: Time from preparation to measurement of the source set to 0. You can give an additional argument like: DecayTime:5day." << "\n";
        }

        // The remaining parameters in the map should be of type params[Nuclide] == Activity in Bq.
        for (itParams = params.begin(); itParams != params.end(); itParams++) {
            if (fHalfLives.find(itParams->first) == fHalfLives.end()) {
                Log::Die("Sc48SourceGen: Unknown nuclide '" + itParams->first + "' given as generator state. \n " + explanation);
            }

            fGenerators[itParams->first] = new DecayChain_Gen();
            fGenerators[itParams->first]->SetState(itParams->first + ":fill:poisson");
            fGenerators[itParams->first]->SetPosState(fActiveVolumeName);
            fGenerators[itParams->first]->SetTimeState(to_string(GetActivityAfterDecayTime(itParams->first, to_double(itParams->second), decayTime)));
        }

        Log::Assert(fGenerators.size() != 0,
                    "Sc48SourceGen: No generators created. Were any nuclides given in the generator state? \n" + explanation);
    }
    else {
        Log::Die("Sc48SourceGen: No state found. \n" + explanation);
    }
}

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

void Sc48SourceGen::ResetTime(double offset)
{
    // Clear next generator
    fNextGen = NULL;

    // Reset time for each generator and determine next generator
    std::map<G4String, DecayChain_Gen*>::iterator itGenerators;
    for(itGenerators = fGenerators.begin(); itGenerators != fGenerators.end(); itGenerators++){
        itGenerators->second->ResetTime(offset);
        if (fNextGen == NULL || itGenerators->second->NextTime() < fNextTime) {
            fNextGen = itGenerators->second;
            fNextTime = itGenerators->second->NextTime();
        }
    }
}

void Sc48SourceGen::SetPosState(G4String /*state*/)
{
    warn << "Warning: /generator/pos/set has no effect for the sc48source generator.\n";
}

void Sc48SourceGen::SetTimeState(G4String /*state*/)
{
    warn << "Warning: /generator/rate/set has no effect for the sc48source generator.\n";
}

} // namespace RAT