lars/geant4
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Import Geant4 10.6.2 source tree

Browse Source
This commit is contained in:
Gabriele Cosmo
2020-05-29 14:54:29 +02:00
parent 8c87fe78c4
commit c02c370437
448 changed files with 23779 additions and 31516 deletions
Download Patch File Download Diff File Expand all files Collapse all files
source/processes/hadronic/models/lend/src/G4LENDGammaCrossSection.cc
+99
View File
@@ -0,0 +1,99 @@
//
// ********************************************************************
// * 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. *
// ********************************************************************
//
////////////////////////////////////////////////////////////////////////////////
// //
// File: G4LENDGammaCrossSection.cc //
// Date: 30 March 2020 //
// Author: Dennis H. Wright //
// //
// Description: cross sections for inelastic scattering of gammas from //
// nuclei including gamma-induced fission. This cross section //
// is very similar to G4LENDCombinedCrossSection except that //
// does not sample elastic or capture reactions since there are //
// no such data for gammas in GND. //
// //
////////////////////////////////////////////////////////////////////////////////
#include "G4LENDGammaCrossSection.hh"
#include "G4LENDInelasticCrossSection.hh"
#include "G4LENDFissionCrossSection.hh"
#include "Randomize.hh"
G4LENDGammaCrossSection::G4LENDGammaCrossSection(G4ParticleDefinition* pd)
:G4LENDCrossSection("LENDGammaCrossSection")
{
proj = pd;
inelasticXS = new G4LENDInelasticCrossSection(pd);
fissionXS = new G4LENDFissionCrossSection(pd);
}
void G4LENDGammaCrossSection::BuildPhysicsTable(const G4ParticleDefinition& pd)
{
inelasticXS->BuildPhysicsTable(pd);
fissionXS->BuildPhysicsTable(pd);
create_used_target_map();
}
G4double
G4LENDGammaCrossSection::GetIsoCrossSection(const G4DynamicParticle* dp,
G4int iZ, G4int iA,
const G4Isotope* isotope,
const G4Element*,
const G4Material* material)
{
G4double XS = 0.0;
XS += inelasticXS->GetIsoCrossSection(dp, iZ, iA, isotope, NULL, material);
XS += fissionXS->GetIsoCrossSection(dp, iZ, iA, isotope, NULL, material);
//G4cout << "G4LENDGammaCrossSection::GetIsoCrossSection "
// << XS/CLHEP::barn << " [barn]" << G4endl;
return XS;
}
G4int G4LENDGammaCrossSection::SelectChannel(const G4DynamicParticle* dp,
G4int iZ, G4int iA,
const G4Isotope* isotope,
const G4Element*,
const G4Material* material)
{
G4int ichannel = -1;
G4double XSs[2];
XSs[0] = inelasticXS->GetIsoCrossSection(dp, iZ, iA, isotope, nullptr, material);
XSs[1] = XSs[0] + fissionXS->GetIsoCrossSection(dp, iZ, iA, isotope, nullptr, material);
G4double total = XSs[1];
G4double random = G4UniformRand();
for (G4int i = 0; i < 2; i++) {
if (random*total <= XSs[i]) {
ichannel = i;
break;
}
}
return ichannel;
}
source/processes/hadronic/models/lend/src/G4LENDGammaModel.cc
+87
View File
@@ -0,0 +1,87 @@
//
// ********************************************************************
// * 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. *
// ********************************************************************
//
////////////////////////////////////////////////////////////////////////////////
// //
// File: G4LENDGammaModel.cc //
// Date: 30 March 2020 //
// Author: Dennis H. Wright //
// //
// Description: model for inelastic scattering of gammas from nuclei //
// including gamma-induced fission. This model is very similar //
// to G4LENDCombinedModel except that it does not sample //
// elastic or capture reactions since there are no such data //
// for gammas in GND. //
// //
////////////////////////////////////////////////////////////////////////////////
#include "G4LENDGammaModel.hh"
#include "G4LENDGammaCrossSection.hh"
#include "G4LENDInelastic.hh"
#include "G4LENDFission.hh"
#include "G4DynamicParticle.hh"
G4LENDGammaModel::G4LENDGammaModel(G4ParticleDefinition* pd)
:G4LENDModel("LENDGammaModel") {
proj = pd;
crossSection = new G4LENDGammaCrossSection(pd);
inelastic = new G4LENDInelastic(pd);
fission = new G4LENDFission(pd);
channels[0] = inelastic;
channels[1] = fission;
}
void G4LENDGammaModel::BuildPhysicsTable(const G4ParticleDefinition& projectile) {
crossSection->BuildPhysicsTable(projectile);
create_used_target_map();
}
G4bool G4LENDGammaModel::HasData(const G4DynamicParticle*, G4int iZ, G4int iA,
G4int iM, const G4Isotope*, const G4Element*,
const G4Material*)
{
G4bool result = false;
if (get_target_from_map(lend_manager->GetNucleusEncoding(iZ, iA, iM) ) != nullptr) result = true;
return result;
}
G4HadFinalState* G4LENDGammaModel::ApplyYourself(const G4HadProjectile& aTrack,
G4Nucleus& aTarg)
{
G4LENDModel* channel = nullptr;
G4int iZ = aTarg.GetZ_asInt();
G4int iA = aTarg.GetA_asInt();
//To pass kinetic energy, need to generate dynamic particle
G4DynamicParticle* dp = new G4DynamicParticle(proj, G4ThreeVector(0.,0.,1.),
aTrack.GetKineticEnergy() );
G4int ichannel = crossSection->SelectChannel(dp, iZ, iA, aTarg.GetIsotope(),
nullptr, aTrack.GetMaterial() );
delete dp;
channel = channels[ichannel];
return channel->ApplyYourself(aTrack, aTarg);
}
source/processes/hadronic/models/lend/src/G4LENDInelastic.cc
+218 -140
View File
@@ -23,166 +23,244 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
#include "G4LENDInelastic.hh"
////////////////////////////////////////////////////////////////////////////////
// //
// File: G4LENDInelastic.cc //
// Date: 24 March 2020 //
// Author: Dennis Wright //
// //
// Description: model for inelastic scattering of neutrons, light ions and //
// gammas at energies of order 20 MeV and lower. //
// This model uses GIDI particle data which are stored mostly //
// in spectrum mode. In this mode, spectra are reproduced //
// for each possible particle type which can result from a //
// given interaction. Unlike Geant4, this is done without //
// consideration of event-by-event conservation rules. //
// Indeed, forcing such conservation on GIDI output products //
// introduces correlations and distortions in the resulting //
// spectra which are not present in the data. //
// //
// In order to use GIDI data within the Geant4 framework, a //
// minimal event-by-event baryon number conservation is //
// enforced which allows deviations of up to 1 GeV without //
// giving warnings. Neither charge, nor energy, nor momentum //
// conservation is enforced. Under this scheme, light //
// fragment (n, p, d, t, alpha) spectra are well reproduced //
// after a large number of events. Charge and energy //
// conservation also approach their event-by-event values in //
// this limit. The mass, charge and energy distributions of //
// large fragments, however, are not expected to reproduce the //
// data very well. This is a result of forcing the crude //
// baryon number conservation and ensuring that the light //
// fragment spectra are correct. //
// //
////////////////////////////////////////////////////////////////////////////////
#include "G4LENDInelastic.hh"
#include "G4SystemOfUnits.hh"
#include "G4Nucleus.hh"
#include "G4IonTable.hh"
#include <algorithm>
#include <random>
G4HadFinalState * G4LENDInelastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aTarg )
G4HadFinalState* G4LENDInelastic::ApplyYourself(const G4HadProjectile& aTrack,
G4Nucleus& aTarg)
{
//return preco->ApplyYourself( aTrack, aTarg );
G4ThreeVector projMom = aTrack.Get4Momentum().vect();
G4double temp = aTrack.GetMaterial()->GetTemperature();
G4ThreeVector proj_p = aTrack.Get4Momentum().vect();
G4int iZ = aTarg.GetZ_asInt();
G4int iA = aTarg.GetA_asInt();
G4int iM = 0;
if (aTarg.GetIsotope() != nullptr) iM = aTarg.GetIsotope()->Getm();
G4double temp = aTrack.GetMaterial()->GetTemperature();
G4double ke = aTrack.GetKineticEnergy();
//G4int iZ = int ( aTarg.GetZ() );
//G4int iA = int ( aTarg.GetN() );
//migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this)
G4int iZ = aTarg.GetZ_asInt();
G4int iA = aTarg.GetA_asInt();
//G4int iM = aTarg.GetM_asInt();
G4int iM = 0;
if ( aTarg.GetIsotope() != NULL ) {
iM = aTarg.GetIsotope()->Getm();
}
//G4cout << "target: Z = " << iZ << " N = " << iA << G4endl;
G4HadFinalState* theResult = &theParticleChange;
theResult->Clear();
G4double ke = aTrack.GetKineticEnergy();
//G4cout << "projectile: KE = " << ke/MeV << " [MeV]" << G4endl;
G4GIDI_target* aGIDITarget =
get_target_from_map(lend_manager->GetNucleusEncoding(iZ, iA, iM) );
if (aGIDITarget == nullptr) {
// G4cout << " No target found " << G4endl;
theParticleChange.Clear();
theParticleChange.SetStatusChange(isAlive);
theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
return &theParticleChange;
}
// return returnUnchanged(aTrack, theResult);
G4HadFinalState* theResult = &theParticleChange;
theResult->Clear();
// Get GIDI final state products for givern target and projectile
G4int loop(0);
G4int loopMax = 1000;
std::vector<G4GIDI_Product>* products;
do {
products = aGIDITarget->getOthersFinalState(ke*MeV, temp, MyRNG, NULL);
loop++;
} while (products == nullptr && loop < loopMax);
G4GIDI_target* aTarget = get_target_from_map( lend_manager->GetNucleusEncoding( iZ , iA , iM ) );
if ( aTarget == NULL ) return returnUnchanged( aTrack , theResult );
// G4LENDInelastic accepts all light fragments and gammas from GIDI (A < 5)
// and removes any heavy fragments which cause large baryon number violation.
// Charge and energy non-conservation still occur, but over a large number
// of events, this improves on average.
std::vector<G4GIDI_Product>* products;
for ( G4int i = 0 ; i != 1024 ; i++ ) {
products = aTarget->getOthersFinalState( ke*MeV, temp, MyRNG, NULL );
if ( products != NULL ) break;
}
//return preco->ApplyYourself( aTrack, aTarg );
if (loop > loopMax - 1) {
// G4cout << " too many loops, return intial state " << G4endl;
G4int iTotZ = iZ + aTrack.GetDefinition()->GetAtomicNumber();
G4int iTotA = iA + aTrack.GetDefinition()->GetAtomicMass();
theParticleChange.Clear();
theParticleChange.SetStatusChange(isAlive);
theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
return &theParticleChange;
if ( products != NULL )
{
//G4cout << "Using LENDModel" << G4endl;
// if (aTrack.GetDefinition() == G4Proton::Proton() ||
// aTrack.GetDefinition() == G4Neutron::Neutron() ) {
// theResult = preco->ApplyYourself(aTrack, aTarg);
// } else {
// theResult = returnUnchanged(aTrack, theResult);
// }
G4ThreeVector psum(0);
G4bool needResidual = true;
int totN = 0;
int totZ = 0;
for ( G4int j = 0; j < int( products->size() ); j++ )
{
} else {
G4int iTotZ = iZ + aTrack.GetDefinition()->GetAtomicNumber();
G4int iTotA = iA + aTrack.GetDefinition()->GetAtomicMass();
G4int jZ = (*products)[j].Z;
G4int jA = (*products)[j].A;
G4int jm = (*products)[j].m;
//TK
//We need coordination LEND *products)[j].m and G4IonTable(Z,A,m)
//Excitation energy of isomer level is might (probably) different each other.
//
//G4cout << "ZA = " << 1000 * (*products)[j].Z + (*products)[j].A << " EK = "
// << (*products)[j].kineticEnergy
// << " px " << (*products)[j].px
// << " py " << (*products)[j].py
// << " pz " << (*products)[j].pz
// << G4endl;
iTotZ -= jZ;
iTotA -= jA;
G4DynamicParticle* theSec = new G4DynamicParticle;
if ( jA == 1 && jZ == 1 )
{
theSec->SetDefinition( G4Proton::Proton() );
totN += 1;
totZ += 1;
}
else if ( jA == 1 && jZ == 0 )
{
theSec->SetDefinition( G4Neutron::Neutron() );
totN += 1;
}
else if ( jZ > 0 )
{
if ( jA != 0 )
{
theSec->SetDefinition( G4IonTable::GetIonTable()->GetIon( jZ , jA , jm ) );
totN += jA;
totZ += jZ;
}
else
{
theSec->SetDefinition( G4IonTable::GetIonTable()->GetIon( jZ , iA+aTrack.GetDefinition()->GetAtomicMass()-totN , jm ) );
iTotZ -= jZ;
iTotA -= iA+aTrack.GetDefinition()->GetAtomicMass()-totN;
needResidual=false;
}
}
else
{
theSec->SetDefinition( G4Gamma::Gamma() );
}
G4ThreeVector p( (*products)[j].px*MeV , (*products)[j].py*MeV , (*products)[j].pz*MeV );
psum += p;
if ( p.mag() == 0 ) p = proj_p - psum;
theSec->SetMomentum( p );
theResult->AddSecondary( theSec );
}
if ( !( iTotZ == 0 && iTotA == 0 ) ) {
if ( iTotZ >= 0 && iTotA > 0 ) {
if ( needResidual ) {
G4DynamicParticle* residual = new G4DynamicParticle;
if ( iTotZ > 0 ) {
residual->SetDefinition( G4IonTable::GetIonTable()->GetIon( iTotZ , iTotA ) );
} else if ( iTotA == 1 ) {
residual->SetDefinition( G4Neutron::Neutron() );
} else {
//G4cout << "Charge or Baryon Number Error #3 iTotZ = " << iTotZ << ", iTotA = " << iTotA << G4endl;
;
}
residual->SetMomentum( proj_p - psum );
theResult->AddSecondary( residual );
} else {
//G4cout << "Charge or Baryon Number Error #1 iTotZ = " << iTotZ << ", iTotA = " << iTotA << G4endl;
;
}
} else {
if ( needResidual ) {
//G4cout << "Charge or Baryon Number Error #2 iTotZ = " << iTotZ << ", iTotA = " << iTotA << G4endl;
;
}
}
}
}
else {
//G4cout << "Using PreCompoundModel" << G4endl;
if ( aTrack.GetDefinition() == G4Proton::Proton() ||
aTrack.GetDefinition() == G4Neutron::Neutron() ) {
theResult = preco->ApplyYourself( aTrack, aTarg );
// Loop over GIDI products and separate light from heavy fragments
G4int GZtot(0);
G4int GAtot(0);
G4int productA(0);
G4int productZ(0);
std::vector<G4int> lightProductIndex;
std::vector<G4int> heavyProductIndex;
for (G4int i = 0; i < int( products->size() ); i++ ) {
productA = (*products)[i].A;
if (productA < 5) {
lightProductIndex.push_back(i);
GZtot += (*products)[i].Z;
GAtot += productA;
} else {
return theResult;
heavyProductIndex.push_back(i);
}
}
delete products;
}
theResult->SetStatusChange( stopAndKill );
// Randomize order of heavies to correct somewhat for sampling bias
// std::random_shuffle(heavyProductIndex.begin(), heavyProductIndex.end() );
// std::cout << " Heavy product index before shuffle : " ;
// for (G4int i = 0; i < int(heavyProductIndex.size() ); i++) std::cout << heavyProductIndex[i] << ", " ;
// std::cout << std::endl;
return theResult;
auto rng = std::default_random_engine {};
std::shuffle(heavyProductIndex.begin(), heavyProductIndex.end(), rng);
// std::cout << " Heavy product index after shuffle : " ;
// for (G4int i = 0; i < int(heavyProductIndex.size() ); i++) std::cout << heavyProductIndex[i] << ", " ;
// std::cout << std::endl;
std::vector<G4int> savedHeavyIndex;
G4int itest(0);
for (G4int i = 0; i < int(heavyProductIndex.size() ); i++) {
itest = heavyProductIndex[i];
productA = (*products)[itest].A;
productZ = (*products)[itest].Z;
if ((GAtot + productA <= iTotA) && (GZtot + productZ <= iTotZ) ) {
savedHeavyIndex.push_back(itest);
GZtot += productZ;
GAtot += productA;
}
}
/*
G4cout << " saved light products = ";
for (G4int k = 0; k < int(lightProductIndex.size() ); k++ ) {
itest = lightProductIndex[k];
G4cout << "(" << (*products)[itest].Z << ", " << (*products)[itest].A << "), ";
}
G4cout << G4endl;
G4cout << " saved heavy products = ";
for (G4int k = 0; k < int(savedHeavyIndex.size() ); k++ ) {
itest = savedHeavyIndex[k];
G4cout << "(" << (*products)[itest].Z << ", " << (*products)[itest].A << "), ";
}
G4cout << G4endl;
*/
// Now convert saved products to Geant4 particles
// Note that, at least for heavy fragments, GIDI masses and Geant4 masses
// have slightly different values.
G4DynamicParticle* theSec = nullptr;
G4ThreeVector Psum;
for (G4int i = 0; i < int(lightProductIndex.size() ); i++) {
itest = lightProductIndex[i];
productZ = (*products)[itest].Z;
productA = (*products)[itest].A;
theSec = new G4DynamicParticle();
if (productA == 1 && productZ == 0) {
theSec->SetDefinition(G4Neutron::Neutron() );
} else if (productA == 1 && productZ == 1) {
theSec->SetDefinition(G4Proton::Proton() );
} else if (productA == 2 && productZ == 1) {
theSec->SetDefinition(G4Deuteron::Deuteron() );
} else if (productA == 3 && productZ == 1) {
theSec->SetDefinition(G4Triton::Triton() );
} else if (productA == 4 && productZ == 2) {
theSec->SetDefinition(G4Alpha::Alpha() );
} else {
theSec->SetDefinition(G4Gamma::Gamma() );
}
G4ThreeVector momentum((*products)[itest].px*MeV,
(*products)[itest].py*MeV,
(*products)[itest].pz*MeV );
Psum += momentum;
theSec->SetMomentum(momentum);
// theResult->AddSecondary(theSec);
theParticleChange.AddSecondary(theSec);
}
G4int productM(0);
for (G4int i = 0; i < int(savedHeavyIndex.size() ); i++) {
itest = savedHeavyIndex[i];
productZ = (*products)[itest].Z;
productA = (*products)[itest].A;
productM = (*products)[itest].m;
theSec = new G4DynamicParticle();
theSec->SetDefinition(G4IonTable::GetIonTable()->GetIon(productZ,
productA,
productM) );
G4ThreeVector momentum((*products)[itest].px*MeV,
(*products)[itest].py*MeV,
(*products)[itest].pz*MeV );
Psum += momentum;
theSec->SetMomentum(momentum);
// theResult->AddSecondary(theSec);
theParticleChange.AddSecondary(theSec);
}
// Create heavy fragment if necessary to try to balance A, Z
// Note: this step is only required to prevent warnings at the process level
// where "catastrophic" non-conservation tolerances are set to 1 GeV.
// The residual generated will not necessarily be the one that would
// occur in the actual reaction.
if (iTotA - GAtot > 1) {
theSec = new G4DynamicParticle();
if (iTotZ == GZtot) {
// Special case when a nucleus of only neutrons is requested
// Violate charge conservation and set Z = 1
// G4cout << " Z = 1, A = "<< iTotA - GAtot << " created " << G4endl;
theSec->SetDefinition(G4IonTable::GetIonTable()->GetIon(1, iTotA-GAtot, 0) );
} else {
theSec->SetDefinition(G4IonTable::GetIonTable()->GetIon(iTotZ-GZtot, iTotA-GAtot, 0) );
}
theSec->SetMomentum(projMom - Psum);
// theResult->AddSecondary(theSec);
theParticleChange.AddSecondary(theSec);
}
} // loop OK
delete products;
// theResult->SetStatusChange( stopAndKill );
theParticleChange.SetStatusChange( stopAndKill );
// return theResult;
return &theParticleChange;
}
source/processes/hadronic/models/lend/src/MCGIDI_outputChannel.cc
+144 -110
View File
@@ -205,122 +205,156 @@ double MCGIDI_outputChannel_getFinalQ( statusMessageReporting *smr, MCGIDI_outpu
/*
************************************************************
*/
int MCGIDI_outputChannel_sampleProductsAtE( statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, MCGIDI_quantitiesLookupModes &modes,
MCGIDI_decaySamplingInfo *decaySamplingInfo, MCGIDI_sampledProductsDatas *productDatas, double *masses_ ) {
int MCGIDI_outputChannel_sampleProductsAtE(statusMessageReporting* smr,
MCGIDI_outputChannel* outputChannel,
MCGIDI_quantitiesLookupModes &modes,
MCGIDI_decaySamplingInfo* decaySamplingInfo,
MCGIDI_sampledProductsDatas* productDatas,
double *masses_ )
{
int i1;
int multiplicity(0);
int secondTwoBody = 0, isDecayChannel = ( outputChannel->reaction == NULL );
double e_in = modes.getProjectileEnergy( );
MCGIDI_product *product;
double phi, p, masses[3];
MCGIDI_distribution *distribution;
MCGIDI_sampledProductsData productData[2];
int i1, multiplicity, secondTwoBody = 0, isDecayChannel = ( outputChannel->reaction == NULL );
double e_in = modes.getProjectileEnergy( );
MCGIDI_product *product;
double phi, p, masses[3];
MCGIDI_distribution *distribution;
MCGIDI_sampledProductsData productData[2];
if (isDecayChannel) {
masses[0] = masses_[0]; /* More work may be needed here. */
masses[1] = masses_[1];
} else {
masses[0] = MCGIDI_reaction_getProjectileMass_MeV( smr, outputChannel->reaction );
masses[1] = MCGIDI_reaction_getTargetMass_MeV( smr, outputChannel->reaction );
}
if( isDecayChannel ) {
masses[0] = masses_[0]; /* More work may be needed here. */
masses[1] = masses_[1]; }
else {
masses[0] = MCGIDI_reaction_getProjectileMass_MeV( smr, outputChannel->reaction );
masses[1] = MCGIDI_reaction_getTargetMass_MeV( smr, outputChannel->reaction );
}
// Loop over all possible final state particles reachable from initial state
// List of these particles (products) was read in from GIDI
// Note: all particles satifying the sampling criteria are included in the
// final state, regardless of charge, energy or baryon number conservation
for (i1 = 0; i1 < outputChannel->numberOfProducts; i1++) {
product = &(outputChannel->products[i1]);
if (product->decayChannel.genre != MCGIDI_channelGenre_undefined_e ) {
if( MCGIDI_outputChannel_sampleProductsAtE(smr, &(product->decayChannel),
modes, decaySamplingInfo,
productDatas, masses ) < 0 ) return( -1 );
} else {
distribution = &(product->distribution);
if( distribution->type == MCGIDI_distributionType_none_e ) continue;
for( i1 = 0; i1 < outputChannel->numberOfProducts; i1++ ) {
product = &(outputChannel->products[i1]);
if( product->decayChannel.genre != MCGIDI_channelGenre_undefined_e ) {
if( MCGIDI_outputChannel_sampleProductsAtE( smr, &(product->decayChannel), modes, decaySamplingInfo, productDatas, masses ) < 0 ) return( -1 ); }
else {
distribution = &(product->distribution);
if( distribution->type == MCGIDI_distributionType_none_e ) continue;
if( !secondTwoBody ) {
if( ( multiplicity = product->multiplicity ) == 0 ) multiplicity = MCGIDI_product_sampleMultiplicity( smr, product, e_in,
decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
while( multiplicity > 0 ) {
if (!secondTwoBody) {
// Sample multiplicity of final state particle at kinetic energy of projectile
// The multiplicity stored in GIDI is a real number whose fractional part is
// compared to a random number to decide what integer value is returned
if ((multiplicity = product->multiplicity) == 0) multiplicity =
MCGIDI_product_sampleMultiplicity(smr, product, e_in,
decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
while (multiplicity > 0) {
multiplicity--;
decaySamplingInfo->pop = product->pop;
decaySamplingInfo->mu = 0;
decaySamplingInfo->Ep = 0;
productData[0].isVelocity = decaySamplingInfo->isVelocity;
productData[0].pop = product->pop;
productData[0].delayedNeutronIndex = product->delayedNeutronIndex;
productData[0].delayedNeutronRate = product->delayedNeutronRate;
productData[0].birthTimeSec = 0;
if( product->delayedNeutronRate > 0 ) {
productData[0].birthTimeSec = -G4Log( decaySamplingInfo->rng( decaySamplingInfo->rngState ) ) / product->delayedNeutronRate;
}
multiplicity--;
decaySamplingInfo->pop = product->pop;
decaySamplingInfo->mu = 0;
decaySamplingInfo->Ep = 0;
productData[0].isVelocity = decaySamplingInfo->isVelocity;
productData[0].pop = product->pop;
productData[0].delayedNeutronIndex = product->delayedNeutronIndex;
productData[0].delayedNeutronRate = product->delayedNeutronRate;
productData[0].birthTimeSec = 0;
if (product->delayedNeutronRate > 0) {
productData[0].birthTimeSec =
-G4Log( decaySamplingInfo->rng( decaySamplingInfo->rngState ) ) / product->delayedNeutronRate;
}
switch( outputChannel->genre ) {
case MCGIDI_channelGenre_twoBody_e :
secondTwoBody = 1;
MCGIDI_angular_sampleMu( smr, distribution->angular, modes, decaySamplingInfo );
if( smr_isOk( smr ) ) {
phi = 2. * M_PI * decaySamplingInfo->rng( decaySamplingInfo->rngState );
MCGIDI_kinetics_2BodyReaction( smr, distribution->angular, e_in, decaySamplingInfo->mu, phi, productData );
if( !smr_isOk( smr ) ) return( -1 );
productData[1].pop = product[1].pop;
productData[1].delayedNeutronIndex = product[1].delayedNeutronIndex;
productData[1].delayedNeutronRate = product->delayedNeutronRate;
productData[1].birthTimeSec = 0;
MCGIDI_sampledProducts_addProduct( smr, productDatas, productData );
if( !smr_isOk( smr ) ) return( -1 );
MCGIDI_sampledProducts_addProduct( smr, productDatas, &(productData[1]) );
if( !smr_isOk( smr ) ) return( -1 );
}
break;
case MCGIDI_channelGenre_uncorrelated_e :
case MCGIDI_channelGenre_sumOfRemaining_e :
masses[2] = MCGIDI_product_getMass_MeV( smr, product );
switch( distribution->type ) {
case MCGIDI_distributionType_uncorrelated_e :
MCGIDI_uncorrelated_sampleDistribution( smr, distribution, modes, decaySamplingInfo );
break;
case MCGIDI_distributionType_energyAngular_e :
MCGIDI_energyAngular_sampleDistribution( smr, distribution, modes, decaySamplingInfo );
break;
case MCGIDI_distributionType_KalbachMann_e :
MCGIDI_KalbachMann_sampleEp( smr, distribution->KalbachMann, modes, decaySamplingInfo );
break;
case MCGIDI_distributionType_angularEnergy_e :
MCGIDI_angularEnergy_sampleDistribution( smr, distribution->angularEnergy, modes, decaySamplingInfo );
break;
default :
printf( "Unknown spectral data form product name = %s, channel genre = %d\n", product->pop->name, outputChannel->genre );
break;
}
break;
case MCGIDI_channelGenre_undefined_e :
printf( "Channel is undefined\n" );
break;
case MCGIDI_channelGenre_twoBodyDecay_e :
printf( "Channel is twoBodyDecay\n" );
break;
case MCGIDI_channelGenre_uncorrelatedDecay_e :
printf( "Channel is uncorrelatedDecay\n" );
break;
default :
printf( "Unsupported channel genre = %d\n", outputChannel->genre );
break;
}
if( !smr_isOk( smr ) ) return( -1 );
if( !secondTwoBody ) {
if( decaySamplingInfo->frame == xDataTOM_frame_centerOfMass ) {
if( MCGIDI_kinetics_COM2Lab( smr, modes, decaySamplingInfo, masses ) != 0 ) return( -1 );
}
productData[0].kineticEnergy = decaySamplingInfo->Ep;
p = std::sqrt( decaySamplingInfo->Ep * ( decaySamplingInfo->Ep + 2. * product->pop->mass_MeV ) );
if( productData[0].isVelocity ) p *= MCGIDI_speedOfLight_cm_sec / std::sqrt( p * p + product->pop->mass_MeV * product->pop->mass_MeV );
productData[0].pz_vz = p * decaySamplingInfo->mu;
p = std::sqrt( 1. - decaySamplingInfo->mu * decaySamplingInfo->mu ) * p;
phi = 2. * M_PI * decaySamplingInfo->rng( decaySamplingInfo->rngState );
productData[0].px_vx = p * std::sin( phi );
productData[0].py_vy = p * std::cos( phi );
MCGIDI_sampledProducts_addProduct( smr, productDatas, productData );
if( !smr_isOk( smr ) ) return( -1 );
}
} // Loop checking, 11.06.2015, T. Koi
switch( outputChannel->genre ) {
case MCGIDI_channelGenre_twoBody_e :
secondTwoBody = 1;
MCGIDI_angular_sampleMu( smr, distribution->angular, modes, decaySamplingInfo );
if (smr_isOk(smr) ) {
phi = 2. * M_PI * decaySamplingInfo->rng( decaySamplingInfo->rngState );
MCGIDI_kinetics_2BodyReaction( smr, distribution->angular, e_in, decaySamplingInfo->mu, phi, productData );
if (!smr_isOk(smr) ) return( -1 );
productData[1].pop = product[1].pop;
productData[1].delayedNeutronIndex = product[1].delayedNeutronIndex;
productData[1].delayedNeutronRate = product->delayedNeutronRate;
productData[1].birthTimeSec = 0;
MCGIDI_sampledProducts_addProduct( smr, productDatas, productData );
if( !smr_isOk( smr ) ) return( -1 );
MCGIDI_sampledProducts_addProduct( smr, productDatas, &(productData[1]) );
if( !smr_isOk( smr ) ) return( -1 );
}
}
}
return( productDatas->numberOfProducts );
break;
case MCGIDI_channelGenre_uncorrelated_e :
case MCGIDI_channelGenre_sumOfRemaining_e :
// Get mass of final state particle, then get its distribution
// masses[0] and masses[1] are incident and target masses
masses[2] = MCGIDI_product_getMass_MeV( smr, product );
switch( distribution->type ) {
case MCGIDI_distributionType_uncorrelated_e :
MCGIDI_uncorrelated_sampleDistribution( smr, distribution, modes, decaySamplingInfo );
break;
case MCGIDI_distributionType_energyAngular_e :
MCGIDI_energyAngular_sampleDistribution( smr, distribution, modes, decaySamplingInfo );
break;
case MCGIDI_distributionType_KalbachMann_e :
MCGIDI_KalbachMann_sampleEp( smr, distribution->KalbachMann, modes, decaySamplingInfo );
break;
case MCGIDI_distributionType_angularEnergy_e :
MCGIDI_angularEnergy_sampleDistribution( smr, distribution->angularEnergy, modes, decaySamplingInfo );
break;
default :
printf( "Unknown spectral data form product name = %s, channel genre = %d\n", product->pop->name, outputChannel->genre );
break;
}
break;
case MCGIDI_channelGenre_undefined_e :
printf( "Channel is undefined\n" );
break;
case MCGIDI_channelGenre_twoBodyDecay_e :
printf( "Channel is twoBodyDecay\n" );
break;
case MCGIDI_channelGenre_uncorrelatedDecay_e :
printf( "Channel is uncorrelatedDecay\n" );
break;
default :
printf( "Unsupported channel genre = %d\n", outputChannel->genre );
break;
}
if (!smr_isOk(smr) ) return( -1 );
if (!secondTwoBody) {
if (decaySamplingInfo->frame == xDataTOM_frame_centerOfMass) {
if (MCGIDI_kinetics_COM2Lab( smr, modes, decaySamplingInfo, masses) != 0 ) return( -1 );
}
// Assign kinematics to final state product
productData[0].kineticEnergy = decaySamplingInfo->Ep;
p = std::sqrt( decaySamplingInfo->Ep * ( decaySamplingInfo->Ep + 2. * product->pop->mass_MeV ) );
if (productData[0].isVelocity) p *= MCGIDI_speedOfLight_cm_sec / std::sqrt( p * p + product->pop->mass_MeV * product->pop->mass_MeV );
productData[0].pz_vz = p * decaySamplingInfo->mu;
p = std::sqrt( 1. - decaySamplingInfo->mu * decaySamplingInfo->mu ) * p;
phi = 2. * M_PI * decaySamplingInfo->rng( decaySamplingInfo->rngState );
productData[0].px_vx = p * std::sin( phi );
productData[0].py_vy = p * std::cos( phi );
MCGIDI_sampledProducts_addProduct( smr, productDatas, productData );
if (!smr_isOk(smr) ) return( -1 );
}
} // while multiplicity
} // if !secondTwoBody
} // if decay channel genre
} // loop over possible final state products
return( productDatas->numberOfProducts );
}
#if defined __cplusplus