#include "G4BGGNucleonElasticXS_MOD.hh"
#include "G4SystemOfUnits.hh"
#include "G4ComponentGGHadronNucleusXsc.hh"
#include "G4NucleonNuclearCrossSection.hh"
#include "G4HadronNucleonXsc.hh"
#include "G4NuclearRadii.hh"
#include "G4Proton.hh"
#include "G4Neutron.hh"
#include "G4NistManager.hh"
#include "G4NuclearRadii.hh"

#include "G4CrossSectionDataSetRegistry.hh"

G4double G4BGGNucleonElasticXS_MOD::theGlauberFacP[93] = {0.0};
G4double G4BGGNucleonElasticXS_MOD::theCoulombFacP[93] = {0.0};
G4double G4BGGNucleonElasticXS_MOD::theGlauberFacN[93] = {0.0};
G4double G4BGGNucleonElasticXS_MOD::theCoulombFacN[93] = {0.0};
G4int    G4BGGNucleonElasticXS_MOD::theA[93] = {0};

#ifdef G4MULTITHREADED
G4Mutex G4BGGNucleonElasticXS_MOD::nucleonElasticXSMutex = G4MUTEX_INITIALIZER;
#endif

G4BGGNucleonElasticXS_MOD::G4BGGNucleonElasticXS_MOD(const G4ParticleDefinition* p)
 : G4VCrossSectionDataSet("BarashenkovGlauberGribov") 
{
  verboseLevel = 0;
  fGlauberEnergy = 91.*GeV;
  fLowEnergy = 14.0*MeV;
  fNucleon = nullptr;
  fGlauber = nullptr;
  fHadron  = nullptr;

  theProton= G4Proton::Proton();
  isProton = (theProton == p);
  isMaster = false;
  SetForAllAtomsAndEnergies(true);

  maxEn = 112900000.0;
  interp.setData(energy,CS);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4BGGNucleonElasticXS_MOD::~G4BGGNucleonElasticXS_MOD()
{
  delete fHadron;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4bool 
G4BGGNucleonElasticXS_MOD::IsElementApplicable(const G4DynamicParticle*, G4int Z,
                                           const G4Material* mat)
{

  return true;

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4bool G4BGGNucleonElasticXS_MOD::IsIsoApplicable(const G4DynamicParticle*, 
                                              G4int Z, G4int,  
                                              const G4Element*,
                                              const G4Material*)
{

  return (1 == Z);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4double
G4BGGNucleonElasticXS_MOD::GetElementCrossSection(const G4DynamicParticle* dp,
                                              G4int ZZ, const G4Material* mat)
{

  G4double cross = 0.0;

  G4String p_name = dp->GetDefinition()->GetParticleName();
  G4String mat_name = mat->GetName();
  G4double ekin = dp->GetKineticEnergy()/eV;

  if (p_name == "proton" and mat_name == "G4_Be" and ZZ == 4 and ekin < maxEn)
  {
    cross = interp(ekin)*barn;
  }
  else
  {
    G4int Z = std::min(ZZ, 92);
    if(1 == Z) {
      cross = 1.0115*GetIsoCrossSection(dp,1,1);
    } else {
      if(ekin <= fLowEnergy) {
        cross = (isProton) ? theCoulombFacP[Z] : theCoulombFacN[Z];
        cross *= CoulombFactor(ekin, Z);
      } else if(ekin > fGlauberEnergy) {
        cross = (isProton) ? theGlauberFacP[Z] : theGlauberFacN[Z];
        cross *= fGlauber->GetElasticGlauberGribov(dp, Z, theA[Z]);
      } else {
        cross = fNucleon->GetElasticCrossSection(dp, Z);
      }
    }
  }

  return cross;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4double
G4BGGNucleonElasticXS_MOD::GetIsoCrossSection(const G4DynamicParticle* dp, 
                                          G4int Z, G4int A, 
                                          const G4Isotope*,
                                          const G4Element*,
                                          const G4Material*)
{
  // this method should be called only for Z = 1
  fHadron->HadronNucleonXscNS(dp->GetDefinition(), theProton, 
                              dp->GetKineticEnergy());
  G4double cross = A*fHadron->GetElasticHadronNucleonXsc();

  return cross;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void G4BGGNucleonElasticXS_MOD::BuildPhysicsTable(const G4ParticleDefinition& p)
{
  if(nullptr != fNucleon) { return; }
  if(&p == theProton || &p == G4Neutron::Neutron()) {
    isProton = (theProton == &p);

  } else {
    G4ExceptionDescription ed;
    ed << "This BGG cross section is applicable only to nucleons and not to " 
       << p.GetParticleName() << G4endl; 
    G4Exception("G4BGGNucleonElasticXS_MOD::BuildPhysicsTable", "had001", 
                FatalException, ed);
    return;
  }

  fNucleon = new G4NucleonNuclearCrossSection();
  fGlauber = new G4ComponentGGHadronNucleusXsc();
  fHadron  = new G4HadronNucleonXsc();

  fNucleon->BuildPhysicsTable(p);

  if(0 == theA[0]) { 
#ifdef G4MULTITHREADED
    G4MUTEXLOCK(&nucleonElasticXSMutex);
    if(0 == theA[0]) { 
#endif
      isMaster = true;
#ifdef G4MULTITHREADED
    }
    G4MUTEXUNLOCK(&nucleonElasticXSMutex);
#endif
  } else {
    return;
  }

  if(isMaster && 0 == theA[0]) {

    theA[0] = theA[1] = 1;
    G4ThreeVector mom(0.0,0.0,1.0);
    G4DynamicParticle dp(theProton, mom, fGlauberEnergy);

    G4NistManager* nist = G4NistManager::Instance();
    G4double csup, csdn;

    if(verboseLevel > 0) {
      G4cout << "### G4BGGNucleonElasticXS_MOD::Initialise for "
             << p.GetParticleName() << G4endl;
    }

    for(G4int iz=2; iz<93; ++iz) {
      G4int A = G4lrint(nist->GetAtomicMassAmu(iz));
      theA[iz] = A;

      csup = fGlauber->GetElasticGlauberGribov(&dp, iz, A);
      csdn = fNucleon->GetElasticCrossSection(&dp, iz);
      theGlauberFacP[iz] = csdn/csup;
    }

    dp.SetDefinition(G4Neutron::Neutron());
    for(G4int iz=2; iz<93; ++iz) {
      csup = fGlauber->GetElasticGlauberGribov(&dp, iz, theA[iz]);
      csdn = fNucleon->GetElasticCrossSection(&dp, iz);
      theGlauberFacN[iz] = csdn/csup;

      if(verboseLevel > 0) { 
        G4cout << "Z= " << iz <<  "  A= " << theA[iz] 
               << " GFactorP= " << theGlauberFacP[iz]
               << " GFactorN= " << theGlauberFacN[iz] << G4endl;
      } 
    }

    theCoulombFacP[0] = theCoulombFacP[1] = 
    theCoulombFacN[0] = theCoulombFacN[1] = 1.0; 
    dp.SetDefinition(theProton);
    dp.SetKineticEnergy(fLowEnergy);
    for(G4int iz=2; iz<93; ++iz) {
      theCoulombFacP[iz] = fNucleon->GetElasticCrossSection(&dp, iz)
        /CoulombFactor(fLowEnergy, iz);
    }
    dp.SetDefinition(G4Neutron::Neutron());
    for(G4int iz=2; iz<93; ++iz) {
      theCoulombFacN[iz] = fNucleon->GetElasticCrossSection(&dp, iz)
        /CoulombFactor(fLowEnergy, iz);

      if(verboseLevel > 0) {
        G4cout << "Z= " << iz <<  "  A= " << theA[iz]
               << " CFactorP= " << theCoulombFacP[iz] 
               << " CFactorN= " << theCoulombFacN[iz] << G4endl; 
      }
    }
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4double G4BGGNucleonElasticXS_MOD::CoulombFactor(G4double kinEnergy, G4int Z)
{
  G4double res= 1.0;
  if(isProton) {
    res = G4NuclearRadii::CoulombFactor(Z, theA[Z], theProton, kinEnergy);
  }
  return res;  
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void G4BGGNucleonElasticXS_MOD::CrossSectionDescription(std::ostream& outFile) const
{
  outFile << "The Barashenkov-Glauber-Gribov cross section handles elastic\n"
          << "scattering of protons and neutrons from nuclei using the\n"
          << "Barashenkov parameterization below 91 GeV and the Glauber-Gribov\n"
          << "parameterization above 91 GeV. n";
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......