#include <G4ParticleHPInelasticCompFS.hh>

Inheritance diagram for G4ParticleHPInelasticCompFS:

Public Member Functions
	G4ParticleHPInelasticCompFS ()

	~G4ParticleHPInelasticCompFS () override

void	Init (G4double A, G4double Z, G4int M, G4String &dirName, G4String &aSFType, G4ParticleDefinition *) override

void	InitGammas (G4double AR, G4double ZR)

G4HadFinalState *	ApplyYourself (const G4HadProjectile &theTrack) override=0

G4ParticleHPFinalState *	New () override=0

G4double	GetXsec (G4double anEnergy) override

G4ParticleHPVector *	GetXsec () override

G4int	SelectExitChannel (G4double eKinetic)

void	CompositeApply (const G4HadProjectile &theTrack, G4ParticleDefinition *aHadron)

void	InitDistributionInitialState (G4ReactionProduct &anIncidentPart, G4ReactionProduct &aTarget, G4int it)

	G4ParticleHPInelasticCompFS (G4ParticleHPInelasticCompFS &)=delete

G4ParticleHPInelasticCompFS &	operator= (const G4ParticleHPInelasticCompFS &right)=delete

Public Member Functions inherited from G4ParticleHPFinalState
	G4ParticleHPFinalState ()

virtual	~G4ParticleHPFinalState ()

void	Init (G4double A, G4double Z, G4String &dirName, G4String &aFSType, G4ParticleDefinition *p)

G4bool	HasXsec ()

G4bool	HasFSData ()

G4bool	HasAnyData ()

void	SetA_Z (G4double anA, G4double aZ, G4int aM=0)

G4double	GetZ ()

G4double	GetN ()

G4double	GetA ()

G4int	GetM ()

void	SetAZMs (G4ParticleHPDataUsed used)

void	SetAZMs (G4double anA, G4double aZ, G4int aM, G4ParticleHPDataUsed used)

void	SetProjectile (G4ParticleDefinition *projectile)

G4ParticleHPFinalState &	operator= (const G4ParticleHPFinalState &right)=delete

	G4ParticleHPFinalState (const G4ParticleHPFinalState &)=delete

Protected Attributes
G4ParticleHPVector *	theXsection [51]

G4ParticleHPEnergyDistribution *	theEnergyDistribution [51]

G4ParticleHPAngular *	theAngularDistribution [51]

G4ParticleHPEnAngCorrelation *	theEnergyAngData [51]

G4ParticleHPPhotonDist *	theFinalStatePhotons [51]

std::vector< G4double >	QI

std::vector< G4int >	LR

G4ParticleHPDeExGammas	theGammas

G4String	gammaPath

Protected Attributes inherited from G4ParticleHPFinalState
G4ParticleDefinition *	theProjectile {nullptr}

G4ParticleHPManager *	fManager

G4IonTable *	ionTable

G4int	theBaseA {0}

G4int	theBaseZ {0}

G4int	theBaseM {0}

G4int	theNDLDataZ {0}

G4int	theNDLDataA {0}

G4int	theNDLDataM {0}

G4int	secID {-1}

G4bool	hasXsec {true}

G4bool	hasFSData {true}

G4bool	hasAnyData {true}

G4ParticleHPNames	theNames

G4Cache< G4HadFinalState * >	theResult

Additional Inherited Members
Protected Member Functions inherited from G4ParticleHPFinalState
void	adjust_final_state (G4LorentzVector)

Detailed Description

Definition at line 51 of file G4ParticleHPInelasticCompFS.hh.

Constructor & Destructor Documentation

◆ G4ParticleHPInelasticCompFS() [1/2]

G4ParticleHPInelasticCompFS::G4ParticleHPInelasticCompFS ( )

Definition at line 72 of file G4ParticleHPInelasticCompFS.cc.

  : G4ParticleHPFinalState()
{
  QI.resize(51);
  LR.resize(51);
  for (G4int i = 0; i < 51; ++i) {
    hasXsec = true;
    theXsection[i] = nullptr;
    theEnergyDistribution[i] = nullptr;
    theAngularDistribution[i] = nullptr;
    theEnergyAngData[i] = nullptr;
    theFinalStatePhotons[i] = nullptr;
    QI[i] = 0.0;
    LR[i] = 0;
  }
}

◆ ~G4ParticleHPInelasticCompFS()

G4ParticleHPInelasticCompFS::~G4ParticleHPInelasticCompFS ( )

override

Definition at line 89 of file G4ParticleHPInelasticCompFS.cc.

{
  for (G4int i = 0; i < 51; ++i) {
    if (theXsection[i] != nullptr) delete theXsection[i];
    if (theEnergyDistribution[i] != nullptr) delete theEnergyDistribution[i];
    if (theAngularDistribution[i] != nullptr) delete theAngularDistribution[i];
    if (theEnergyAngData[i] != nullptr) delete theEnergyAngData[i];
    if (theFinalStatePhotons[i] != nullptr) delete theFinalStatePhotons[i];
  }
}

◆ G4ParticleHPInelasticCompFS() [2/2]

G4ParticleHPInelasticCompFS::G4ParticleHPInelasticCompFS ( G4ParticleHPInelasticCompFS & )

delete

Member Function Documentation

◆ ApplyYourself()

G4HadFinalState * G4ParticleHPInelasticCompFS::ApplyYourself ( const G4HadProjectile & theTrack )

overridepure virtual

Reimplemented from G4ParticleHPFinalState.

Implemented in G4ParticleHPAInelasticFS, G4ParticleHPDInelasticFS, G4ParticleHPHe3InelasticFS, G4ParticleHPNInelasticFS, G4ParticleHPPInelasticFS, and G4ParticleHPTInelasticFS.

◆ CompositeApply()

void G4ParticleHPInelasticCompFS::CompositeApply	(	const G4HadProjectile &	theTrack,
		G4ParticleDefinition *	aHadron )

Definition at line 254 of file G4ParticleHPInelasticCompFS.cc.

{
  // prepare neutron
  if (theResult.Get() == nullptr) theResult.Put(new G4HadFinalState);
  theResult.Get()->Clear();
  G4double eKinetic = theTrack.GetKineticEnergy();
  const G4HadProjectile* hadProjectile = &theTrack;
  G4ReactionProduct incidReactionProduct(hadProjectile->GetDefinition());
  incidReactionProduct.SetMomentum(hadProjectile->Get4Momentum().vect());
  incidReactionProduct.SetKineticEnergy(eKinetic);
 
  // prepare target
  for (G4int i = 0; i < 50; ++i) {
    if (theXsection[i] != nullptr) {
      break;
    }
  }
 
  G4double targetMass = G4NucleiProperties::GetNuclearMass(theBaseA, theBaseZ);
#ifdef G4VERBOSE
  if (fManager->GetDEBUG())
    G4cout << "G4ParticleHPInelasticCompFS::CompositeApply A=" << theBaseA << " Z="
       << theBaseZ << " incident " << hadProjectile->GetDefinition()->GetParticleName()
       << G4endl;
#endif
  G4ReactionProduct theTarget;
  G4Nucleus aNucleus;
  // G4ThreeVector neuVelo =
  // (1./hadProjectile->GetDefinition()->GetPDGMass())*incidReactionProduct.GetMomentum(); theTarget
  // = aNucleus.GetBiasedThermalNucleus( targetMass/hadProjectile->GetDefinition()->GetPDGMass() ,
  // neuVelo, theTrack.GetMaterial()->GetTemperature()); G4Nucleus::GetBiasedThermalNucleus requests
  // normalization of mass and velocity in neutron mass
  G4ThreeVector neuVelo = incidReactionProduct.GetMomentum() / CLHEP::neutron_mass_c2;
  theTarget = aNucleus.GetBiasedThermalNucleus(targetMass / CLHEP::neutron_mass_c2,
                                               neuVelo, theTrack.GetMaterial()->GetTemperature());
 
  theTarget.SetDefinition(G4IonTable::GetIonTable()->GetIon(theBaseZ, theBaseA, 0.0));
 
  // prepare the residual mass
  G4double residualMass = 0;
  G4int residualZ = theBaseZ +
    G4lrint((theProjectile->GetPDGCharge() - aDefinition->GetPDGCharge())/CLHEP::eplus);
  G4int residualA = theBaseA + theProjectile->GetBaryonNumber() - aDefinition->GetBaryonNumber();
  residualMass = G4NucleiProperties::GetNuclearMass(residualA, residualZ);
 
  // prepare energy in target rest frame
  G4ReactionProduct boosted;
  boosted.Lorentz(incidReactionProduct, theTarget);
  eKinetic = boosted.GetKineticEnergy();
 
  // select exit channel for composite FS class.
  G4int it = SelectExitChannel(eKinetic);
 
  // E. Mendoza (2018) -- to use JENDL/AN-2005
  if (theEnergyDistribution[it] == nullptr && theAngularDistribution[it] == nullptr
      && theEnergyAngData[it] == nullptr)
  {
    if (theEnergyDistribution[50] != nullptr || theAngularDistribution[50] != nullptr
        || theEnergyAngData[50] != nullptr)
    {
      it = 50;
    }
  }
 
  // set target and neutron in the relevant exit channel
  InitDistributionInitialState(incidReactionProduct, theTarget, it);
 
  //---------------------------------------------------------------------//
  // Hook for NRESP71MODEL
  if (fManager->GetUseNRESP71Model() && eKinetic < 20 * CLHEP::MeV) {
    if (theBaseZ == 6)  // If the reaction is with Carbon...
    {
      if (theProjectile == G4Neutron::Definition()) {
        if (use_nresp71_model(aDefinition, it, theTarget, boosted)) return;
      }
    }
  }
  //---------------------------------------------------------------------//
 
  G4ReactionProductVector* thePhotons = nullptr;
  G4ReactionProductVector* theParticles = nullptr;
  G4ReactionProduct aHadron;
  aHadron.SetDefinition(aDefinition);  // what if only cross-sections exist ==> Na 23 11 @@@@
  G4double availableEnergy = incidReactionProduct.GetKineticEnergy()
                             + incidReactionProduct.GetMass() - aHadron.GetMass()
                             + (targetMass - residualMass);
 
  if (availableEnergy < 0) {
    availableEnergy = 0;
  }
  G4int nothingWasKnownOnHadron = 0;
  G4double eGamm = 0;
  G4int iLevel = -1;
  // max gamma energy and index
  G4int imaxEx = theGammas.GetNumberOfLevels() - 1;
 
  // without photon has it = 0
  if (50 == it) {
    // Excitation level is not determined
    aHadron.SetKineticEnergy(availableEnergy * residualMass / (aHadron.GetMass() + residualMass));
 
    // TK add safty 100909
    G4double p2 =
      (aHadron.GetTotalEnergy() * aHadron.GetTotalEnergy() - aHadron.GetMass() * aHadron.GetMass());
    G4double p = (p2 > 0.0) ? std::sqrt(p2) : 0.0;
    aHadron.SetMomentum(p * incidReactionProduct.GetMomentum() /
            incidReactionProduct.GetTotalMomentum());
  }
  else {
    iLevel = imaxEx;
  }
 
  if (theAngularDistribution[it] != nullptr)  // MF4
  {
    if (theEnergyDistribution[it] != nullptr)  // MF5
    {
      //************************************************************
      /*
            aHadron.SetKineticEnergy(theEnergyDistribution[it]->Sample(eKinetic, dummy));
            G4double eSecN = aHadron.GetKineticEnergy();
      */
      //************************************************************
      // EMendoza --> maximum allowable energy should be taken into account.
      G4double dqi = 0.0;
      if (QI[it] < 0 || 849 < QI[it])
        dqi = QI[it];  // For backword compatibility QI introduced since G4NDL3.15
      G4double MaxEne = eKinetic + dqi;
      G4double eSecN = 0.;
 
      G4int icounter = 0;
      G4int icounter_max = 1024;
      G4int dummy = 0;
      do {
        ++icounter;
        if (icounter > icounter_max) {
          G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
                 << __FILE__ << "." << G4endl;
          break;
        }
        eSecN = theEnergyDistribution[it]->Sample(eKinetic, dummy);
      } while (eSecN > MaxEne);  // Loop checking, 11.05.2015, T. Koi
      aHadron.SetKineticEnergy(eSecN);
      //************************************************************
      eGamm = eKinetic - eSecN;
      for (iLevel = imaxEx; iLevel >= 0; --iLevel) {
        if (theGammas.GetLevelEnergy(iLevel) < eGamm) break;
      }
      if (iLevel < imaxEx && iLevel >= 0) {
        if (G4UniformRand() > 0.5) {
          ++iLevel;
        }
      }
    }
    else {
      G4double eExcitation = 0;
      for (iLevel = imaxEx; iLevel >= 0; --iLevel) {
        if (theGammas.GetLevelEnergy(iLevel) < eKinetic) break;
      }
 
      // Use QI value for calculating excitation energy of residual.
      G4bool useQI = false;
      G4double dqi = QI[it];
      if (dqi < 0 || 849 < dqi) useQI = true;  // Former libraries do not have values in this range
 
      if (useQI) {
        eExcitation = std::max(0., QI[0] - QI[it]);  // Bug fix 2333
 
        // Re-evaluate iLevel based on this eExcitation
        iLevel = 0;
        G4bool find = false;
        const G4double level_tolerance = 1.0 * CLHEP::keV;
 
        // VI: the first level is ground
        if (0 < imaxEx) {
          for (iLevel = 1; iLevel <= imaxEx; ++iLevel) {
            G4double elevel = theGammas.GetLevelEnergy(iLevel);
            if (std::abs(eExcitation - elevel) < level_tolerance) {
              find = true;
              break;
            }
            if (eExcitation < elevel) {
              find = true;
              iLevel = std::max(iLevel - 1, 0);
              break;
            }
          }
 
          // If proper level cannot be found, use the maximum level
          if (!find) iLevel = imaxEx;
        }
      }
 
      if (fManager->GetDEBUG() && eKinetic - eExcitation < 0) {
        throw G4HadronicException(
          __FILE__, __LINE__,
          "SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
      }
      if (eKinetic - eExcitation < 0) eExcitation = 0;
      if (iLevel != -1) aHadron.SetKineticEnergy(eKinetic - eExcitation);
    }
    theAngularDistribution[it]->SampleAndUpdate(aHadron);
 
    if (theFinalStatePhotons[it] == nullptr) {
      thePhotons = theGammas.GetDecayGammas(iLevel);
      eGamm -= theGammas.GetLevelEnergy(iLevel);
    }
  }
  else if (theEnergyAngData[it] != nullptr)  // MF6
  {
    theParticles = theEnergyAngData[it]->Sample(eKinetic);
 
    // Adjust A and Z in the case of miss much between selected data and target nucleus
    if (theParticles != nullptr) {
      G4int sumA = 0;
      G4int sumZ = 0;
      G4int maxA = 0;
      G4int jAtMaxA = 0;
      for (G4int j = 0; j != (G4int)theParticles->size(); ++j) {
    auto ptr = theParticles->at(j);
    G4int barnum = ptr->GetDefinition()->GetBaryonNumber();
        if (barnum > maxA) {
          maxA = barnum;
          jAtMaxA = j;
        }
        sumA += barnum;
        sumZ += G4lrint(ptr->GetDefinition()->GetPDGCharge()/CLHEP::eplus);
      }
      G4int dA = theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber() - sumA;
      G4int dZ = theBaseZ +
    G4lrint(hadProjectile->GetDefinition()->GetPDGCharge()/CLHEP::eplus) - sumZ;
      if (dA < 0 || dZ < 0) {
        G4int newA = theParticles->at(jAtMaxA)->GetDefinition()->GetBaryonNumber() + dA;
        G4int newZ = 
      G4lrint(theParticles->at(jAtMaxA)->GetDefinition()->GetPDGCharge()/CLHEP::eplus) + dZ;
        G4ParticleDefinition* pd = ionTable->GetIon(newZ, newA);
        theParticles->at(jAtMaxA)->SetDefinition(pd);
      }
    }
  }
  else {
    // @@@ what to do, if we have photon data, but no info on the hadron itself
    nothingWasKnownOnHadron = 1;
  }
 
  if (theFinalStatePhotons[it] != nullptr) {
    // the photon distributions are in the Nucleus rest frame.
    // TK residual rest frame
    G4ReactionProduct boosted_tmp;
    boosted_tmp.Lorentz(incidReactionProduct, theTarget);
    G4double anEnergy = boosted_tmp.GetKineticEnergy();
    thePhotons = theFinalStatePhotons[it]->GetPhotons(anEnergy);
    G4double aBaseEnergy = theFinalStatePhotons[it]->GetLevelEnergy();
    G4double testEnergy = 0;
    if (thePhotons != nullptr && !thePhotons->empty()) {
      aBaseEnergy -= (*thePhotons)[0]->GetTotalEnergy();
    }
    if (theFinalStatePhotons[it]->NeedsCascade()) {
      while (aBaseEnergy > 0.01 * CLHEP::keV)  // Loop checking, 11.05.2015, T. Koi
      {
        // cascade down the levels
        G4bool foundMatchingLevel = false;
        G4int closest = 2;
        G4double deltaEold = -1;
        for (G4int j = 1; j < it; ++j) {
          if (theFinalStatePhotons[j] != nullptr) {
            testEnergy = theFinalStatePhotons[j]->GetLevelEnergy();
          }
          else {
            testEnergy = 0;
          }
          G4double deltaE = std::abs(testEnergy - aBaseEnergy);
          if (deltaE < 0.1 * CLHEP::keV) {
            G4ReactionProductVector* theNext = theFinalStatePhotons[j]->GetPhotons(anEnergy);
            if (thePhotons != nullptr) thePhotons->push_back(theNext->operator[](0));
            aBaseEnergy = testEnergy - theNext->operator[](0)->GetTotalEnergy();
            delete theNext;
            foundMatchingLevel = true;
            break;  // ===>
          }
          if (theFinalStatePhotons[j] != nullptr && (deltaE < deltaEold || deltaEold < 0.)) {
            closest = j;
            deltaEold = deltaE;
          }
        }  // <=== the break goes here.
        if (!foundMatchingLevel) {
          G4ReactionProductVector* theNext = theFinalStatePhotons[closest]->GetPhotons(anEnergy);
          if (thePhotons != nullptr) thePhotons->push_back(theNext->operator[](0));
          aBaseEnergy = aBaseEnergy - theNext->operator[](0)->GetTotalEnergy();
          delete theNext;
        }
      }
    }
  }
 
  if (thePhotons != nullptr) {
    for (auto const & p : *thePhotons) {
      // back to lab
      p->Lorentz(*p, -1. * theTarget);
    }
  }
  if (nothingWasKnownOnHadron != 0) {
    // In this case, hadron should be isotropic in CM
    // Next 12 lines are Emilio's replacement
    // G4double QM=(incidReactionProduct.GetMass()+targetMass)-(aHadron.GetMass()+residualMass);
    // G4double eExcitation = QM-QI[it];
    // G4double eExcitation = QI[0] - QI[it];  // Fix of bug #1838
    // if(eExcitation<20*CLHEP::keV){eExcitation=0;}
 
    G4double eExcitation = std::max(0., QI[0] - QI[it]);  // Fix of bug #2333
 
    two_body_reaction(&incidReactionProduct, &theTarget, &aHadron, eExcitation);
    if (thePhotons == nullptr && eExcitation > 0) {
      for (iLevel = imaxEx; iLevel >= 0; --iLevel) {
        if (theGammas.GetLevelEnergy(iLevel) < eExcitation + 5 * keV) break;  // 5 keV tolerance
      }
      thePhotons = theGammas.GetDecayGammas(iLevel);
    }
  }
 
  // fill the result
  // Beware - the recoil is not necessarily in the particles...
  // Can be calculated from momentum conservation?
  // The idea is that the particles ar emitted forst, and the gammas only once the
  // recoil is on the residual; assumption is that gammas do not contribute to
  // the recoil.
  // This needs more design @@@
 
  G4bool needsSeparateRecoil = false;
  G4int totalBaryonNumber = 0;
  G4int totalCharge = 0;
  G4ThreeVector totalMomentum(0);
  if (theParticles != nullptr) {
    const G4ParticleDefinition* aDef;
    for (std::size_t ii0 = 0; ii0 < theParticles->size(); ++ii0) {
      aDef = (*theParticles)[ii0]->GetDefinition();
      totalBaryonNumber += aDef->GetBaryonNumber();
      totalCharge += G4lrint(aDef->GetPDGCharge()/CLHEP::eplus);
      totalMomentum += (*theParticles)[ii0]->GetMomentum();
    }
    if (totalBaryonNumber
        != theBaseA +  hadProjectile->GetDefinition()->GetBaryonNumber())
    {
      needsSeparateRecoil = true;
      residualA = theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber()
    - totalBaryonNumber;
      residualZ = theBaseZ +
    G4lrint((hadProjectile->GetDefinition()->GetPDGCharge() - totalCharge)/CLHEP::eplus);
    }
  }
 
  std::size_t nPhotons = 0;
  if (thePhotons != nullptr) {
    nPhotons = thePhotons->size();
  }
 
  G4DynamicParticle* theSec;
 
  if (theParticles == nullptr) {
    theSec = new G4DynamicParticle;
    theSec->SetDefinition(aHadron.GetDefinition());
    theSec->SetMomentum(aHadron.GetMomentum());
    theResult.Get()->AddSecondary(theSec, secID);
#ifdef G4VERBOSE
    if (fManager->GetDEBUG())
      G4cout << " G4ParticleHPInelasticCompFS::BaseApply  add secondary1 "
             << theSec->GetParticleDefinition()->GetParticleName()
             << " E= " << theSec->GetKineticEnergy() << " NSECO "
             << theResult.Get()->GetNumberOfSecondaries() << G4endl;
#endif
 
    aHadron.Lorentz(aHadron, theTarget);
    G4ReactionProduct theResidual;
    theResidual.SetDefinition(ionTable->GetIon(residualZ, residualA, 0));
    theResidual.SetKineticEnergy(aHadron.GetKineticEnergy() * aHadron.GetMass()
                                 / theResidual.GetMass());
 
    // 080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #6
    // theResidual.SetMomentum(-1.*aHadron.GetMomentum());
    G4ThreeVector incidentNeutronMomentum = incidReactionProduct.GetMomentum();
    theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
 
    theResidual.Lorentz(theResidual, -1. * theTarget);
    G4ThreeVector totalPhotonMomentum(0, 0, 0);
    if (thePhotons != nullptr) {
      for (std::size_t i = 0; i < nPhotons; ++i) {
        totalPhotonMomentum += (*thePhotons)[i]->GetMomentum();
      }
    }
    theSec = new G4DynamicParticle;
    theSec->SetDefinition(theResidual.GetDefinition());
    theSec->SetMomentum(theResidual.GetMomentum() - totalPhotonMomentum);
    theResult.Get()->AddSecondary(theSec, secID);
#ifdef G4VERBOSE
    if (fManager->GetDEBUG())
      G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary2 "
             << theSec->GetParticleDefinition()->GetParticleName()
             << " E= " << theSec->GetKineticEnergy() << " NSECO "
             << theResult.Get()->GetNumberOfSecondaries() << G4endl;
#endif
  }
  else {
    for (std::size_t i0 = 0; i0 < theParticles->size(); ++i0) {
      theSec = new G4DynamicParticle;
      theSec->SetDefinition((*theParticles)[i0]->GetDefinition());
      theSec->SetMomentum((*theParticles)[i0]->GetMomentum());
      theResult.Get()->AddSecondary(theSec, secID);
#ifdef G4VERBOSE
      if (fManager->GetDEBUG())
        G4cout << " G4ParticleHPInelasticCompFS::BaseApply add secondary3 "
               << theSec->GetParticleDefinition()->GetParticleName()
               << " E= " << theSec->GetKineticEnergy() << " NSECO "
               << theResult.Get()->GetNumberOfSecondaries() << G4endl;
#endif
      delete (*theParticles)[i0];
    }
    delete theParticles;
    if (needsSeparateRecoil && residualZ != 0) {
      G4ReactionProduct theResidual;
      theResidual.SetDefinition(ionTable->GetIon(residualZ, residualA, 0));
      G4double resiualKineticEnergy = theResidual.GetMass() * theResidual.GetMass();
      resiualKineticEnergy += totalMomentum * totalMomentum;
      resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.GetMass();
      theResidual.SetKineticEnergy(resiualKineticEnergy);
 
      // 080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4
      // theResidual.SetMomentum(-1.*totalMomentum);
      // G4ThreeVector incidentNeutronMomentum = incidReactionProduct.GetMomentum();
      // theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
      // 080717 TK Comment still do NOT include photon's mometum which produce by thePhotons
      theResidual.SetMomentum(incidReactionProduct.GetMomentum() + theTarget.GetMomentum()
                              - totalMomentum);
 
      theSec = new G4DynamicParticle;
      theSec->SetDefinition(theResidual.GetDefinition());
      theSec->SetMomentum(theResidual.GetMomentum());
      theResult.Get()->AddSecondary(theSec, secID);
#ifdef G4VERBOSE
      if (fManager->GetDEBUG())
        G4cout << " G4ParticleHPInelasticCompFS::BaseApply add secondary4 "
               << theSec->GetParticleDefinition()->GetParticleName()
               << " E= " << theSec->GetKineticEnergy() << " NSECO "
               << theResult.Get()->GetNumberOfSecondaries() << G4endl;
#endif
    }
  }
  if (thePhotons != nullptr) {
    for (std::size_t i = 0; i < nPhotons; ++i) {
      theSec = new G4DynamicParticle;
      // Bug reported Chao Zhang ([email protected]), Dongming Mei([email protected]) Feb. 25,
      // 2009 theSec->SetDefinition(G4Gamma::Gamma());
      theSec->SetDefinition((*thePhotons)[i]->GetDefinition());
      // But never cause real effect at least with G4NDL3.13 TK
      theSec->SetMomentum((*thePhotons)[i]->GetMomentum());
      theResult.Get()->AddSecondary(theSec, secID);
#ifdef G4VERBOSE
      if (fManager->GetDEBUG())
        G4cout << " G4ParticleHPInelasticCompFS::BaseApply add secondary5 "
               << theSec->GetParticleDefinition()->GetParticleName()
               << " E= " << theSec->GetKineticEnergy() << " NSECO "
               << theResult.Get()->GetNumberOfSecondaries() << G4endl;
#endif
 
      delete thePhotons->operator[](i);
    }
    // some garbage collection
    delete thePhotons;
  }
 
  G4ParticleDefinition* targ_pd = ionTable->GetIon(theBaseZ, theBaseA, 0.0);
  G4LorentzVector targ_4p_lab(
    theTarget.GetMomentum(),
    std::sqrt(targ_pd->GetPDGMass() * targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2()));
  G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
  G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
  adjust_final_state(init_4p_lab);
 
  // clean up the primary neutron
  theResult.Get()->SetStatusChange(stopAndKill);
}

Referenced by G4ParticleHPAInelasticFS::ApplyYourself(), G4ParticleHPDInelasticFS::ApplyYourself(), G4ParticleHPHe3InelasticFS::ApplyYourself(), G4ParticleHPNInelasticFS::ApplyYourself(), G4ParticleHPPInelasticFS::ApplyYourself(), and G4ParticleHPTInelasticFS::ApplyYourself().

◆ GetXsec() [1/2]

G4ParticleHPVector * G4ParticleHPInelasticCompFS::GetXsec ( )

inlineoverridevirtual

Reimplemented from G4ParticleHPFinalState.

Definition at line 71 of file G4ParticleHPInelasticCompFS.hh.

71{ return theXsection[50]; }

Referenced by SelectExitChannel().

◆ GetXsec() [2/2]

G4double G4ParticleHPInelasticCompFS::GetXsec ( G4double anEnergy )

inlineoverridevirtual

Reimplemented from G4ParticleHPFinalState.

Definition at line 66 of file G4ParticleHPInelasticCompFS.hh.

  {
    return std::max(0., theXsection[50]->GetY(anEnergy));
  }

◆ Init()

void G4ParticleHPInelasticCompFS::Init	(	G4double	A,
		G4double	Z,
		G4int	M,
		G4String &	dirName,
		G4String &	aSFType,
		G4ParticleDefinition *	)

overridevirtual

Implements G4ParticleHPFinalState.

Reimplemented in G4ParticleHPNInelasticFS, G4ParticleHPPInelasticFS, and G4ParticleHPTInelasticFS.

Definition at line 129 of file G4ParticleHPInelasticCompFS.cc.

{
  gammaPath = fManager->GetNeutronHPPath() + "/Inelastic/Gammas/";
  G4String tString = dirName;
  SetA_Z(A, Z, M);
  G4bool dbool;
  G4ParticleHPDataUsed aFile =
    theNames.GetName(theBaseA, theBaseZ, M, tString, aFSType, dbool);
  SetAZMs(aFile);
  G4String filename = aFile.GetName();
#ifdef G4VERBOSE
  if (fManager->GetDEBUG())
    G4cout << " G4ParticleHPInelasticCompFS::Init FILE " << filename << G4endl;
#endif
 
  SetAZMs(A, Z, M, aFile);
 
  if (!dbool || (theBaseZ <= 2 && (theNDLDataZ != theBaseZ || theNDLDataA != theBaseA)))
  {
#ifdef G4VERBOSE
    if (fManager->GetDEBUG())
      G4cout << "Skipped = " << filename << " " << A << " " << Z << G4endl;
#endif
    hasAnyData = false;
    hasFSData = false;
    hasXsec = false;
    return;
  }
  std::istringstream theData(std::ios::in);
  fManager->GetDataStream(filename, theData);
  if (!theData)  //"!" is a operator of ios
  {
    hasAnyData = false;
    hasFSData = false;
    hasXsec = false;
    return;
  }
  // here we go
  G4int infoType, dataType, dummy;
  G4int sfType, it;
  hasFSData = false;
  while (theData >> infoType)  // Loop checking, 11.05.2015, T. Koi
  {
    hasFSData = true;
    theData >> dataType;
    theData >> sfType >> dummy;
    it = 50;
    if (sfType >= 600 || (sfType < 100 && sfType >= 50))
      it = sfType % 50;
    if (dataType == 3)
    {
      G4double dqi;
      G4int ilr;
      theData >> dqi >> ilr;
 
      QI[it] = dqi * CLHEP::eV;
      LR[it] = ilr;
      theXsection[it] = new G4ParticleHPVector;
      G4int total;
      theData >> total;
      theXsection[it]->Init(theData, total, CLHEP::eV);
    }
    else if (dataType == 4) {
      theAngularDistribution[it] = new G4ParticleHPAngular;
      theAngularDistribution[it]->Init(theData);
    }
    else if (dataType == 5) {
      theEnergyDistribution[it] = new G4ParticleHPEnergyDistribution;
      theEnergyDistribution[it]->Init(theData);
    }
    else if (dataType == 6) {
      theEnergyAngData[it] = new G4ParticleHPEnAngCorrelation(theProjectile);
      //      G4cout << this << " CompFS theEnergyAngData " << it << theEnergyAngData[it] << G4endl;
      theEnergyAngData[it]->Init(theData);
    }
    else if (dataType == 12) {
      theFinalStatePhotons[it] = new G4ParticleHPPhotonDist;
      theFinalStatePhotons[it]->InitMean(theData);
    }
    else if (dataType == 13) {
      theFinalStatePhotons[it] = new G4ParticleHPPhotonDist;
      theFinalStatePhotons[it]->InitPartials(theData, theXsection[50]);
    }
    else if (dataType == 14) {
      theFinalStatePhotons[it]->InitAngular(theData);
    }
    else if (dataType == 15) {
      theFinalStatePhotons[it]->InitEnergies(theData);
    }
    else {
      G4ExceptionDescription ed;
      ed << "Z=" << theBaseZ << " A=" << theBaseA << " dataType=" << dataType
     << " projectile: " << theProjectile->GetParticleName();
      G4Exception("G4ParticleHPInelasticCompFS::Init", "hadr01", JustWarning,
          ed, "Data-type unknown");
    }
  }
}

Referenced by G4ParticleHPAInelasticFS::Init(), G4ParticleHPDInelasticFS::Init(), G4ParticleHPHe3InelasticFS::Init(), G4ParticleHPNInelasticFS::Init(), G4ParticleHPPInelasticFS::Init(), and G4ParticleHPTInelasticFS::Init().

◆ InitDistributionInitialState()

void G4ParticleHPInelasticCompFS::InitDistributionInitialState	(	G4ReactionProduct &	anIncidentPart,
		G4ReactionProduct &	aTarget,
		G4int	it )

Definition at line 100 of file G4ParticleHPInelasticCompFS.cc.

{
  if (theAngularDistribution[it] != nullptr) {
    theAngularDistribution[it]->SetTarget(aTarget);
    theAngularDistribution[it]->SetProjectileRP(inPart);
  }
 
  if (theEnergyAngData[it] != nullptr) {
    theEnergyAngData[it]->SetTarget(aTarget);
    theEnergyAngData[it]->SetProjectileRP(inPart);
  }
}

Referenced by CompositeApply().

◆ InitGammas()

void G4ParticleHPInelasticCompFS::InitGammas	(	G4double	AR,
		G4double	ZR )

Definition at line 114 of file G4ParticleHPInelasticCompFS.cc.

{
  G4int Z = G4lrint(ZR);
  G4int A = G4lrint(AR);
  std::ostringstream ost;
  ost << gammaPath << "z" << Z << ".a" << A;
  G4String aName = ost.str();
  std::ifstream from(aName, std::ios::in);
 
  if (!from) return;  // no data found for this isotope
  std::ifstream theGammaData(aName, std::ios::in);
 
  theGammas.Init(theGammaData);
}

Referenced by G4ParticleHPAInelasticFS::Init(), G4ParticleHPDInelasticFS::Init(), G4ParticleHPHe3InelasticFS::Init(), G4ParticleHPNInelasticFS::Init(), G4ParticleHPPInelasticFS::Init(), and G4ParticleHPTInelasticFS::Init().

◆ New()

G4ParticleHPFinalState * G4ParticleHPInelasticCompFS::New ( )

overridepure virtual

Implements G4ParticleHPFinalState.

Implemented in G4ParticleHPAInelasticFS, G4ParticleHPDInelasticFS, G4ParticleHPHe3InelasticFS, G4ParticleHPNInelasticFS, G4ParticleHPPInelasticFS, and G4ParticleHPTInelasticFS.

◆ operator=()

G4ParticleHPInelasticCompFS & G4ParticleHPInelasticCompFS::operator= ( const G4ParticleHPInelasticCompFS & right )

delete

◆ SelectExitChannel()

G4int G4ParticleHPInelasticCompFS::SelectExitChannel ( G4double eKinetic )

Definition at line 229 of file G4ParticleHPInelasticCompFS.cc.

{
  G4double running[50];
  running[0] = 0;
  G4int i;
  for (i = 0; i < 50; ++i) {
    if (i != 0) running[i] = running[i - 1];
    if (theXsection[i] != nullptr) {
      running[i] += std::max(0., theXsection[i]->GetXsec(eKinetic));
    }
  }
  G4double random = G4UniformRand();
  G4double sum = running[49];
  G4int it = 50;
  if (0 != sum) {
    G4int i0;
    for (i0 = 0; i0 < 50; ++i0) {
      it = i0;
      if (random < running[i0] / sum) break;
    }
  }
  return it;
}