// from ROOT
#include "TMath.h"
#include "TString.h"

// from hydra
#include "hkalmixture.h"
#include "hkaldef.h"
#include "htool.h"
using namespace Kalman;

#include <iostream>
#include <iomanip>
using namespace std;

//_HADES_CLASS_DESCRIPTION
///////////////////////////////////////////////////////////////////////////////
//
// Class: HKalMixture
//
//
//
///////////////////////////////////////////////////////////////////////////////


//  -----------------------------------
//  Ctors and Dtor
//  -----------------------------------

HKalMixture::HKalMixture() : TMixture() {
    fDensityMixt    = NULL;
    fExciteEnerMixt = NULL;
    fRadLengthMixt  = NULL;
}

HKalMixture::HKalMixture(const char *name, const char *title, Int_t nmixt)
: TMixture(name, title, nmixt) {
    // Make a mixture that consists of nmixt different materials.
    // To add materials use the defineElement() method. After all materials
    // have been added, the atomic and nuclear properties of the mixture have
    // to be calculated with the calcMean() function.
    // name:   material name
    // title:  title
    // nmixt:  number of materials in the mixture

    if (nmixt == 0) {
        fDensityMixt    = NULL;
        fExciteEnerMixt = NULL;
        fRadLengthMixt  = NULL;
        return;
    }
    Int_t nm = TMath::Abs(nmixt);
    fDensityMixt    = new Float_t[nm];
    fExciteEnerMixt = new Float_t[nm];
    fRadLengthMixt  = new Float_t[nm];
}

HKalMixture::~HKalMixture() {
    delete [] fDensityMixt;
    delete [] fExciteEnerMixt;
    delete [] fRadLengthMixt;
    fDensityMixt    = NULL;
    fExciteEnerMixt = NULL;
    fRadLengthMixt  = NULL;
}

//  -----------------------------------
//  Implementation of public methods
//  -----------------------------------

void HKalMixture::calcMassFracFromVolFrac(Float_t w[], Int_t n, const Float_t rho[], const Float_t vol[]) {
    // This function will calculate the mass fractions of a mixture given the components' densities and
    // volume fractions. Only works for ideal mixtures.
    // w: the mass weights (return value)
    // n: the number of components in the mixture. All arrays must have at least n elements.
    // rho: the densities of the components
    // vol: the volume fractions

    Float_t rhoMixt = 0;
    // Density of mixture.
    for(Int_t i = 0; i < n; i++) {
        rhoMixt += rho[i] * vol[i];
    }
    // Calculate mass fraction in mixture from volume fractions.
    // Formula from Geant manual
    for(Int_t i = 0; i < n; i++) {
        w[i] = vol[i] * rho[i] / rhoMixt;
    }
}


void HKalMixture::calcMixture() {
    // Calculates properties of the mixture:
    // atomic mass A, atomic number Z, density, radiation length and mean excitation energy.

    fA          = 0;
    fExciteEner = 0.;
    fDensity    = 0;
    fRadLength  = 0;
    fZ          = 0;
    Float_t ZA  = 0.;
    Float_t lnExEnerWeight = 0.;

    // Effective A, Z and density are done as in Geant 3.16, chapter CONS110.
    // The mean excitation energy is calculated using Bragg's rule:
    // S.M. Seltzer and M.J. Berger, Int. J. of Applied Rad. 33, 1189 (1982).
    //
    // Calculating density of a compound:
    // rho = sum(m_i)   / V                   with V = sum(V_i) = sum(m_i/rho_i)
    //     = sum(m_i)   / sum(m_i/rho_i)      m_i = p_i * M
    //     = sum(p_i*M) / sum(p_i*M/rho_i)
    //     = sum(p_i)*M / (sum(p_i/rho_i)*M)
    //     = sum(p_i)   / sum(p_i/rho_i)      sum(pi) = 1
    //     = 1 / (sum(p_i/rho_i)
    // p_i = m_i/M  mass fraction of component i

    for(Int_t i = 0; i < GetNmixt(); i++) {
        if(fWmixt[i] != 0.) {
            fA             += fWmixt[i] * fAmixt[i];
            fZ             += fWmixt[i] * fZmixt[i];
            // Add weighted radiation lengths of components in 1/(g/mm^2)
            fRadLength     += fWmixt[i] / (fDensityMixt[i] * fRadLengthMixt[i]);
            fDensity       += fWmixt[i] / fDensityMixt[i];
            ZA             += fWmixt[i] * fZmixt[i] / fAmixt[i];
            lnExEnerWeight += TMath::Log(fExciteEnerMixt[i]) * fWmixt[i] * fZmixt[i] / fAmixt[i];
        }
    }
    if(fDensity != 0.) {
        fDensity = 1. / fDensity;
    }
    if(fRadLength != 0.) {
        //fRadLength *= fDensity; // total radiation length in mm
        fRadLength = 1. / (fDensity * fRadLength);
    }
    fExciteEner = TMath::Exp(lnExEnerWeight/ZA);
}

void HKalMixture::defineElement(Int_t n, const HKalMixture &mat, Float_t w) {
    defineElement(n, mat.GetA(), mat.GetZ(), mat.GetDensity(), mat.GetExciteEner(), mat.GetRadLength(), w);
}

void HKalMixture::defineElement(Int_t n, const Float_t mat[5], Float_t w) {
    defineElement(n, mat[kMatIdxA], mat[kMatIdxZ], mat[kMatIdxDensity], mat[kMatIdxExEner], mat[kMatIdxRadLength], w);
}

void HKalMixture::defineElement(Int_t n, Float_t a, Float_t z, Float_t dense, Float_t exener,
        Float_t radl, Float_t w) {
    if (n < 0 || n >= TMath::Abs(fNmixt)) {
        Error("defineElement()", Form("Invalid element index %i.", n));
        return;
    }

    DefineElement(n, a, z, w); // parent's method
    fDensityMixt   [n] = dense;
    fExciteEnerMixt[n] = exener;
    fRadLengthMixt [n] = radl;

    //cout<<"define Element "<<n<<" of "<<fNmixt<<endl;
    //this->print("mixt, elements");
}

void HKalMixture::print(const char *opt) const {
    cout<<"**** Print content of material "<<GetName()<<" ****"<<endl;
    TString stropt(opt);
    if(stropt.Contains("mix", TString::kIgnoreCase) || stropt.IsNull()) {
        cout<<"Mixture properties:"<<endl;
        cout<<"A                      = "<<GetA()<<endl;
        cout<<"Z                      = "<<GetZ()<<endl;
        cout<<"density                = "<<GetDensity()<<endl;
        cout<<"mean excitation energy = "<<GetExciteEner()<<endl;
        cout<<"radiation length       = "<<GetRadLength()<<endl;
    }
    if(stropt.Contains("elements", TString::kIgnoreCase)) {
        for(Int_t i = 0; i < GetNmixt(); i++) {
            cout<<"Properties of material "<<i<<" in mixture:"<<endl;
            cout<<"A                      = "<<setprecision(3)<<fAmixt[i]<<endl;
            cout<<"Z                      = "<<setprecision(3)<<fZmixt[i]<<endl;
            cout<<"density                = "<<setprecision(8)<<fDensityMixt[i]<<endl;
            cout<<"radiation length       = "<<setprecision(8)<<fRadLengthMixt[i]<<endl;
            cout<<"mean excitation energy = "<<setprecision(8)<<fExciteEnerMixt[i]<<endl;
            cout<<"material weight        = "<<setprecision(3)<<fWmixt[i]<<endl;
        }
    }
    cout << resetiosflags(ios::showpoint) << endl;

    cout<<"**** End print of mixture ****"<<endl;
}