//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// ABLAXX statistical de-excitation model
// Pekka Kaitaniemi, HIP (translation)
// Christelle Schmidt, IPNL (fission code)
// Davide Mancusi, CEA (contact person INCL/ABLA)
// Aatos Heikkinen, HIP (project coordination)
//
#define ABLAXX_IN_GEANT4_MODE 1

#include "globals.hh"

// Data structures needed by ABLA evaporation code.

#ifndef G4AblaDataDefs_hh
#define G4AblaDataDefs_hh 1

#ifdef ABLAXX_IN_GEANT4_MODE
#include "globals.hh"
#else
#include "G4INCLGeant4Compat.hh"
#endif

#include <cmath>

// ABLA

class G4Nevent {
public:
  G4Nevent() {};
  ~G4Nevent() {};
  
  G4int ii;
};

// ABLA
#define PACESIZEROWS 500
#define PACESIZECOLS 500
/**
 * Masses.
 */

class G4Pace {

public:
  G4Pace() {};

  ~G4Pace() {};
  
  G4double dm[PACESIZEROWS][PACESIZECOLS];
};

#define EC2SUBROWS 154
#define EC2SUBCOLS 99
  /**
   *
   */

class G4Ec2sub {
public:
  G4Ec2sub() {};

  ~G4Ec2sub() {};

  G4double ecnz[EC2SUBROWS][EC2SUBCOLS]; 

  /**
   * Dump the contents of the ecnz data table.
   */
  void dump() {
    for(G4int i = 0; i < EC2SUBROWS; i++) {
      for(G4int j = 0; j < EC2SUBCOLS; j++) {
	//G4cout << ecnz[i][j] << " ";
      }
      //      G4cout << G4endl;
    }
  }
};

class G4Ald {
public:
  /**
   * 
   */
  G4Ald()
    :av(0.0), as(0.0), ak(0.0), optafan(0.0)
  {};
  ~G4Ald() {};
  
  G4double av,as,ak,optafan;
};

#define ECLDROWS 154
#define ECLDCOLS 99
/**
 * Shell corrections and deformations.
 */

class G4Ecld {

public:
  G4Ecld() {};
  ~G4Ecld() {};

  /**
   * Ground state shell correction frldm for a spherical ground state.
   */
  G4double ecgnz[ECLDROWS][ECLDCOLS];

  /**
   * Shell correction for the saddle point (now: == 0).
   */
  G4double ecfnz[ECLDROWS][ECLDCOLS];

  /**
   * Difference between deformed ground state and ldm value.
   */
  G4double vgsld[ECLDROWS][ECLDCOLS]; 

  /**
   * Alpha ground state deformation (this is not beta2!)       
   * beta2 = std::sqrt(5/(4pi)) * alpha 
   */
  G4double alpha[ECLDROWS][ECLDCOLS];
};

class G4Fiss {
  /**
   * Options and parameters for fission channel.
   */

public:
  G4Fiss()
    :akap(0.0), bet(0.0), homega(0.0), koeff(0.0), ifis(0.0),
     optshp(0), optxfis(0), optles(0), optcol(0)
  {};
  ~G4Fiss() {};
  
  G4double akap,bet,homega,koeff,ifis;
  G4int optshp, optxfis,optles,optcol;
};

#define FBROWS 101
#define FBCOLS 161
/**
 * Fission barriers.
 */
 
class G4Fb {
  
public:
  G4Fb() {};
  ~G4Fb() {;}
  
  //  G4double efa[FBROWS][FBCOLS];
  G4double efa[FBCOLS][FBROWS];
};

/**
 * Options
 */

class G4Opt {

public:
  G4Opt()
    :optemd(0), optcha(0), eefac(0.0)
  {};
  ~G4Opt() {};
  
  G4int optemd,optcha;
  G4double eefac;                                  
};

#define EENUCSIZE 2002
#define XHESIZE 50
class G4Eenuc {
public:
  G4Eenuc() {
    for(G4int i = 0; i < EENUCSIZE; ++i) {
      she[i] = 0.0;
    }
    for(G4int i = 0; i < XHESIZE; ++i) {
      for(G4int j = 0; j < EENUCSIZE; ++j) {
	xhe[i][j] = 0.0;
      }
    }
  };
  ~G4Eenuc() {};
  
  G4double she[EENUCSIZE],xhe[XHESIZE][EENUCSIZE];                                            
};

//#define VOLANTSIZE 200
#define VOLANTSIZE 301
/**
 * Evaporation and fission output data.
 */

class G4Volant {
  
public:
  G4Volant()
  {
    clear();
  }

  ~G4Volant() {};

  void clear()
  {
    for(G4int i = 0; i < VOLANTSIZE; i++) {
      copied[i] = false;
      acv[i] = 0;
      zpcv[i] = 0;
      pcv[i] = 0;
      xcv[i] = 0;
      ycv[i] = 0;
      zcv[i] = 0;
      iv = 0;
    }
  }

  G4double getTotalMass()
  {
    G4double total = 0.0;
    for(G4int i = 0; i <= iv; i++) {
      total += acv[i];
    }
    return total;
  }

  void dump()
  {
    G4double totA = 0.0, totZ = 0.0, totP = 0.0;
    //    G4cout <<"i \t ACV \t ZPCV \t PCV" << G4endl; 
    for(G4int i = 0; i <= iv; i++) {
      if(i == 0 && acv[i] != 0) {
	//	G4cout <<"G4Volant: Particle stored at index " << i << G4endl;
      }
      totA += acv[i];
      totZ += zpcv[i];
      totP += pcv[i];
      //      G4cout << "volant" << i << "\t" << acv[i] << " \t " << zpcv[i] << " \t " << pcv[i] << G4endl;
    }
    //    G4cout <<"Particle count index (iv) = " << iv << G4endl;
    //    G4cout <<"ABLA Total: A = " << totA << " Z = " << totZ <<  " momentum = " << totP << G4endl;
  }

  G4double acv[VOLANTSIZE],zpcv[VOLANTSIZE],pcv[VOLANTSIZE],xcv[VOLANTSIZE];
  G4double ycv[VOLANTSIZE],zcv[VOLANTSIZE];
  G4bool copied[VOLANTSIZE];
  G4int iv; 
};

#define VARNTPSIZE 301
class G4VarNtp {
public:
  G4VarNtp() {
    clear();
  };

  ~G4VarNtp() {};

  /**
   * Clear and initialize all variables and arrays.
   */
  void clear() {
    particleIndex = 0;
    projType = 0;
    projEnergy = 0.0;
    targetA = 0;
    targetZ = 0;
    masp = 0.0; mzsp = 0.0; exsp = 0.0; mrem = 0.0;
    // To be deleted?
    spectatorA = 0;
    spectatorZ = 0;
    spectatorEx = 0.0;
    spectatorM = 0.0;
    spectatorT = 0.0;
    spectatorP1 = 0.0;
    spectatorP2 = 0.0;
    spectatorP3 = 0.0;
    massini = 0;
    mzini = 0;
    exini = 0;
    pcorem = 0;
    mcorem = 0;
    pxrem = 0;
    pyrem = 0;
    pzrem = 0;
    erecrem = 0;
    mulncasc = 0;
    mulnevap = 0;
    mulntot = 0;
    bimpact = 0.0;
    jremn = 0;
    kfis = 0;
    estfis = 0;
    izfis = 0;
    iafis = 0;
    ntrack = 0;
    needsFermiBreakup = false;
    for(G4int i = 0; i < VARNTPSIZE; i++) {
      itypcasc[i] = 0;
      avv[i] = 0;
      zvv[i] = 0;
      enerj[i] = 0.0;
      plab[i] = 0.0;
      tetlab[i] = 0.0;
      philab[i] = 0.0;
      full[i] = false;
    }
  }

  /**
   * Add a particle to the INCL/ABLA final output.
   */
  void addParticle(G4double A, G4double Z, G4double E, G4double P, G4double theta, G4double phi) {
    if(full[particleIndex]) {
      //      G4cout <<"A = " << Z << " Z = " << Z << G4endl;
    } else {
      avv[particleIndex] = (int) A;
      zvv[particleIndex] = (int) Z;
      enerj[particleIndex] = E;
      plab[particleIndex] = P;
      tetlab[particleIndex] = theta;
      philab[particleIndex] = phi;
      full[particleIndex] = true;
      ntrack = particleIndex + 1;
      particleIndex++;
    }
  }

  /**
   * Baryon number conservation check.
   */
  G4int getTotalBaryonNumber() {
    G4int baryonNumber = 0;
    for(G4int i = 0; i < ntrack; i++) {
      if(avv[i] > 0) {
	baryonNumber += avv[i];
      }
    }
    return baryonNumber;
  }

  /**
   * Return total energy.
   */
  G4double getTotalEnergy() {
    G4double energy = 0.0;
    for(G4int i = 0; i < ntrack; i++) {
      energy += std::sqrt(std::pow(plab[i], 2) + std::pow(getMass(i), 2)); // E^2 = p^2 + m^2
    }

    return energy;
  }

  /**
   * Return total three momentum.
   */
  G4double getTotalThreeMomentum() {
    G4double momentum = 0;
    for(G4int i = 0; i < ntrack; i++) {
      momentum += plab[i];
    }
    return momentum;
  }

  G4double getMomentumSum() {
    G4double momentum = 0;
    for(G4int i = 0; i < ntrack; i++) {
      momentum += plab[i];
    }
    return momentum;
  }

  G4double getMass(G4int particle) {
    const G4double protonMass = 938.272;
    const G4double neutronMass = 939.565;
    const G4double pionMass = 139.57;

    G4double mass = 0.0;
    if(avv[particle] ==  1 && zvv[particle] ==  1) mass = protonMass;
    if(avv[particle] ==  1 && zvv[particle] ==  0) mass = neutronMass;
    if(avv[particle] == -1)                        mass = pionMass;
    if(avv[particle] > 1)
      mass = avv[particle] * protonMass + zvv[particle] * neutronMass;
    return mass;
  }

  /**
   * Dump debugging output.
   */
  void dump() {
    G4int nProton = 0, nNeutron = 0;
    G4int nPiPlus = 0, nPiZero = 0, nPiMinus = 0;
    G4int nH2 = 0, nHe3 = 0, nAlpha = 0;
    G4int nFragments = 0;
    G4int nParticles = 0;
    for(G4int i = 0; i < ntrack; i++) {
      nParticles++;
      if(avv[i] ==  1 && zvv[i] ==  1) nProton++;  // Count multiplicities
      if(avv[i] ==  1 && zvv[i] ==  0) nNeutron++;
      if(avv[i] == -1 && zvv[i] ==  1) nPiPlus++;
      if(avv[i] == -1 && zvv[i] ==  0) nPiZero++;
      if(avv[i] == -1 && zvv[i] == -1) nPiMinus++;
      if(avv[i] ==  2 && zvv[i] ==  1) nH2++;
      if(avv[i] ==  3 && zvv[i] ==  2) nHe3++;
      if(avv[i] ==  4 && zvv[i] ==  2) nAlpha++;
      if(                zvv[i] >   2) nFragments++;
    }
  }

  /**
   * Projectile type.
   */
  G4int projType;

  /**
   * Projectile energy.
   */
  G4double projEnergy;

  /**
   * Target mass number.
   */
  G4int targetA;

  /**
   * Target charge number.
   */
  G4int targetZ;

  /**
   * Projectile spectator A, Z, Eex;
   */
  G4double masp, mzsp, exsp, mrem;

  /**
   * Spectator nucleus mass number for light ion projectile support.
   */
  G4int spectatorA;

  /**
   * Spectator nucleus charge number for light ion projectile support.
   */
  G4int spectatorZ;

  /**
   * Spectator nucleus excitation energy for light ion projectile support.
   */
  G4double spectatorEx;

  /**
   * Spectator nucleus mass.
   */
  G4double spectatorM;

  /**
   * Spectator nucleus kinetic energy.
   */
  G4double spectatorT;

  /**
   * Spectator nucleus momentum x-component.
   */
  G4double spectatorP1;

  /**
   * Spectator nucleus momentum y-component.
   */
  G4double spectatorP2;

  /**
   * Spectator nucleus momentum z-component.
   */
  G4double spectatorP3;

  /**
   * A of the remnant.
   */
  G4double massini;

  /**
   * Z of the remnant.
   */
  G4double mzini;

  /**
   * Excitation energy.
   */
  G4double exini;

  G4double pcorem, mcorem, pxrem, pyrem, pzrem, erecrem;

  /**
   * Cascade n multip.
   */
  G4int mulncasc;

  /**
   * Evaporation n multip.
   */
  G4int mulnevap;

  /**
   * Total n multip.
   */
  G4int mulntot;

  /**
   * Impact parameter.
   */
  G4double bimpact;

  /**
   * Remnant Intrinsic Spin.
   */
  G4int jremn;

  /**
   * Fission 1/0=Y/N.
   */
  G4int kfis;

  /**
   * Excit energy at fis.
   */
  G4double estfis;

  /**
   * Z of fiss nucleus.
   */
  G4int izfis;

  /**
   * A of fiss nucleus.
   */
  G4int iafis;

  /**
   * Number of particles.
   */
  G4int ntrack;

  /**
   * The state of the index:
   * true = reserved
   * false = free
   */
  G4bool full[VARNTPSIZE];

  /**
   * Does this nucleus require Fermi break-up treatment? Only
   * applicable when used together with Geant4.
   * true = do fermi break-up (and skip ABLA part)
   * false = use ABLA
   */
  G4bool needsFermiBreakup;

  /**
   * emitted in cascade (0) or evaporation (1).
   */
  G4int itypcasc[VARNTPSIZE];

  
  /**
   * A (-1 for pions).
   */
  G4int avv[VARNTPSIZE];

  /**
   * Z
   */
  G4int zvv[VARNTPSIZE];

  /**
   * Kinetic energy.
   */
  G4double enerj[VARNTPSIZE];

  /**
   * Momentum.
   */
  G4double plab[VARNTPSIZE];

  /**
   * Theta angle.
   */
  G4double tetlab[VARNTPSIZE];

  /**
   * Phi angle.
   */
  G4double philab[VARNTPSIZE];

private:
  G4int particleIndex;
};

#endif