////////////////////////////////////////////////////////////////////
/// \class RAT::DecayChain
///
/// \brief Sets decay chain isotope kinematic, etc information
///
/// \author Joe Formaggio and Jason Detwiler -- not on SNO+
/// \author Logan Sibley <lsibley@ualberta.ca> -- contact person
///
/// REVISION HISTORY:\n
///    04 Jan 2006 : WGS - Drag this code into the 21st century: use CLHEP,
///                  use strings, use STL. \n
///    02 Nov 2009 : Aleksandra Bialek - modifications \n
///    28 May 2010 : Logan Sibley - changes to accommodate screen option \n
///    12 Aug 2012 : Logan Sibley - positron emission/electron capture now
///                  handled correctly \n
///    23 Sep 2012 : Logan Sibley - gamma decay now taken into account and is
///                  used in place of the previous change (12 Aug 2010) to
///                  handle electron capture. Positron emission still okay. \n
///    03 Jan 2013 : Logan Sibley - Now adds parent info into DS \n
///    11 Nov 2014 : Matthew Strait - Don't set a variable right before it
///                  goes out of scope
///
///
/// \details
///
////////////////////////////////////////////////////////////////////
#ifndef __RAT_DecayChain__
#define __RAT_DecayChain__

#include <RAT/BetaFunction.hh>

#include <CLHEP/Vector/LorentzVector.h>

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <map>

namespace RAT{

static const double N_A = 6.022e+23;

class DecayChain
{
  public:
    // Structure to hold (particle ID, four-vector):
    typedef struct {
      int ID;
      CLHEP::HepLorentzVector vector;
    }
    ParticleInfo_t;

    DecayChain();
    DecayChain(const std::string Name);
    ~DecayChain();

    void Reset();
    void Show();

    void AddElement(const std::string Name, int iDecay = DecayBeta, int iLocation = 1, double tau = 1., double wt = 1.);
    void RemoveElement(int iBranch);

    void SetElementName(int iBranch, const std::string Name);
    void SetElementNumber(int iBranch);
    void SetElementBranch(int iBranch, int iLocation);
    void SetElementWt(int iBranch, double wt = 1.);
    void SetElementDecay(int iBranch, int iDecay);
    void SetLifetime(int iBranch, double tau);

    bool ReadInputFile(const std::string dName);

    const std::string GetElementName(int iBranch);
    int GetElementBranch(int iBranch);
    int GetElementDecay(int iBranch);
    double GetElementWt(int iBranch);
    double GetLifetime(int iBranch);

    void GenerateDecayElement(int iBranch);
    void GenerateDecayElement(const std::string iElement);

    void GenerateFullChain(const std::string iElement);
    void GenerateFullChain(int iBranch = 0);

    void GenerateDecayChain(const std::string iNameStart);
    void GenerateDecayChain(int iStart = 0);

    const std::string GetParentName(int iEvent);
    double GetParentMass(int iEvent);
    double GetParentCharge(int iEvent);
    const std::string GetDaughterName(int iEvent);
    int GetEventID(int iEvent);
    double GetEventEnergy(int iEvent);
    double GetEventTime(int iEvent);

    bool CheckParentMatch(int iEvent);

    void SetDecaySequence();
    void SetElementProb(int iBranch, double prob);
    void SetEventTime(int iEvent, double time);
    double GetElementProb(int iBranch);

    int  CheckChainLevel(int iElement);
    bool CheckInChain(int iElement, int jElement);

    void SetParticleInfo(int iTag, int jTag, int kTag);

    ParticleInfo_t AddDaughterInfo(int iTag);
    ParticleInfo_t GetParticleInfo(int iEvent);

    double GetRandomNumber(double rmin = 0., double rmax = 1.);
    inline void SetInMiddleChain(bool inMiddle=false) {fInMiddleChain = inMiddle;};
    inline void SetAlphaDecayStart(bool alphaStart=false) {fAlphaDecayStart = alphaStart;};
    inline void SetGammaDecayStart(bool gammaStart=false) {fGammaDecayStart = gammaStart;};
    inline void SetScreenCorr(bool screenCorr=false) {fSetScreenCorr = screenCorr;};
    inline void SetIsVerbose(bool iSet=false) {fIsVerbose = iSet;};
    inline void SetInputFileName(const std::string Name) {fInputFileName = Name;};
    inline void SetChainName(const std::string Name) {fChainName = Name;};
    inline void SetIsEquilibrium(bool iEq = true) {fIsEquilibrium = iEq;};
    inline void SelectParent(const std::string iElement = "ALL") {fSelectParent = iElement;};
    inline void SetNGenerated(int nP) {fNGenerated = nP;};
    inline void IncludeDaughter(bool iSelect = true){fIncludeDaughter = iSelect;};

    inline bool GetIsEquilibrium() {return fIsEquilibrium;};
    inline bool GetInMiddleChain() {return fInMiddleChain;};
    inline bool GetGammaDecayStart() {return fGammaDecayStart;};
    inline int GetScreenCorr() {return fSetScreenCorr;};
    inline int GetNumberOfElements() {return fNumberOfElements;};
    inline int GetElementNumber() {return fElementNumber;};
    inline int GetNGenerated() {return fNGenerated;};

    inline const std::string GetInputFileName() {return fInputFileName;};
    inline const std::string GetChainName() {return fChainName;};

  private:
    inline void SetParentName(int iParent, const std::string name) {fParentName[iParent] = name;};
    inline void SetParentMass(int iParent, double A) {fParentMass[iParent] = A;};
    inline void SetParentCharge(int iParent, double Z) {fParentCharge[iParent] = Z;};
    inline void SetDaughterName(int iDaughter, const std::string name) {fDaughterName[iDaughter] = name;};

    bool fIsVerbose;
    bool fIsEquilibrium;
    bool fIsChainElement;
    bool fInMiddleChain;  // to start chain at the isotope defined in *.mac
    bool fAlphaDecayStart;// to start chain at the isotope with alpha decay
    bool fGammaDecayStart;// to start chain at the isotope with gamma decay
    bool fSetScreenCorr;  // to set screening correction

    std::string fChainName;

    int fNumberOfElements;
    int fElementNumber;
    std::map<size_t, std::string> fElementName;
    std::map<size_t, int> fElementBranch;
    std::map<size_t, int> fElementDecay;
    std::map<size_t, double> fElementWeight;
    std::map<size_t, double> fElementLife;
    std::map<size_t, double> fElementProb;
    double fTstart;

    std::map<size_t, BetaFunction*> fElement;

    bool fIncludeDaughter;

    std::string fInputFileName;
    std::string fSelectParent;

    int fNGenerated;
    std::map<size_t, std::string> fParentName;
    std::map<size_t, double> fParentMass;
    std::map<size_t, double> fParentCharge;
    std::map<size_t, std::string> fDaughterName;
    std::map<size_t, int> fParticleID;
    std::map<size_t, double> fParticleEnergy;
    std::map<size_t, double> fParticleTime;

    std::map<size_t, ParticleInfo_t> fParticleInfo;

};

}
#endif