G4ParticleHPContAngularPar.cc

Go to the documentation of this file.

//
// ********************************************************************
// * 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.          *
// ********************************************************************
//
// neutron_hp -- source file
// J.P. Wellisch, Nov-1996
// A prototype of the low energy neutron transport model.
//
// 09-May-06 fix in Sample by T. Koi
// 080318 Fix Compilation warnings - gcc-4.3.0 by T. Koi
//        (This fix has a real effect to the code.)
// 080409 Fix div0 error with G4FPE by T. Koi
// 080612 Fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #1
// 080714 Limiting the sum of energy of secondary particles by T. Koi
// 080801 Fix div0 error wiht G4FPE and memory leak by T. Koi
// 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
//
// P. Arce, June-2014 Conversion neutron_hp to particle_hp
//
// June-2019 - E. Mendoza --> redefinition of the residual mass to consider incident particles
// different than neutrons.
//
// V. Ivanchenko, July-2023 Basic revision of particle HP classes
//
 
#include "G4ParticleHPContAngularPar.hh"
 
#include "G4ParticleDefinition.hh"
#include "G4Alpha.hh"
#include "G4Deuteron.hh"
#include "G4Electron.hh"
#include "G4Gamma.hh"
#include "G4He3.hh"
#include "G4IonTable.hh"
#include "G4Neutron.hh"
#include "G4NucleiProperties.hh"
#include "G4ParticleHPKallbachMannSyst.hh"
#include "G4ParticleHPLegendreStore.hh"
#include "G4ParticleHPManager.hh"
#include "G4ParticleHPVector.hh"
#include "G4PhysicalConstants.hh"
#include "G4Positron.hh"
#include "G4Proton.hh"
#include "G4SystemOfUnits.hh"
#include "G4Triton.hh"
 
#include <set>
#include <vector>
 
G4ParticleHPContAngularPar::G4ParticleHPContAngularPar(const G4ParticleDefinition* p)
{
  theProjectile = (nullptr == p) ? G4Neutron::Neutron() : p;
  toBeCached v;
  fCache.Put(v);
  if (G4ParticleHPManager::GetInstance()->GetDoNotAdjustFinalState()) adjustResult = false;
}
 
G4ParticleHPContAngularPar::G4ParticleHPContAngularPar(G4ParticleHPContAngularPar& val)
{
  theEnergy = val.theEnergy;
  nEnergies = val.nEnergies;
  nDiscreteEnergies = val.nDiscreteEnergies;
  nAngularParameters = val.nAngularParameters;
  theProjectile = val.theProjectile;
  theManager = val.theManager;
  theInt = val.theInt;
  adjustResult = val.adjustResult;
  theMinEner = val.theMinEner;
  theMaxEner = val.theMaxEner;
  theEnergiesTransformed = val.theEnergiesTransformed;
  theDiscreteEnergies = val.theDiscreteEnergies;
  theDiscreteEnergiesOwn = val.theDiscreteEnergiesOwn;
  toBeCached v;
  fCache.Put(v);
  const std::size_t esize = nEnergies > 0 ? nEnergies : 1;
  theAngular = new G4ParticleHPList[esize];
  for (G4int ie = 0; ie < nEnergies; ++ie) {
    theAngular[ie].SetLabel(val.theAngular[ie].GetLabel());
    for (G4int ip = 0; ip < nAngularParameters; ++ip) {
      theAngular[ie].SetValue(ip, val.theAngular[ie].GetValue(ip));
    }
  }
}
 
G4ParticleHPContAngularPar::~G4ParticleHPContAngularPar()
{
  delete[] theAngular;
}
 
void G4ParticleHPContAngularPar::Init(std::istream& aDataFile, const G4ParticleDefinition* p)
{
  adjustResult = true;
  if (G4ParticleHPManager::GetInstance()->GetDoNotAdjustFinalState()) adjustResult = false;
 
  theProjectile = (nullptr == p) ? G4Neutron::Neutron() : p;
 
  aDataFile >> theEnergy >> nEnergies >> nDiscreteEnergies >> nAngularParameters;
  theEnergy *= eV;
  const std::size_t esize = nEnergies > 0 ? nEnergies : 1;
  theAngular = new G4ParticleHPList[esize];
  G4double sEnergy;
  for (G4int i = 0; i < nEnergies; ++i) {
    aDataFile >> sEnergy;
    sEnergy *= eV;
    theAngular[i].SetLabel(sEnergy);
    theAngular[i].Init(aDataFile, nAngularParameters, 1.);
    theMinEner = std::min(theMinEner, sEnergy);
    theMaxEner = std::max(theMaxEner, sEnergy);
  }
}
 
G4ReactionProduct* G4ParticleHPContAngularPar::Sample(G4double anEnergy, G4double massCode,
                                                      G4double /*targetMass*/, G4int angularRep,
                                                      G4int /*interpolE*/)
{
  // The following line is needed because it may change between runs by UI command
  adjustResult = true;
  if (G4ParticleHPManager::GetInstance()->GetDoNotAdjustFinalState()) adjustResult = false;
 
  auto result = new G4ReactionProduct;
  auto Z = static_cast<G4int>(massCode / 1000);
  auto A = static_cast<G4int>(massCode - 1000 * Z);
  if (massCode == 0) {
    result->SetDefinition(G4Gamma::Gamma());
  }
  else if (A == 0) {
    result->SetDefinition(G4Electron::Electron());
    if (Z == 1) result->SetDefinition(G4Positron::Positron());
  }
  else if (A == 1) {
    result->SetDefinition(G4Neutron::Neutron());
    if (Z == 1) result->SetDefinition(G4Proton::Proton());
  }
  else if (A == 2) {
    result->SetDefinition(G4Deuteron::Deuteron());
  }
  else if (A == 3) {
    result->SetDefinition(G4Triton::Triton());
    if (Z == 2) result->SetDefinition(G4He3::He3());
  }
  else if (A == 4) {
    result->SetDefinition(G4Alpha::Alpha());
    if (Z != 2)
      throw G4HadronicException(__FILE__, __LINE__,
                                "G4ParticleHPContAngularPar: Unknown ion case 1");
  }
  else {
    result->SetDefinition(G4IonTable::GetIonTable()->GetIon(Z, A, 0));
  }
 
  G4int i(0);
  G4int it(0);
  G4double fsEnergy(0);
  G4double cosTh(0);
  /*
  G4cout << "G4ParticleHPContAngularPar::Sample E=" << anEnergy <<" Z=" << Z << " A=" << A
         << " angularRep=" << angularRep << " Nd=" << nDiscreteEnergies 
         << " Ne=" << nEnergies << G4endl;
  */
  if (angularRep == 1) {
    if (nDiscreteEnergies != 0) {
      // 1st check remaining_energy
      // if this is the first set it. (How?)
      if (fCache.Get().fresh) {
        // Discrete Lines, larger energies come first
        // Continues Emssions, low to high                                 LAST
        fCache.Get().remaining_energy =
          std::max(theAngular[0].GetLabel(), theAngular[nEnergies - 1].GetLabel());
        fCache.Get().fresh = false;
      }
 
      // Cheating for small remaining_energy
      // Temporary solution
      if (nDiscreteEnergies == nEnergies) {
        fCache.Get().remaining_energy =
          std::max(fCache.Get().remaining_energy,
                   theAngular[nDiscreteEnergies - 1].GetLabel());  // Minimum Line
      }
      else {
        G4double cont_min = 0.0;
        for (G4int j = nDiscreteEnergies; j < nEnergies; ++j) {
          cont_min = theAngular[j].GetLabel();
          if (theAngular[j].GetValue(0) != 0.0) break;
        }
        fCache.Get().remaining_energy = std::max(
          fCache.Get().remaining_energy, std::min(theAngular[nDiscreteEnergies - 1].GetLabel(),
                                                   cont_min));  // Minimum Line or grid
      }
 
      G4double random = G4UniformRand();
      auto running = new G4double[nEnergies + 1];
      running[0] = 0.0;
 
      G4double delta;
      for (G4int j = 0; j < nDiscreteEnergies; ++j) {
        delta = 0.0;
        if (theAngular[j].GetLabel() <= fCache.Get().remaining_energy)
          delta = theAngular[j].GetValue(0);
        running[j + 1] = running[j] + delta;
      }
 
      G4double tot_prob_DIS = std::max(running[nDiscreteEnergies], 0.0);
 
      G4double delta1;
      for (G4int j = nDiscreteEnergies; j < nEnergies; ++j) {
        delta1 = 0.0;
        G4double e_low = 0.0;
        G4double e_high = 0.0;
        if (theAngular[j].GetLabel() <= fCache.Get().remaining_energy)
          delta1 = theAngular[j].GetValue(0);
 
        // To calculate Prob. e_low and e_high should be in eV
        // There are two cases:
        // 1: theAngular[nDiscreteEnergies].GetLabel() != 0.0
        //    delta1 should be used between j-1 and j
        //    At j = nDiscreteEnergies (the first) e_low should be set explicitly
        if (theAngular[j].GetLabel() != 0) {
          if (j == nDiscreteEnergies) {
            e_low = 0.0 / eV;
          }
          else {
            if (j < 1) j = 1;  // Protection against evaluation of arrays at index j-1
            e_low = theAngular[j - 1].GetLabel() / eV;
          }
          e_high = theAngular[j].GetLabel() / eV;
        }
 
        // 2: theAngular[nDiscreteEnergies].GetLabel() == 0.0
        //    delta1 should be used between j and j+1
        if (theAngular[j].GetLabel() == 0.0) {
          e_low = theAngular[j].GetLabel() / eV;
          if (j != nEnergies - 1) {
            e_high = theAngular[j + 1].GetLabel() / eV;
          }
          else {
            e_high = theAngular[j].GetLabel() / eV;
          }
        }
 
        running[j + 1] = running[j] + ((e_high - e_low) * delta1);
      }
      G4double tot_prob_CON = std::max(running[nEnergies] - running[nDiscreteEnergies], 0.0);
 
      // Give up in the pathological case of null probabilities
      if (tot_prob_DIS == 0.0 && tot_prob_CON == 0.0) {
        delete[] running;
    return result;
      }
      // Normalize random
      random *= (tot_prob_DIS + tot_prob_CON);
      // 2nd Judge Discrete or not
 
      // This should be relatively close to 1  For safty
      if (random <= (tot_prob_DIS / (tot_prob_DIS + tot_prob_CON))
          || nDiscreteEnergies == nEnergies)
      {
        // Discrete Emission
        for (G4int j = 0; j < nDiscreteEnergies; ++j) {
          // Here we should use i+1
          if (random < running[j + 1]) {
            it = j;
            break;
          }
        }
        fsEnergy = theAngular[it].GetLabel();
 
        G4ParticleHPLegendreStore theStore(1);
        theStore.Init(0, fsEnergy, nAngularParameters);
        for (G4int j = 0; j < nAngularParameters; ++j) {
          theStore.SetCoeff(0, j, theAngular[it].GetValue(j));
        }
        // use it to sample.
        cosTh = theStore.SampleMax(fsEnergy);
        // Done
      }
      else {
        // Continuous emission
        for (G4int j = nDiscreteEnergies; j < nEnergies; ++j) {
          // Here we should use i
          if (random < running[j]) {
            it = j;
            break;
          }
        }
 
        if (it < 1) it = 1;  // Protection against evaluation of arrays at index it-1
 
        G4double x1 = running[it - 1];
        G4double x2 = running[it];
 
        G4double y1 = 0.0;
        if (it != nDiscreteEnergies) y1 = theAngular[it - 1].GetLabel();
        G4double y2 = theAngular[it].GetLabel();
 
        fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it), random, x1, x2, y1, y2);
 
        G4ParticleHPLegendreStore theStore(2);
        theStore.Init(0, y1, nAngularParameters);
        theStore.Init(1, y2, nAngularParameters);
        theStore.SetManager(theManager);
        G4int itt;
        for (G4int j = 0; j < nAngularParameters; ++j) {
          itt = it;
          if (it == nDiscreteEnergies) itt = it + 1;
          // "This case "it-1" has data for Discrete, so we will use an extrpolated values it and
          // it+1
          theStore.SetCoeff(0, j, theAngular[itt - 1].GetValue(j));
          theStore.SetCoeff(1, j, theAngular[itt].GetValue(j));
        }
        // use it to sample.
        cosTh = theStore.SampleMax(fsEnergy);
 
        // Done
      }
 
      // The remaining energy needs to be lowered by the photon energy in *any* case.
      // Otherwise additional photons with too high energy will be produced - therefore the
      // adjustResult condition has been removed
      fCache.Get().remaining_energy -= fsEnergy;
      delete[] running;
 
      // end (nDiscreteEnergies != 0) branch
    }
    else {
      // Only continue, TK will clean up
      if (fCache.Get().fresh) {
        fCache.Get().remaining_energy = theAngular[nEnergies - 1].GetLabel();
        fCache.Get().fresh = false;
      }
 
      G4double random = G4UniformRand();
      auto running = new G4double[nEnergies];
      running[0] = 0;
      G4double weighted = 0;
      for (i = 1; i < nEnergies; i++) {
        running[i] = running[i - 1];
        if (fCache.Get().remaining_energy >= theAngular[i].GetLabel()) {
          running[i] += theInt.GetBinIntegral(
            theManager.GetScheme(i - 1), theAngular[i - 1].GetLabel(), theAngular[i].GetLabel(),
            theAngular[i - 1].GetValue(0), theAngular[i].GetValue(0));
          weighted += theInt.GetWeightedBinIntegral(
            theManager.GetScheme(i - 1), theAngular[i - 1].GetLabel(), theAngular[i].GetLabel(),
            theAngular[i - 1].GetValue(0), theAngular[i].GetValue(0));
        }
      }
 
      // Cache the mean energy in this distribution
      if (nEnergies == 1 || running[nEnergies - 1] == 0) {
        fCache.Get().currentMeanEnergy = 0.0;
      }
      else {
        fCache.Get().currentMeanEnergy = weighted / running[nEnergies - 1];
      }
 
      if (nEnergies == 1) it = 0;
      if (running[nEnergies - 1] != 0) {
        for (i = 1; i < nEnergies; i++) {
          it = i;
          if (random < running[i] / running[nEnergies - 1]) break;
        }
      }
 
      if (running[nEnergies - 1] == 0) it = 0;
      if (it < nDiscreteEnergies || it == 0) {
        if (it == 0) {
          fsEnergy = theAngular[it].GetLabel();
          G4ParticleHPLegendreStore theStore(1);
          theStore.Init(0, fsEnergy, nAngularParameters);
          for (i = 0; i < nAngularParameters; i++) {
            theStore.SetCoeff(0, i, theAngular[it].GetValue(i));
          }
          // use it to sample.
          cosTh = theStore.SampleMax(fsEnergy);
        }
        else {
          G4double e1, e2;
          e1 = theAngular[it - 1].GetLabel();
          e2 = theAngular[it].GetLabel();
          fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it), random,
                                        running[it - 1] / running[nEnergies - 1],
                                        running[it] / running[nEnergies - 1], e1, e2);
          // fill a Legendrestore
          G4ParticleHPLegendreStore theStore(2);
          theStore.Init(0, e1, nAngularParameters);
          theStore.Init(1, e2, nAngularParameters);
          for (i = 0; i < nAngularParameters; i++) {
            theStore.SetCoeff(0, i, theAngular[it - 1].GetValue(i));
            theStore.SetCoeff(1, i, theAngular[it].GetValue(i));
          }
          // use it to sample.
          theStore.SetManager(theManager);
          cosTh = theStore.SampleMax(fsEnergy);
        }
      }
      else {  // continuum contribution
        G4double x1 = running[it - 1] / running[nEnergies - 1];
        G4double x2 = running[it] / running[nEnergies - 1];
        G4double y1 = theAngular[it - 1].GetLabel();
        G4double y2 = theAngular[it].GetLabel();
        fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it), random, x1, x2, y1, y2);
        G4ParticleHPLegendreStore theStore(2);
        theStore.Init(0, y1, nAngularParameters);
        theStore.Init(1, y2, nAngularParameters);
        theStore.SetManager(theManager);
        for (i = 0; i < nAngularParameters; i++) {
          theStore.SetCoeff(0, i, theAngular[it - 1].GetValue(i));
          theStore.SetCoeff(1, i, theAngular[it].GetValue(i));
        }
        // use it to sample.
        cosTh = theStore.SampleMax(fsEnergy);
      }
      delete[] running;
 
      // The remaining energy needs to be lowered by the photon energy in
      // *any* case.  Otherwise additional photons with too much energy will be
      // produced - therefore the  adjustResult condition has been removed
 
      fCache.Get().remaining_energy -= fsEnergy;
      // end if (nDiscreteEnergies != 0)
    }
    // end of (angularRep == 1) branch
  }
  else if (angularRep == 2) {
    // first get the energy (already the right for this incoming energy)
    G4int j;
    auto running = new G4double[nEnergies];
    running[0] = 0;
    G4double weighted = 0;
    for (j = 1; j < nEnergies; ++j) {
      if (j != 0) running[j] = running[j - 1];
      running[j] += theInt.GetBinIntegral(theManager.GetScheme(j - 1), theAngular[j - 1].GetLabel(),
                                          theAngular[j].GetLabel(), theAngular[j - 1].GetValue(0),
                                          theAngular[j].GetValue(0));
      weighted += theInt.GetWeightedBinIntegral(
        theManager.GetScheme(j - 1), theAngular[j - 1].GetLabel(), theAngular[j].GetLabel(),
        theAngular[j - 1].GetValue(0), theAngular[j].GetValue(0));
    }
 
    // Cache the mean energy in this distribution
    if (nEnergies == 1)
      fCache.Get().currentMeanEnergy = 0.0;
    else
      fCache.Get().currentMeanEnergy = weighted / running[nEnergies - 1];
 
    G4int itt(0);
    G4double randkal = G4UniformRand();
    for (j = 1; j < nEnergies; ++j) {
      itt = j;
      if (randkal*running[nEnergies - 1] < running[j]) break;
    }
 
    // Interpolate the secondary energy
    G4double x, x1, x2, y1, y2;
    if (itt == 0) itt = 1;
    x = randkal * running[nEnergies - 1];
    x1 = running[itt - 1];
    x2 = running[itt];
    G4double compoundFraction;
    // interpolate energy
    y1 = theAngular[itt - 1].GetLabel();
    y2 = theAngular[itt].GetLabel();
    fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(itt - 1), x, x1, x2, y1, y2);
 
    // For theta, interpolate the compoundFractions
    G4double cLow = theAngular[itt - 1].GetValue(1);
    G4double cHigh = theAngular[itt].GetValue(1);
    compoundFraction = theInt.Interpolate(theManager.GetScheme(itt), fsEnergy, y1, y2, cLow, cHigh);
 
    if (compoundFraction > 1.0)
      compoundFraction = 1.0;  // Protection against unphysical interpolation
 
    delete[] running;
 
    // get cosTh
    G4double incidentEnergy = anEnergy;
    G4double incidentMass = theProjectile->GetPDGMass();
    G4double productEnergy = fsEnergy;
    G4double productMass = result->GetMass();
    auto targetZ = G4int(fCache.Get().theTargetCode / 1000);
    auto targetA = G4int(fCache.Get().theTargetCode - 1000 * targetZ);
 
    // To correspond to natural composition (-nat-) data files.
    if (targetA == 0) targetA = G4int(fCache.Get().theTarget->GetMass() / amu_c2 + 0.5);
    G4double targetMass = fCache.Get().theTarget->GetMass();
    auto incidentA = G4int(incidentMass / amu_c2 + 0.5);
    auto incidentZ = G4int(theProjectile->GetPDGCharge() + 0.5);
    G4int residualA = targetA + incidentA - A;
    G4int residualZ = targetZ + incidentZ - Z;
    G4double residualMass = G4NucleiProperties::GetNuclearMass(residualA, residualZ);
 
    G4ParticleHPKallbachMannSyst theKallbach(
      compoundFraction, incidentEnergy, incidentMass, productEnergy, productMass, residualMass,
      residualA, residualZ, targetMass, targetA, targetZ, incidentA, incidentZ, A, Z);
    cosTh = theKallbach.Sample(anEnergy);
    // end (angularRep == 2) branch
  }
  else if (angularRep > 10 && angularRep < 16) {
    G4double random = G4UniformRand();
    auto running = new G4double[nEnergies];
    running[0] = 0;
    G4double weighted = 0;
    for (i = 1; i < nEnergies; ++i) {
      if (i != 0) running[i] = running[i - 1];
      running[i] += theInt.GetBinIntegral(theManager.GetScheme(i - 1), theAngular[i - 1].GetLabel(),
                                          theAngular[i].GetLabel(), theAngular[i - 1].GetValue(0),
                                          theAngular[i].GetValue(0));
      weighted += theInt.GetWeightedBinIntegral(
        theManager.GetScheme(i - 1), theAngular[i - 1].GetLabel(), theAngular[i].GetLabel(),
        theAngular[i - 1].GetValue(0), theAngular[i].GetValue(0));
    }
 
    // Cache the mean energy in this distribution
    if (nEnergies == 1)
      fCache.Get().currentMeanEnergy = 0.0;
    else
      fCache.Get().currentMeanEnergy = weighted / running[nEnergies - 1];
 
    if (nEnergies == 1) it = 0;
    for (i = 1; i < nEnergies; i++) {
      it = i;
      if (random < running[i] / running[nEnergies - 1]) break;
    }
 
    if (it < nDiscreteEnergies || it == 0) {
      if (it == 0) {
        fsEnergy = theAngular[0].GetLabel();
        G4ParticleHPVector theStore;
        G4int aCounter = 0;
        for (G4int j = 1; j < nAngularParameters; j += 2) {
          theStore.SetX(aCounter, theAngular[0].GetValue(j));
          theStore.SetY(aCounter, theAngular[0].GetValue(j + 1));
          aCounter++;
        }
        G4InterpolationManager aMan;
        aMan.Init(angularRep - 10, nAngularParameters - 1);
        theStore.SetInterpolationManager(aMan);
        cosTh = theStore.Sample();
      }
      else {
        fsEnergy = theAngular[it].GetLabel();
        G4ParticleHPVector theStore;
        G4InterpolationManager aMan;
        aMan.Init(angularRep - 10, nAngularParameters - 1);
        theStore.SetInterpolationManager(aMan);  // Store interpolates f(costh)
        G4InterpolationScheme currentScheme = theManager.GetInverseScheme(it);
        G4int aCounter = 0;
        for (G4int j = 1; j < nAngularParameters; j += 2) {
          theStore.SetX(aCounter, theAngular[it].GetValue(j));
          theStore.SetY(aCounter, theInt.Interpolate(currentScheme, random,
                                                     running[it - 1] / running[nEnergies - 1],
                                                     running[it] / running[nEnergies - 1],
                                                     theAngular[it - 1].GetValue(j + 1),
                                                     theAngular[it].GetValue(j + 1)));
          ++aCounter;
        }
        cosTh = theStore.Sample();
      }
    }
    else {
      G4double x1 = running[it - 1] / running[nEnergies - 1];
      G4double x2 = running[it] / running[nEnergies - 1];
      G4double y1 = theAngular[it - 1].GetLabel();
      G4double y2 = theAngular[it].GetLabel();
      fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it), random, x1, x2, y1, y2);
      G4ParticleHPVector theBuff1;
      G4ParticleHPVector theBuff2;
      G4InterpolationManager aMan;
      aMan.Init(angularRep - 10, nAngularParameters - 1);
 
      G4int j;
      for (i = 0, j = 1; i < nAngularParameters; i++, j += 2) {
        theBuff1.SetX(i, theAngular[it - 1].GetValue(j));
        theBuff1.SetY(i, theAngular[it - 1].GetValue(j + 1));
        theBuff2.SetX(i, theAngular[it].GetValue(j));
        theBuff2.SetY(i, theAngular[it].GetValue(j + 1));
      }
 
      G4ParticleHPVector theStore;
      theStore.SetInterpolationManager(aMan);  // Store interpolates f(costh)
      x1 = y1;
      x2 = y2;
      G4double x, y;
      for (i = 0; i < theBuff1.GetVectorLength(); i++) {
        x = theBuff1.GetX(i);  // costh binning identical
        y1 = theBuff1.GetY(i);
        y2 = theBuff2.GetY(i);
        y = theInt.Interpolate(theManager.GetScheme(it), fsEnergy, theAngular[it - 1].GetLabel(),
                               theAngular[it].GetLabel(), y1, y2);
        theStore.SetX(i, x);
        theStore.SetY(i, y);
      }
      cosTh = theStore.Sample();
    }
    delete[] running;
  }
  else {
    throw G4HadronicException(__FILE__, __LINE__,
                              "G4ParticleHPContAngularPar::Sample: Unknown angular representation");
  }
  //G4cout << "  Efin=" << fsEnergy << G4endl;
  result->SetKineticEnergy(fsEnergy);
 
  G4double phi = twopi * G4UniformRand();
  if(cosTh > 1.0) { cosTh = 1.0; }
  else if (cosTh < -1.0) { cosTh = -1.0; }
  G4double sinth = std::sqrt((1.0 - cosTh)*(1.0 + cosTh));
  G4double mtot = result->GetTotalMomentum();
  G4ThreeVector tempVector(mtot * sinth * std::cos(phi), mtot * sinth * std::sin(phi), mtot * cosTh);
  result->SetMomentum(tempVector);
  return result;
}
 
void G4ParticleHPContAngularPar::PrepareTableInterpolation()
{
  // Discrete energies: store own energies in a map for faster searching
  //
  // The data files sometimes have identical discrete energies (likely typos)
  // which would lead to overwriting the already existing index and hence
  // creating a hole in the lookup table.
  // No attempt is made here to correct for the energies - rather an epsilon
  // is subtracted from the energy in order to uniquely identify the line
 
  for (G4int ie = 0; ie < nDiscreteEnergies; ie++) {
    // check if energy is already present and subtract epsilon if that's the case
    G4double myE = theAngular[ie].GetLabel();
    while (theDiscreteEnergiesOwn.find(myE) != theDiscreteEnergiesOwn.end()) {
      myE -= 1e-6;
    }
    theDiscreteEnergiesOwn[myE] = ie;
  }
  return;
}
 
void G4ParticleHPContAngularPar::BuildByInterpolation(G4double anEnergy,
                                                      G4InterpolationScheme aScheme,
                                                      G4ParticleHPContAngularPar& angpar1,
                                                      G4ParticleHPContAngularPar& angpar2)
{
  G4int ie, ie1, ie2, ie1Prev, ie2Prev;
  // Only rebuild the interpolation table if there is a new interaction.
  // For several subsequent samplings of final state particles in the same
  // interaction the existing table should be used
  if (!fCache.Get().fresh) return;
 
  // Make copies of angpar1 and angpar2. Since these are given by reference
  // it can not be excluded that one of them is "this". Hence this code uses
  // potentially the old "this" for creating the new this - which leads to
  // memory corruption if the old is not stored as separarte object for lookup
  const G4ParticleHPContAngularPar copyAngpar1(angpar1), copyAngpar2(angpar2);
 
  nAngularParameters = copyAngpar1.nAngularParameters;
  theManager = copyAngpar1.theManager;
  theEnergy = anEnergy;
  theMinEner = DBL_MAX;  // min and max will be re-calculated after interpolation
  theMaxEner = -DBL_MAX;
 
  // The two discrete sets must be merged. A vector holds the temporary data to
  // be copied to the array in the end.  Since the G4ParticleHPList class
  // contains pointers, can't simply assign elements of this type. Each member
  // needs to call the explicit Set() method instead.
 
  // First, average probabilities for those lines that are in both sets
  const std::map<G4double, G4int> discEnerOwn1 = copyAngpar1.GetDiscreteEnergiesOwn();
  const std::map<G4double, G4int> discEnerOwn2 = copyAngpar2.GetDiscreteEnergiesOwn();
  std::map<G4double, G4int>::const_iterator itedeo1;
  std::map<G4double, G4int>::const_iterator itedeo2;
  std::vector<G4ParticleHPList*> vAngular(discEnerOwn1.size());
  G4double discEner1;
  for (itedeo1 = discEnerOwn1.cbegin(); itedeo1 != discEnerOwn1.cend(); ++itedeo1) {
    discEner1 = itedeo1->first;
    if (discEner1 < theMinEner) {
      theMinEner = discEner1;
    }
    if (discEner1 > theMaxEner) {
      theMaxEner = discEner1;
    }
    ie1 = itedeo1->second;
    itedeo2 = discEnerOwn2.find(discEner1);
    if (itedeo2 == discEnerOwn2.cend()) {
      ie2 = -1;
    }
    else {
      ie2 = itedeo2->second;
    }
    vAngular[ie1] = new G4ParticleHPList();
    vAngular[ie1]->SetLabel(copyAngpar1.theAngular[ie1].GetLabel());
    G4double val1, val2;
    for (G4int ip = 0; ip < nAngularParameters; ++ip) {
      val1 = copyAngpar1.theAngular[ie1].GetValue(ip);
      if (ie2 != -1) {
        val2 = copyAngpar2.theAngular[ie2].GetValue(ip);
      }
      else {
        val2 = 0.;
      }
      G4double value = theInt.Interpolate(aScheme, anEnergy, copyAngpar1.theEnergy,
                                          copyAngpar2.theEnergy, val1, val2);
      vAngular[ie1]->SetValue(ip, value);
    }
  }  // itedeo1 loop
 
  // Add the ones in set2 but not in set1
  std::vector<G4ParticleHPList*>::const_iterator itv;
  G4double discEner2;
  for (itedeo2 = discEnerOwn2.cbegin(); itedeo2 != discEnerOwn2.cend(); ++itedeo2) {
    discEner2 = itedeo2->first;
    ie2 = itedeo2->second;
    G4bool notFound = true;
    itedeo1 = discEnerOwn1.find(discEner2);
    if (itedeo1 != discEnerOwn1.cend()) {
      notFound = false;
    }
    if (notFound) {
      // not yet in list
      if (discEner2 < theMinEner) {
        theMinEner = discEner2;
      }
      if (discEner2 > theMaxEner) {
        theMaxEner = discEner2;
      }
      // find position to insert
      G4bool isInserted = false;
      ie = 0;
      for (itv = vAngular.cbegin(); itv != vAngular.cend(); ++itv, ++ie) {
        if (discEner2 > (*itv)->GetLabel()) {
          itv = vAngular.insert(itv, new G4ParticleHPList);
          (*itv)->SetLabel(copyAngpar2.theAngular[ie2].GetLabel());
          isInserted = true;
          break;
        }
      }
      if (!isInserted) {
        ie = (G4int)vAngular.size();
        vAngular.push_back(new G4ParticleHPList);
        vAngular[ie]->SetLabel(copyAngpar2.theAngular[ie2].GetLabel());
        isInserted = true;
      }
 
      G4double val1, val2;
      for (G4int ip = 0; ip < nAngularParameters; ++ip) {
        val1 = 0;
        val2 = copyAngpar2.theAngular[ie2].GetValue(ip);
        G4double value = theInt.Interpolate(aScheme, anEnergy, copyAngpar1.theEnergy,
                                            copyAngpar2.theEnergy, val1, val2);
        vAngular[ie]->SetValue(ip, value);
      }
    }  // end if(notFound)
  }  // end loop on itedeo2
 
  // Store new discrete list
  nDiscreteEnergies = (G4int)vAngular.size();
  delete[] theAngular;
  theAngular = nullptr;
  if (nDiscreteEnergies > 0) {
    theAngular = new G4ParticleHPList[nDiscreteEnergies];
  }
  theDiscreteEnergiesOwn.clear();
  theDiscreteEnergies.clear();
  for (ie = 0; ie < nDiscreteEnergies; ++ie) {
    theAngular[ie].SetLabel(vAngular[ie]->GetLabel());
    for (G4int ip = 0; ip < nAngularParameters; ++ip) {
      theAngular[ie].SetValue(ip, vAngular[ie]->GetValue(ip));
    }
    theDiscreteEnergiesOwn[theAngular[ie].GetLabel()] = ie;
    theDiscreteEnergies.insert(theAngular[ie].GetLabel());
  }
 
  // The continuous energies need to be made from scratch like the discrete
  // ones. Therefore the re-assignemnt of theAngular needs to be done
  // after the continuous energy set is also finalized. Only then the
  // total number of nEnergies is known and the array can be allocated.
 
  // Get minimum and maximum energy interpolating
  // Don't use theMinEner or theMaxEner here, since the transformed energies
  // need the interpolated range from the original Angpar
  G4double interMinEner = copyAngpar1.GetMinEner()
                          + (theEnergy - copyAngpar1.GetEnergy())
                              * (copyAngpar2.GetMinEner() - copyAngpar1.GetMinEner())
                              / (copyAngpar2.GetEnergy() - copyAngpar1.GetEnergy());
  G4double interMaxEner = copyAngpar1.GetMaxEner()
                          + (theEnergy - copyAngpar1.GetEnergy())
                              * (copyAngpar2.GetMaxEner() - copyAngpar1.GetMaxEner())
                              / (copyAngpar2.GetEnergy() - copyAngpar1.GetEnergy());
 
  // Loop to energies of new set
  theEnergiesTransformed.clear();
 
  G4int nEnergies1 = copyAngpar1.GetNEnergies();
  G4int nDiscreteEnergies1 = copyAngpar1.GetNDiscreteEnergies();
  G4double minEner1 = copyAngpar1.GetMinEner();
  G4double maxEner1 = copyAngpar1.GetMaxEner();
  G4int nEnergies2 = copyAngpar2.GetNEnergies();
  G4int nDiscreteEnergies2 = copyAngpar2.GetNDiscreteEnergies();
  G4double minEner2 = copyAngpar2.GetMinEner();
  G4double maxEner2 = copyAngpar2.GetMaxEner();
 
  // First build the list of transformed energies normalized
  // to the new min max by assuming that the min-max range of
  // each set would be scalable to the new, interpolated min
  // max range
 
  G4double e1(0.);
  G4double eTNorm1(0.);
  for (ie1 = nDiscreteEnergies1; ie1 < nEnergies1; ++ie1) {
    e1 = copyAngpar1.theAngular[ie1].GetLabel();
    eTNorm1 = (e1 - minEner1);
    if (maxEner1 != minEner1) eTNorm1 /= (maxEner1 - minEner1);
    if (eTNorm1 >= 0 && eTNorm1 <= 1) theEnergiesTransformed.insert(eTNorm1);
  }
 
  G4double e2(0.);
  G4double eTNorm2(0.);
  for (ie2 = nDiscreteEnergies2; ie2 < nEnergies2; ++ie2) {
    e2 = copyAngpar2.theAngular[ie2].GetLabel();
    eTNorm2 = (e2 - minEner2);
    if (maxEner2 != minEner2) eTNorm2 /= (maxEner2 - minEner2);
    if (eTNorm2 >= 0 && eTNorm2 <= 1) theEnergiesTransformed.insert(eTNorm2);
  }
 
  // Now the list of energies is complete
  nEnergies = nDiscreteEnergies + (G4int)theEnergiesTransformed.size();
 
  // Create final array of angular parameters
  const std::size_t esize = nEnergies > 0 ? nEnergies : 1;
  auto theNewAngular = new G4ParticleHPList[esize];
 
  // Copy discrete energies and interpolated parameters to new array
 
  if (theAngular != nullptr) {
    for (ie = 0; ie < nDiscreteEnergies; ++ie) {
      theNewAngular[ie].SetLabel(theAngular[ie].GetLabel());
      for (G4int ip = 0; ip < nAngularParameters; ++ip) {
        theNewAngular[ie].SetValue(ip, theAngular[ie].GetValue(ip));
      }
    }
    delete[] theAngular;
  }
  theAngular = theNewAngular;
 
  // Interpolate the continuous energies for new array
  auto iteet = theEnergiesTransformed.begin();
 
  G4double e1Interp(0.);
  G4double e2Interp(0.);
  for (ie = nDiscreteEnergies; ie < nEnergies; ++ie, ++iteet) {
    G4double eT = (*iteet);
 
    //--- Use eT1 = eT: Get energy and parameters of copyAngpar1 for this eT
    e1Interp = (maxEner1 - minEner1) * eT + minEner1;
    //----- Get parameter value corresponding to this e1Interp
    for (ie1 = nDiscreteEnergies1; ie1 < nEnergies1; ++ie1) {
      if ((copyAngpar1.theAngular[ie1].GetLabel() - e1Interp) > 1.E-10 * e1Interp) break;
    }
    ie1Prev = ie1 - 1;
    if (ie1 == 0) ++ie1Prev;
    if (ie1 == nEnergies1) {
      ie1--;
      ie1Prev = ie1;
    }
 
    //--- Use eT2 = eT: Get energy and parameters of copyAngpar2 for this eT
    e2Interp = (maxEner2 - minEner2) * eT + minEner2;
    //----- Get parameter value corresponding to this e2Interp
    for (ie2 = nDiscreteEnergies2; ie2 < nEnergies2; ++ie2) {
      if ((copyAngpar2.theAngular[ie2].GetLabel() - e2Interp) > 1.E-10 * e2Interp) break;
    }
    ie2Prev = ie2 - 1;
    if (ie2 == 0) ++ie2Prev;
    if (ie2 == nEnergies2) {
      ie2--;
      ie2Prev = ie2;
    }
 
    //---- Energy corresponding to energy transformed
    G4double eN = (interMaxEner - interMinEner) * eT + interMinEner;
 
    theAngular[ie].SetLabel(eN);
    if (eN < theMinEner) {
      theMinEner = eN;
    }
    if (eN > theMaxEner) {
      theMaxEner = eN;
    }
 
    G4double val1(0.);
    G4double val2(0.);
    G4double value(0.);
    for (G4int ip = 0; ip < nAngularParameters; ++ip) {
      val1 = theInt.Interpolate2(
               theManager.GetScheme(ie), e1Interp, copyAngpar1.theAngular[ie1Prev].GetLabel(),
               copyAngpar1.theAngular[ie1].GetLabel(), copyAngpar1.theAngular[ie1Prev].GetValue(ip),
               copyAngpar1.theAngular[ie1].GetValue(ip))
             * (maxEner1 - minEner1);
      val2 = theInt.Interpolate2(
               theManager.GetScheme(ie), e2Interp, copyAngpar2.theAngular[ie2Prev].GetLabel(),
               copyAngpar2.theAngular[ie2].GetLabel(), copyAngpar2.theAngular[ie2Prev].GetValue(ip),
               copyAngpar2.theAngular[ie2].GetValue(ip))
             * (maxEner2 - minEner2);
 
      value = theInt.Interpolate(aScheme, anEnergy, copyAngpar1.theEnergy, copyAngpar2.theEnergy,
                                 val1, val2);
      if (interMaxEner != interMinEner) {
        value /= (interMaxEner - interMinEner);
      }
      else if (value != 0) {
        throw G4HadronicException(__FILE__, __LINE__,
                                  "G4ParticleHPContAngularPar::PrepareTableInterpolation "
                                  "interMaxEner == interMinEner and  value != 0.");
      }
      theAngular[ie].SetValue(ip, value);
    }
  }  // end loop on nDiscreteEnergies
 
  for (itv = vAngular.cbegin(); itv != vAngular.cend(); ++itv)
    delete (*itv);
}
 
void G4ParticleHPContAngularPar::Dump() const
{
  G4cout << theEnergy << " " << nEnergies << " " << nDiscreteEnergies << " " << nAngularParameters
         << G4endl;
 
  for (G4int ii = 0; ii < nEnergies; ++ii)
    theAngular[ii].Dump();
}