lars/geant4
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
e2d2f9810a9c69a246ddee17ae224bf9d5ac3453
geant4/source/processes/hadronic/util/src/G4ReactionDynamics.cc
T
Gabriele Cosmo e2d2f9810a Import Geant4 9.6.0 source tree
2016-06-09 17:01:34 +02:00

4074 lines
154 KiB
C++
Raw Blame History

//
// ********************************************************************
// * 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. *
// ********************************************************************
//
//
//
// Hadronic Process: Reaction Dynamics
// original by H.P. Wellisch
// modified by J.L. Chuma, TRIUMF, 19-Nov-1996
// Last modified: 27-Mar-1997
// modified by H.P. Wellisch, 24-Apr-97
// H.P. Wellisch, 25.Apr-97: Side of current and target particle taken into account
// H.P. Wellisch, 29.Apr-97: Bug fix in NuclearReaction. (pseudo1 was without energy)
// J.L. Chuma, 30-Apr-97: Changed return value for GenerateXandPt. It seems possible
// that GenerateXandPt could eliminate some secondaries, but
// still finish its calculations, thus we would not want to
// then use TwoCluster to again calculate the momenta if vecLen
// was less than 6.
// J.L. Chuma, 10-Jun-97: Modified NuclearReaction. Was not creating ReactionProduct's
// with the new operator, thus they would be meaningless when
// going out of scope.
// J.L. Chuma, 20-Jun-97: Modified GenerateXandPt and TwoCluster to fix units problems
// J.L. Chuma, 23-Jun-97: Modified ProduceStrangeParticlePairs to fix units problems
// J.L. Chuma, 26-Jun-97: Modified ProduceStrangeParticlePairs to fix array indices
// which were sometimes going out of bounds
// J.L. Chuma, 04-Jul-97: Many minor modifications to GenerateXandPt and TwoCluster
// J.L. Chuma, 06-Aug-97: Added original incident particle, before Fermi motion and
// evaporation effects are included, needed for self absorption
// and corrections for single particle spectra (shower particles)
// logging stopped 1997
// J. Allison, 17-Jun-99: Replaced a min function to get correct behaviour on DEC.
#include <iostream>
#include <signal.h>
#include "G4ReactionDynamics.hh"
#include "G4AntiProton.hh"
#include "G4AntiNeutron.hh"
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "Randomize.hh"
#include "G4HadReentrentException.hh"
// #include "DumpFrame.hh"
/* G4double GetQValue(G4ReactionProduct * aSec)
{
double QValue=0;
if(aSec->GetDefinition()->GetParticleType() == "baryon")
{
if(aSec->GetDefinition()->GetBaryonNumber() < 0)
{
QValue = aSec->GetTotalEnergy();
QValue += G4Neutron::Neutron()->GetPDGMass();
}
else
{
G4double ss = 0;
ss +=aSec->GetDefinition()->GetPDGMass();
if(aSec->GetDefinition() == G4Proton::Proton())
{
ss -=G4Proton::Proton()->GetPDGMass();
}
else
{
ss -=G4Neutron::Neutron()->GetPDGMass();
}
ss += aSec->GetKineticEnergy();
QValue = ss;
}
}
else if(aSec->GetDefinition()->GetPDGEncoding() == 0)
{
QValue = aSec->GetKineticEnergy();
}
else
{
QValue = aSec->GetTotalEnergy();
}
return QValue;
}
*/
G4bool G4ReactionDynamics::GenerateXandPt(
G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
G4int &vecLen,
G4ReactionProduct &modifiedOriginal, // Fermi motion & evap. effects included
const G4HadProjectile *originalIncident, // the original incident particle
G4ReactionProduct &currentParticle,
G4ReactionProduct &targetParticle,
const G4DynamicParticle* originalTarget,
const G4Nucleus &targetNucleus,
G4bool &incidentHasChanged,
G4bool &targetHasChanged,
G4bool leadFlag,
G4ReactionProduct &leadingStrangeParticle )
{
//
// derived from original FORTRAN code GENXPT by H. Fesefeldt (11-Oct-1987)
//
// Generation of X- and PT- values for incident, target, and all secondary particles
// A simple single variable description E D3S/DP3= F(Q) with
// Q^2 = (M*X)^2 + PT^2 is used. Final state kinematic is produced
// by an FF-type iterative cascade method
//
// internal units are GeV
//
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
// Protection in case no secondary has been created; cascades down to two-body.
if(vecLen == 0) return false;
G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
G4ParticleDefinition *aProton = G4Proton::Proton();
G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
G4int i, l;
G4bool veryForward = false;
const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
targetMass*targetMass +
2.0*targetMass*etOriginal ); // GeV
G4double currentMass = currentParticle.GetMass()/GeV;
targetMass = targetParticle.GetMass()/GeV;
//
// randomize the order of the secondary particles
// note that the current and target particles are not affected
//
for( i=0; i<vecLen; ++i )
{
G4int itemp = G4int( G4UniformRand()*vecLen );
G4ReactionProduct pTemp = *vec[itemp];
*vec[itemp] = *vec[i];
*vec[i] = pTemp;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
if( currentMass == 0.0 && targetMass == 0.0 )
{
// Target and projectile have annihilated. Replace them with the first
// two secondaries in the list. Current particle KE is maintained.
G4double ek = currentParticle.GetKineticEnergy();
G4ThreeVector mom = currentParticle.GetMomentum();
currentParticle = *vec[0];
targetParticle = *vec[1];
for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
G4ReactionProduct *temp = vec[vecLen-1];
delete temp;
temp = vec[vecLen-2];
delete temp;
vecLen -= 2;
currentMass = currentParticle.GetMass()/GeV;
targetMass = targetParticle.GetMass()/GeV;
incidentHasChanged = true;
targetHasChanged = true;
currentParticle.SetKineticEnergy( ek );
currentParticle.SetMomentum( mom );
veryForward = true;
}
const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
const G4double protonMass = aProton->GetPDGMass()/MeV;
if (originalIncident->GetDefinition()->GetParticleSubType() == "kaon"
&& G4UniformRand() >= 0.7) {
G4ReactionProduct temp = currentParticle;
currentParticle = targetParticle;
targetParticle = temp;
incidentHasChanged = true;
targetHasChanged = true;
currentMass = currentParticle.GetMass()/GeV;
targetMass = targetParticle.GetMass()/GeV;
}
const G4double afc = std::min( 0.75,
0.312+0.200*std::log(std::log(centerofmassEnergy*centerofmassEnergy))+
std::pow(centerofmassEnergy*centerofmassEnergy,1.5)/6000.0 );
G4double freeEnergy = centerofmassEnergy-currentMass-targetMass;
G4double forwardEnergy = freeEnergy/2.;
G4int forwardCount = 1; // number of particles in forward hemisphere
G4double backwardEnergy = freeEnergy/2.;
G4int backwardCount = 1; // number of particles in backward hemisphere
if(veryForward)
{
if(currentParticle.GetSide()==-1)
{
forwardEnergy += currentMass;
forwardCount --;
backwardEnergy -= currentMass;
backwardCount ++;
}
if(targetParticle.GetSide()!=-1)
{
backwardEnergy += targetMass;
backwardCount --;
forwardEnergy -= targetMass;
forwardCount ++;
}
}
for( i=0; i<vecLen; ++i )
{
if( vec[i]->GetSide() == -1 )
{
++backwardCount;
backwardEnergy -= vec[i]->GetMass()/GeV;
} else {
++forwardCount;
forwardEnergy -= vec[i]->GetMass()/GeV;
}
}
//
// Add particles from intranuclear cascade.
// nuclearExcitationCount = number of new secondaries produced by nuclear excitation
// extraCount = number of nucleons within these new secondaries
//
G4double xtarg;
if( centerofmassEnergy < (2.0+G4UniformRand()) )
xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount+vecLen+2)/2.0;
else
xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount);
if( xtarg <= 0.0 )xtarg = 0.01;
G4int nuclearExcitationCount = Poisson( xtarg );
if(atomicWeight<1.0001) nuclearExcitationCount = 0;
G4int extraNucleonCount = 0;
G4double extraNucleonMass = 0.0;
if( nuclearExcitationCount > 0 )
{
const G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
const G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
G4int momentumBin = 0;
while( (momentumBin < 6) &&
(modifiedOriginal.GetTotalMomentum()/GeV > psup[momentumBin]) )
++momentumBin;
momentumBin = std::min( 5, momentumBin );
//
// NOTE: in GENXPT, these new particles were given negative codes
// here I use NewlyAdded = true instead
//
for( i=0; i<nuclearExcitationCount; ++i )
{
G4ReactionProduct * pVec = new G4ReactionProduct();
if( G4UniformRand() < nucsup[momentumBin] )
{
if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
pVec->SetDefinition( aProton );
else
pVec->SetDefinition( aNeutron );
pVec->SetSide( -2 ); // -2 means backside nucleon
++extraNucleonCount;
backwardEnergy += pVec->GetMass()/GeV;
extraNucleonMass += pVec->GetMass()/GeV;
}
else
{
G4double ran = G4UniformRand();
if( ran < 0.3181 )
pVec->SetDefinition( aPiPlus );
else if( ran < 0.6819 )
pVec->SetDefinition( aPiZero );
else
pVec->SetDefinition( aPiMinus );
pVec->SetSide( -1 ); // backside particle, but not a nucleon
}
pVec->SetNewlyAdded( true ); // true is the same as IPA(i)<0
vec.SetElement( vecLen++, pVec );
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
backwardEnergy -= pVec->GetMass()/GeV;
++backwardCount;
}
}
//
// assume conservation of kinetic energy in forward & backward hemispheres
//
G4int is, iskip;
while( forwardEnergy <= 0.0 ) // must eliminate a particle from the forward side
{
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
iskip = G4int(G4UniformRand()*forwardCount) + 1; // 1 <= iskip <= forwardCount
is = 0;
G4int forwardParticlesLeft = 0;
for( i=(vecLen-1); i>=0; --i )
{
if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
{
forwardParticlesLeft = 1;
if( ++is == iskip )
{
forwardEnergy += vec[i]->GetMass()/GeV;
for( G4int j=i; j<(vecLen-1); j++ )*vec[j] = *vec[j+1]; // shift up
--forwardCount;
delete vec[vecLen-1];
if( --vecLen == 0 )return false; // all the secondaries have been eliminated
break; // --+
} // |
} // |
} // break goes down to here
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
if( forwardParticlesLeft == 0 )
{
G4int iremove = -1;
for (i = 0; i < vecLen; i++) {
if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
iremove = i;
break;
}
}
if (iremove == -1) {
for (i = 0; i < vecLen; i++) {
if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
iremove = i;
break;
}
}
}
if (iremove == -1) iremove = 0;
forwardEnergy += vec[iremove]->GetMass()/GeV;
if (vec[iremove]->GetSide() > 0) --forwardCount;
for (i = iremove; i < vecLen-1; i++) *vec[i] = *vec[i+1];
delete vec[vecLen-1];
vecLen--;
if (vecLen == 0) return false; // all secondaries have been eliminated
break;
}
} // while
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
while( backwardEnergy <= 0.0 ) // must eliminate a particle from the backward side
{
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
iskip = G4int(G4UniformRand()*backwardCount) + 1; // 1 <= iskip <= backwardCount
is = 0;
G4int backwardParticlesLeft = 0;
for( i=(vecLen-1); i>=0; --i )
{
if( vec[i]->GetSide() < 0 && vec[i]->GetMayBeKilled())
{
backwardParticlesLeft = 1;
if( ++is == iskip ) // eliminate the i'th particle
{
if( vec[i]->GetSide() == -2 )
{
--extraNucleonCount;
extraNucleonMass -= vec[i]->GetMass()/GeV;
backwardEnergy -= vec[i]->GetTotalEnergy()/GeV;
}
backwardEnergy += vec[i]->GetTotalEnergy()/GeV;
for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
--backwardCount;
delete vec[vecLen-1];
if( --vecLen == 0 )return false; // all the secondaries have been eliminated
break;
}
}
}
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
if( backwardParticlesLeft == 0 )
{
G4int iremove = -1;
for (i = 0; i < vecLen; i++) {
if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
iremove = i;
break;
}
}
if (iremove == -1) {
for (i = 0; i < vecLen; i++) {
if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
iremove = i;
break;
}
}
}
if (iremove == -1) iremove = 0;
backwardEnergy += vec[iremove]->GetMass()/GeV;
if (vec[iremove]->GetSide() > 0) --backwardCount;
for (i = iremove; i < vecLen-1; i++) *vec[i] = *vec[i+1];
delete vec[vecLen-1];
vecLen--;
if (vecLen == 0) return false; // all secondaries have been eliminated
break;
}
} // while
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
//
// define initial state vectors for Lorentz transformations
// the pseudoParticles have non-standard masses, hence the "pseudo"
//
G4ReactionProduct pseudoParticle[10];
for( i=0; i<10; ++i )pseudoParticle[i].SetZero();
pseudoParticle[0].SetMass( mOriginal*GeV );
pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
pseudoParticle[0].SetTotalEnergy(
std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
pseudoParticle[1].SetMass( protonMass*MeV ); // this could be targetMass
pseudoParticle[1].SetTotalEnergy( protonMass*MeV );
pseudoParticle[3].SetMass( protonMass*(1+extraNucleonCount)*MeV );
pseudoParticle[3].SetTotalEnergy( protonMass*(1+extraNucleonCount)*MeV );
pseudoParticle[8].SetMomentum( 1.0*GeV, 0.0, 0.0 );
pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
G4double dndl[20];
//
// main loop for 4-momentum generation
// see Pitha-report (Aachen) for a detailed description of the method
//
G4double aspar, pt, et, x, pp, pp1, rmb, wgt;
G4int innerCounter, outerCounter;
G4bool eliminateThisParticle, resetEnergies, constantCrossSection;
G4double forwardKinetic = 0.0, backwardKinetic = 0.0;
//
// process the secondary particles in reverse order
// the incident particle is Done after the secondaries
// nucleons, including the target, in the backward hemisphere are also Done later
//
G4double binl[20] = {0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.11,1.25,
1.43,1.67,2.0,2.5,3.33,5.00,10.00};
G4int backwardNucleonCount = 0; // number of nucleons in backward hemisphere
G4double totalEnergy, kineticEnergy, vecMass;
for( i=(vecLen-1); i>=0; --i )
{
G4double phi = G4UniformRand()*twopi;
if( vec[i]->GetNewlyAdded() ) // added from intranuclear cascade
{
if( vec[i]->GetSide() == -2 ) // is a nucleon
{
if( backwardNucleonCount < 18 )
{
if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
for(G4int j=0; j<vecLen; j++) delete vec[j];
vecLen = 0;
throw G4HadReentrentException(__FILE__, __LINE__,
"G4ReactionDynamics::GenerateXandPt : a pion has been counted as a backward nucleon");
}
vec[i]->SetSide( -3 );
++backwardNucleonCount;
continue;
}
}
}
//
// set pt and phi values, they are changed somewhat in the iteration loop
// set mass parameter for lambda fragmentation model
//
vecMass = vec[i]->GetMass()/GeV;
G4double ran = -std::log(1.0-G4UniformRand())/3.5;
if( vec[i]->GetSide() == -2 ) // backward nucleon
{
if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon" ||
vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
aspar = 0.75;
pt = std::sqrt( std::pow( ran, 1.7 ) );
} else { // vec[i] must be a proton, neutron,
aspar = 0.20; // lambda, sigma, xsi, or ion
pt = std::sqrt( std::pow( ran, 1.2 ) );
}
} else { // not a backward nucleon
if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
aspar = 0.75;
pt = std::sqrt( std::pow( ran, 1.7 ) );
} else if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
aspar = 0.70;
pt = std::sqrt( std::pow( ran, 1.7 ) );
} else { // vec[i] must be a proton, neutron,
aspar = 0.65; // lambda, sigma, xsi, or ion
pt = std::sqrt( std::pow( ran, 1.5 ) );
}
}
pt = std::max( 0.001, pt );
vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
for( G4int j=0; j<20; ++j )binl[j] = j/(19.*pt);
if( vec[i]->GetSide() > 0 )
et = pseudoParticle[0].GetTotalEnergy()/GeV;
else
et = pseudoParticle[1].GetTotalEnergy()/GeV;
dndl[0] = 0.0;
//
// start of outer iteration loop
//
outerCounter = 0;
eliminateThisParticle = true;
resetEnergies = true;
while( ++outerCounter < 3 )
{
for( l=1; l<20; ++l )
{
x = (binl[l]+binl[l-1])/2.;
pt = std::max( 0.001, pt );
if( x > 1.0/pt )
dndl[l] += dndl[l-1]; // changed from just = on 02 April 98
else
dndl[l] = et * aspar/std::sqrt( std::pow((1.+aspar*x*aspar*x),3) )
* (binl[l]-binl[l-1]) / std::sqrt( pt*x*et*pt*x*et + pt*pt + vecMass*vecMass )
+ dndl[l-1];
}
innerCounter = 0;
vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
//
// start of inner iteration loop
//
while( ++innerCounter < 7 )
{
ran = G4UniformRand()*dndl[19];
l = 1;
while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
l = std::min( 19, l );
x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
if( vec[i]->GetSide() < 0 )x *= -1.;
vec[i]->SetMomentum( x*et*GeV ); // set the z-momentum
totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
vec[i]->SetTotalEnergy( totalEnergy*GeV );
kineticEnergy = vec[i]->GetKineticEnergy()/GeV;
if( vec[i]->GetSide() > 0 ) // forward side
{
if( (forwardKinetic+kineticEnergy) < 0.95*forwardEnergy )
{
pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
forwardKinetic += kineticEnergy;
pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
pseudoParticle[6].SetMomentum( 0.0 ); // set the z-momentum
phi = pseudoParticle[6].Angle( pseudoParticle[8] );
if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
phi += pi + normal()*pi/12.0;
if( phi > twopi )phi -= twopi;
if( phi < 0.0 )phi = twopi - phi;
outerCounter = 2; // leave outer loop
eliminateThisParticle = false; // don't eliminate this particle
resetEnergies = false;
break; // leave inner loop
}
if( innerCounter > 5 )break; // leave inner loop
if( backwardEnergy >= vecMass ) // switch sides
{
vec[i]->SetSide( -1 );
forwardEnergy += vecMass;
backwardEnergy -= vecMass;
++backwardCount;
}
} else { // backward side
if( extraNucleonCount > 19 ) // commented out to duplicate ?bug? in GENXPT
x = 0.999;
G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
if( (backwardKinetic+kineticEnergy) < xxx*backwardEnergy )
{
pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
backwardKinetic += kineticEnergy;
pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
pseudoParticle[6].SetMomentum( 0.0 ); // set the z-momentum
phi = pseudoParticle[6].Angle( pseudoParticle[8] );
if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
phi += pi + normal() * pi / 12.0;
if( phi > twopi )phi -= twopi;
if( phi < 0.0 )phi = twopi - phi;
outerCounter = 2; // leave outer loop
eliminateThisParticle = false; // don't eliminate this particle
resetEnergies = false;
break; // leave inner loop
}
if( innerCounter > 5 )break; // leave inner loop
if( forwardEnergy >= vecMass ) // switch sides
{
vec[i]->SetSide( 1 );
forwardEnergy -= vecMass;
backwardEnergy += vecMass;
backwardCount--;
}
}
G4ThreeVector momentum = vec[i]->GetMomentum();
vec[i]->SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
pt *= 0.9;
dndl[19] *= 0.9;
} // closes inner loop
if( resetEnergies )
{
// if we get to here, the inner loop has been Done 6 Times
// reset the kinetic energies of previously Done particles, if they are lighter
// than protons and in the forward hemisphere
// then continue with outer loop
//
forwardKinetic = 0.0;
backwardKinetic = 0.0;
pseudoParticle[4].SetZero();
pseudoParticle[5].SetZero();
for( l=i+1; l<vecLen; ++l )
{
if (vec[l]->GetSide() > 0 ||
vec[l]->GetDefinition()->GetParticleSubType() == "kaon" ||
vec[l]->GetDefinition()->GetParticleSubType() == "pi") {
G4double tempMass = vec[l]->GetMass()/MeV;
totalEnergy = 0.95*vec[l]->GetTotalEnergy()/MeV + 0.05*tempMass;
totalEnergy = std::max( tempMass, totalEnergy );
vec[l]->SetTotalEnergy( totalEnergy*MeV );
pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) );
pp1 = vec[l]->GetMomentum().mag()/MeV;
if( pp1 < 1.0e-6*GeV )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi2 = twopi*G4UniformRand();
vec[l]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
pp*sintheta*std::sin(phi2)*MeV,
pp*costheta*MeV ) ;
} else {
vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );
}
G4double px = vec[l]->GetMomentum().x()/MeV;
G4double py = vec[l]->GetMomentum().y()/MeV;
pt = std::max( 1.0, std::sqrt( px*px + py*py ) )/GeV;
if( vec[l]->GetSide() > 0 )
{
forwardKinetic += vec[l]->GetKineticEnergy()/GeV;
pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
} else {
backwardKinetic += vec[l]->GetKineticEnergy()/GeV;
pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
}
} // if pi, K or forward
} // for l
} // if resetEnergies
} // closes outer loop
if( eliminateThisParticle && vec[i]->GetMayBeKilled()) // not enough energy, eliminate this particle
{
if( vec[i]->GetSide() > 0 )
{
--forwardCount;
forwardEnergy += vecMass;
} else {
if( vec[i]->GetSide() == -2 )
{
--extraNucleonCount;
extraNucleonMass -= vecMass;
backwardEnergy -= vecMass;
}
--backwardCount;
backwardEnergy += vecMass;
}
for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
G4ReactionProduct *temp = vec[vecLen-1];
delete temp;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
if( --vecLen == 0 )return false; // all the secondaries have been eliminated
pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
pseudoParticle[6].SetMomentum( 0.0 ); // set z-momentum
phi = pseudoParticle[6].Angle( pseudoParticle[8] );
if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
phi += pi + normal() * pi / 12.0;
if( phi > twopi )phi -= twopi;
if( phi < 0.0 )phi = twopi - phi;
}
} // closes main for loop
//
// for the incident particle: it was placed in the forward hemisphere
// set pt and phi values, they are changed somewhat in the iteration loop
// set mass parameter for lambda fragmentation model
//
G4double phi = G4UniformRand()*twopi;
G4double ran = -std::log(1.0-G4UniformRand());
if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
aspar = 0.60;
pt = std::sqrt( std::pow( ran/6.0, 1.7 ) );
} else if (currentParticle.GetDefinition()->GetParticleSubType() == "kaon") {
aspar = 0.50;
pt = std::sqrt( std::pow( ran/5.0, 1.4 ) );
} else {
aspar = 0.40;
pt = std::sqrt( std::pow( ran/4.0, 1.2 ) );
}
for( G4int j=0; j<20; ++j )binl[j] = j/(19.*pt);
currentParticle.SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
et = pseudoParticle[0].GetTotalEnergy()/GeV;
dndl[0] = 0.0;
vecMass = currentParticle.GetMass()/GeV;
for( l=1; l<20; ++l )
{
x = (binl[l]+binl[l-1])/2.;
if( x > 1.0/pt )
dndl[l] += dndl[l-1]; // changed from just = on 02 April 98
else
dndl[l] = aspar/std::sqrt( std::pow((1.+sqr(aspar*x)),3) ) *
(binl[l]-binl[l-1]) * et / std::sqrt( pt*x*et*pt*x*et + pt*pt + vecMass*vecMass ) +
dndl[l-1];
}
ran = G4UniformRand()*dndl[19];
l = 1;
while( (ran>dndl[l]) && (l<20) )l++;
l = std::min( 19, l );
x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
currentParticle.SetMomentum( x*et*GeV ); // set the z-momentum
if( forwardEnergy < forwardKinetic )
totalEnergy = vecMass + 0.04*std::fabs(normal());
else
totalEnergy = vecMass + forwardEnergy - forwardKinetic;
currentParticle.SetTotalEnergy( totalEnergy*GeV );
pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
pp1 = currentParticle.GetMomentum().mag()/MeV;
if( pp1 < 1.0e-6*GeV )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi2 = twopi*G4UniformRand();
currentParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
pp*sintheta*std::sin(phi2)*MeV,
pp*costheta*MeV ) ;
} else {
currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
}
pseudoParticle[4] = pseudoParticle[4] + currentParticle;
//
// Current particle now finished
//
// Begin target particle
//
if( backwardNucleonCount < 18 )
{
targetParticle.SetSide( -3 );
++backwardNucleonCount;
}
else
{
// set pt and phi values, they are changed somewhat in the iteration loop
// set mass parameter for lambda fragmentation model
//
vecMass = targetParticle.GetMass()/GeV;
ran = -std::log(1.0-G4UniformRand());
aspar = 0.40;
pt = std::max( 0.001, std::sqrt( std::pow( ran/4.0, 1.2 ) ) );
targetParticle.SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
for( G4int j=0; j<20; ++j )binl[j] = (j-1.)/(19.*pt);
et = pseudoParticle[1].GetTotalEnergy()/GeV;
dndl[0] = 0.0;
outerCounter = 0;
eliminateThisParticle = true; // should never eliminate the target particle
resetEnergies = true;
while( ++outerCounter < 3 ) // start of outer iteration loop
{
for( l=1; l<20; ++l )
{
x = (binl[l]+binl[l-1])/2.;
if( x > 1.0/pt )
dndl[l] += dndl[l-1]; // changed from just = on 02 April 98
else
dndl[l] = aspar/std::sqrt( std::pow((1.+aspar*x*aspar*x),3) ) *
(binl[l]-binl[l-1])*et / std::sqrt( pt*x*et*pt*x*et+pt*pt+vecMass*vecMass ) +
dndl[l-1];
}
innerCounter = 0;
while( ++innerCounter < 7 ) // start of inner iteration loop
{
l = 1;
ran = G4UniformRand()*dndl[19];
while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
l = std::min( 19, l );
x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
if( targetParticle.GetSide() < 0 )x *= -1.;
targetParticle.SetMomentum( x*et*GeV ); // set the z-momentum
totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
targetParticle.SetTotalEnergy( totalEnergy*GeV );
if( targetParticle.GetSide() < 0 )
{
if( extraNucleonCount > 19 )x=0.999;
G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
if( (backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy )
{
pseudoParticle[5] = pseudoParticle[5] + targetParticle;
backwardKinetic += totalEnergy - vecMass;
pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
pseudoParticle[6].SetMomentum( 0.0 ); // set z-momentum
phi = pseudoParticle[6].Angle( pseudoParticle[8] );
if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
phi += pi + normal() * pi / 12.0;
if( phi > twopi )phi -= twopi;
if( phi < 0.0 )phi = twopi - phi;
outerCounter = 2; // leave outer loop
eliminateThisParticle = false; // don't eliminate this particle
resetEnergies = false;
break; // leave inner loop
}
if( innerCounter > 5 )break; // leave inner loop
if( forwardEnergy >= vecMass ) // switch sides
{
targetParticle.SetSide( 1 );
forwardEnergy -= vecMass;
backwardEnergy += vecMass;
--backwardCount;
}
G4ThreeVector momentum = targetParticle.GetMomentum();
targetParticle.SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
pt *= 0.9;
dndl[19] *= 0.9;
}
else // target has gone to forward side
{
if( forwardEnergy < forwardKinetic )
totalEnergy = vecMass + 0.04*std::fabs(normal());
else
totalEnergy = vecMass + forwardEnergy - forwardKinetic;
targetParticle.SetTotalEnergy( totalEnergy*GeV );
pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
pp1 = targetParticle.GetMomentum().mag()/MeV;
if( pp1 < 1.0e-6*GeV )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi2 = twopi*G4UniformRand();
targetParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
pp*sintheta*std::sin(phi2)*MeV,
pp*costheta*MeV ) ;
}
else
targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
pseudoParticle[4] = pseudoParticle[4] + targetParticle;
outerCounter = 2; // leave outer loop
eliminateThisParticle = false; // don't eliminate this particle
resetEnergies = false;
break; // leave inner loop
}
} // closes inner loop
if( resetEnergies )
{
// if we get to here, the inner loop has been Done 6 Times
// reset the kinetic energies of previously Done particles, if they are lighter
// than protons and in the forward hemisphere
// then continue with outer loop
forwardKinetic = backwardKinetic = 0.0;
pseudoParticle[4].SetZero();
pseudoParticle[5].SetZero();
for( l=0; l<vecLen; ++l ) // changed from l=1 on 02 April 98
{
if (vec[l]->GetSide() > 0 ||
vec[l]->GetDefinition()->GetParticleSubType() == "kaon" ||
vec[l]->GetDefinition()->GetParticleSubType() == "pi") {
G4double tempMass = vec[l]->GetMass()/GeV;
totalEnergy =
std::max( tempMass, 0.95*vec[l]->GetTotalEnergy()/GeV + 0.05*tempMass );
vec[l]->SetTotalEnergy( totalEnergy*GeV );
pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) )*GeV;
pp1 = vec[l]->GetMomentum().mag()/MeV;
if( pp1 < 1.0e-6*GeV )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi2 = twopi*G4UniformRand();
vec[l]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
pp*sintheta*std::sin(phi2)*MeV,
pp*costheta*MeV ) ;
}
else
vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );
pt = std::max( 0.001*GeV, std::sqrt( sqr(vec[l]->GetMomentum().x()/MeV) +
sqr(vec[l]->GetMomentum().y()/MeV) ) )/GeV;
if( vec[l]->GetSide() > 0)
{
forwardKinetic += vec[l]->GetKineticEnergy()/GeV;
pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
} else {
backwardKinetic += vec[l]->GetKineticEnergy()/GeV;
pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
}
} // if pi, K or forward
} // for l
} // if (resetEnergies)
} // closes outer loop
}
//
// Target particle finished.
//
// Now produce backward nucleons with a cluster model
//
pseudoParticle[6].Lorentz( pseudoParticle[3], pseudoParticle[2] );
pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
if( backwardNucleonCount == 1 ) // target particle is the only backward nucleon
{
G4double ekin =
std::min( backwardEnergy-backwardKinetic, centerofmassEnergy/2.0-protonMass/GeV );
if( ekin < 0.04 )ekin = 0.04 * std::fabs( normal() );
vecMass = targetParticle.GetMass()/GeV;
totalEnergy = ekin+vecMass;
targetParticle.SetTotalEnergy( totalEnergy*GeV );
pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
pp1 = pseudoParticle[6].GetMomentum().mag()/MeV;
if( pp1 < 1.0e-6*GeV )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi2 = twopi*G4UniformRand();
targetParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
pp*sintheta*std::sin(phi2)*MeV,
pp*costheta*MeV ) ;
} else {
targetParticle.SetMomentum( pseudoParticle[6].GetMomentum() * (pp/pp1) );
}
pseudoParticle[5] = pseudoParticle[5] + targetParticle;
}
else // more than one backward nucleon
{
const G4double cpar[] = { 0.6, 0.6, 0.35, 0.15, 0.10 };
const G4double gpar[] = { 2.6, 2.6, 1.80, 1.30, 1.20 };
// Replaced the following min function to get correct behaviour on DEC.
// G4int tempCount = std::min( 5, backwardNucleonCount ) - 1;
G4int tempCount;
if (backwardNucleonCount < 5)
{
tempCount = backwardNucleonCount;
}
else
{
tempCount = 5;
}
tempCount--;
//cout << "backwardNucleonCount " << backwardNucleonCount << G4endl;
//cout << "tempCount " << tempCount << G4endl;
G4double rmb0 = 0.0;
if( targetParticle.GetSide() == -3 )
rmb0 += targetParticle.GetMass()/GeV;
for( i=0; i<vecLen; ++i )
{
if( vec[i]->GetSide() == -3 )rmb0 += vec[i]->GetMass()/GeV;
}
rmb = rmb0 + std::pow(-std::log(1.0-G4UniformRand()),cpar[tempCount]) / gpar[tempCount];
totalEnergy = pseudoParticle[6].GetTotalEnergy()/GeV;
vecMass = std::min( rmb, totalEnergy );
pseudoParticle[6].SetMass( vecMass*GeV );
pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
pp1 = pseudoParticle[6].GetMomentum().mag()/MeV;
if( pp1 < 1.0e-6*GeV )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi2 = twopi*G4UniformRand();
pseudoParticle[6].SetMomentum( -pp*sintheta*std::cos(phi2)*MeV,
-pp*sintheta*std::sin(phi2)*MeV,
-pp*costheta*MeV ) ;
}
else
pseudoParticle[6].SetMomentum( pseudoParticle[6].GetMomentum() * (-pp/pp1) );
G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV; // tempV contains the backward nucleons
tempV.Initialize( backwardNucleonCount );
G4int tempLen = 0;
if( targetParticle.GetSide() == -3 )tempV.SetElement( tempLen++, &targetParticle );
for( i=0; i<vecLen; ++i )
{
if( vec[i]->GetSide() == -3 )tempV.SetElement( tempLen++, vec[i] );
}
if( tempLen != backwardNucleonCount )
{
G4cerr << "tempLen is not the same as backwardNucleonCount" << G4endl;
G4cerr << "tempLen = " << tempLen;
G4cerr << ", backwardNucleonCount = " << backwardNucleonCount << G4endl;
G4cerr << "targetParticle side = " << targetParticle.GetSide() << G4endl;
G4cerr << "currentParticle side = " << currentParticle.GetSide() << G4endl;
for( i=0; i<vecLen; ++i )
G4cerr << "particle #" << i << " side = " << vec[i]->GetSide() << G4endl;
G4Exception("G4ReactionDynamics::GenerateXandPt", "601",
FatalException, "Mismatch in nucleon count");
}
constantCrossSection = true;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
if( tempLen >= 2 )
{
wgt = GenerateNBodyEvent(
pseudoParticle[6].GetMass(), constantCrossSection, tempV, tempLen );
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
if( targetParticle.GetSide() == -3 )
{
targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
// tempV contains the real stuff
pseudoParticle[5] = pseudoParticle[5] + targetParticle;
}
for( i=0; i<vecLen; ++i )
{
if( vec[i]->GetSide() == -3 )
{
vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
}
}
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
}
//
// Lorentz transformation in lab system
//
if( vecLen == 0 )return false; // all the secondaries have been eliminated
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[1] );
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
//
// leadFlag will be true
// iff original particle is at least as heavy as K+ and not a proton or
// neutron AND if incident particle is at least as heavy as K+ and it is
// not a proton or neutron leadFlag is set to the incident particle
// or
// target particle is at least as heavy as K+ and it is not a proton or
// neutron leadFlag is set to the target particle
//
G4bool leadingStrangeParticleHasChanged = true;
if( leadFlag )
{
if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
leadingStrangeParticleHasChanged = false;
if( leadingStrangeParticleHasChanged &&
( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() ) )
leadingStrangeParticleHasChanged = false;
if( leadingStrangeParticleHasChanged )
{
for( i=0; i<vecLen; i++ )
{
if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
{
leadingStrangeParticleHasChanged = false;
break;
}
}
}
if( leadingStrangeParticleHasChanged )
{
G4bool leadTest =
(leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "pi");
G4bool targetTest =
(targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
targetParticle.GetDefinition()->GetParticleSubType() == "pi");
// following modified by JLC 22-Oct-97
if( (leadTest&&targetTest) || !(leadTest||targetTest) ) // both true or both false
{
targetParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
targetHasChanged = true;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
else
{
currentParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
incidentHasChanged = false;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
}
} // end of if( leadFlag )
// Get number of final state nucleons and nucleons remaining in
// target nucleus
std::pair<G4int, G4int> finalStateNucleons =
GetFinalStateNucleons(originalTarget, vec, vecLen);
G4int protonsInFinalState = finalStateNucleons.first;
G4int neutronsInFinalState = finalStateNucleons.second;
G4int numberofFinalStateNucleons =
protonsInFinalState + neutronsInFinalState;
if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
originalIncident->GetDefinition()->GetPDGMass() <
G4Lambda::Lambda()->GetPDGMass())
numberofFinalStateNucleons++;
numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
G4int PinNucleus = std::max(0,
targetNucleus.GetZ_asInt() - protonsInFinalState);
G4int NinNucleus = std::max(0,
targetNucleus.GetN_asInt() - neutronsInFinalState);
pseudoParticle[3].SetMomentum( 0.0, 0.0, pOriginal*GeV );
pseudoParticle[3].SetMass( mOriginal*GeV );
pseudoParticle[3].SetTotalEnergy(
std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
G4int diff = 0;
if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
if(numberofFinalStateNucleons == 1) diff = 0;
pseudoParticle[4].SetMomentum( 0.0, 0.0, 0.0 );
pseudoParticle[4].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV );
pseudoParticle[4].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV );
G4double theoreticalKinetic =
pseudoParticle[3].GetTotalEnergy()/MeV +
pseudoParticle[4].GetTotalEnergy()/MeV -
currentParticle.GetMass()/MeV -
targetParticle.GetMass()/MeV;
G4double simulatedKinetic =
currentParticle.GetKineticEnergy()/MeV + targetParticle.GetKineticEnergy()/MeV;
pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
pseudoParticle[3].Lorentz( pseudoParticle[3], pseudoParticle[5] );
pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[5] );
pseudoParticle[7].SetZero();
pseudoParticle[7] = pseudoParticle[7] + currentParticle;
pseudoParticle[7] = pseudoParticle[7] + targetParticle;
for( i=0; i<vecLen; ++i )
{
pseudoParticle[7] = pseudoParticle[7] + *vec[i];
simulatedKinetic += vec[i]->GetKineticEnergy()/MeV;
theoreticalKinetic -= vec[i]->GetMass()/MeV;
}
if( vecLen <= 16 && vecLen > 0 )
{
// must create a new set of ReactionProducts here because GenerateNBody will
// modify the momenta for the particles, and we don't want to do this
//
G4ReactionProduct tempR[130];
tempR[0] = currentParticle;
tempR[1] = targetParticle;
for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
tempV.Initialize( vecLen+2 );
G4int tempLen = 0;
for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
constantCrossSection = true;
wgt = GenerateNBodyEvent( pseudoParticle[3].GetTotalEnergy()/MeV+
pseudoParticle[4].GetTotalEnergy()/MeV,
constantCrossSection, tempV, tempLen );
if (wgt == -1) {
G4double Qvalue = 0;
for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
wgt = GenerateNBodyEvent( Qvalue/MeV,
constantCrossSection, tempV, tempLen );
}
if(wgt>-.5)
{
theoreticalKinetic = 0.0;
for( i=0; i<tempLen; ++i )
{
pseudoParticle[6].Lorentz( *tempV[i], pseudoParticle[4] );
theoreticalKinetic += pseudoParticle[6].GetKineticEnergy()/MeV;
}
}
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
//
// Make sure, that the kinetic energies are correct
//
if( simulatedKinetic != 0.0 )
{
wgt = (theoreticalKinetic)/simulatedKinetic;
theoreticalKinetic = currentParticle.GetKineticEnergy()/MeV * wgt;
simulatedKinetic = theoreticalKinetic;
currentParticle.SetKineticEnergy( theoreticalKinetic*MeV );
pp = currentParticle.GetTotalMomentum()/MeV;
pp1 = currentParticle.GetMomentum().mag()/MeV;
if( pp1 < 1.0e-6*GeV )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi2 = twopi*G4UniformRand();
currentParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
pp*sintheta*std::sin(phi2)*MeV,
pp*costheta*MeV ) ;
}
else
{
currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
}
theoreticalKinetic = targetParticle.GetKineticEnergy()/MeV * wgt;
targetParticle.SetKineticEnergy( theoreticalKinetic*MeV );
simulatedKinetic += theoreticalKinetic;
pp = targetParticle.GetTotalMomentum()/MeV;
pp1 = targetParticle.GetMomentum().mag()/MeV;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
if( pp1 < 1.0e-6*GeV )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi2 = twopi*G4UniformRand();
targetParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
pp*sintheta*std::sin(phi2)*MeV,
pp*costheta*MeV ) ;
} else {
targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
}
for( i=0; i<vecLen; ++i )
{
theoreticalKinetic = vec[i]->GetKineticEnergy()/MeV * wgt;
simulatedKinetic += theoreticalKinetic;
vec[i]->SetKineticEnergy( theoreticalKinetic*MeV );
pp = vec[i]->GetTotalMomentum()/MeV;
pp1 = vec[i]->GetMomentum().mag()/MeV;
if( pp1 < 1.0e-6*GeV )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi2 = twopi*G4UniformRand();
vec[i]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
pp*sintheta*std::sin(phi2)*MeV,
pp*costheta*MeV ) ;
}
else
vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
}
}
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
Rotate( numberofFinalStateNucleons, pseudoParticle[3].GetMomentum(),
modifiedOriginal, originalIncident, targetNucleus,
currentParticle, targetParticle, vec, vecLen );
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
//
// add black track particles
// the total number of particles produced is restricted to 198
// this may have influence on very high energies
//
if( atomicWeight >= 1.5 )
{
// npnb is number of proton/neutron black track particles
// ndta is the number of deuterons, tritons, and alphas produced
// epnb is the kinetic energy available for proton/neutron black track particles
// edta is the kinetic energy available for deuteron/triton/alpha particles
//
G4int npnb = 0;
G4int ndta = 0;
G4double epnb, edta;
if (veryForward) {
epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
} else {
epnb = targetNucleus.GetPNBlackTrackEnergy();
edta = targetNucleus.GetDTABlackTrackEnergy();
}
const G4double pnCutOff = 0.001;
const G4double dtaCutOff = 0.001;
const G4double kineticMinimum = 1.e-6;
const G4double kineticFactor = -0.010;
G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules
const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
if( epnb >= pnCutOff )
{
npnb = Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
if( numberofFinalStateNucleons + npnb > atomicWeight )
npnb = G4int(atomicWeight+0.00001 - numberofFinalStateNucleons);
npnb = std::min( npnb, 127-vecLen );
}
if( edta >= dtaCutOff )
{
ndta = Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
ndta = std::min( ndta, 127-vecLen );
}
if (npnb == 0 && ndta == 0) npnb = 1;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum,
kineticFactor, modifiedOriginal,
PinNucleus, NinNucleus, targetNucleus,
vec, vecLen);
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
//if( centerofmassEnergy <= (4.0+G4UniformRand()) )
// MomentumCheck( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
//
// calculate time delay for nuclear reactions
//
if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
else
currentParticle.SetTOF( 1.0 );
return true;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
void G4ReactionDynamics::SuppressChargedPions(
G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
G4int &vecLen,
const G4ReactionProduct &modifiedOriginal,
G4ReactionProduct &currentParticle,
G4ReactionProduct &targetParticle,
const G4Nucleus &targetNucleus,
G4bool &incidentHasChanged,
G4bool &targetHasChanged )
{
// this code was originally in the fortran code TWOCLU
//
// suppress charged pions, for various reasons
//
G4double mOriginal = modifiedOriginal.GetMass()/GeV;
G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
G4double cmEnergy = std::sqrt( mOriginal*mOriginal + targetMass*targetMass +
2.0*targetMass*etOriginal );
G4double eAvailable = cmEnergy - mOriginal - targetMass;
for (G4int i = 0; i < vecLen; i++) eAvailable -= vec[i]->GetMass()/GeV;
const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/GeV;
G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
G4ParticleDefinition* aPiZero = G4PionZero::PionZero();
G4ParticleDefinition *aProton = G4Proton::Proton();
G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
G4double piMass = aPiPlus->GetPDGMass()/GeV;
G4double nucleonMass = aNeutron->GetPDGMass()/GeV;
const G4bool antiTest =
modifiedOriginal.GetDefinition() != G4AntiProton::AntiProton() &&
modifiedOriginal.GetDefinition() != G4AntiNeutron::AntiNeutron() &&
modifiedOriginal.GetDefinition() != G4AntiLambda::AntiLambda() &&
modifiedOriginal.GetDefinition() != G4AntiSigmaPlus::AntiSigmaPlus() &&
modifiedOriginal.GetDefinition() != G4AntiSigmaMinus::AntiSigmaMinus() &&
modifiedOriginal.GetDefinition() != G4AntiXiZero::AntiXiZero() &&
modifiedOriginal.GetDefinition() != G4AntiXiMinus::AntiXiMinus() &&
modifiedOriginal.GetDefinition() != G4AntiOmegaMinus::AntiOmegaMinus();
if( antiTest && (
currentParticle.GetDefinition() == aPiPlus ||
currentParticle.GetDefinition() == aPiZero ||
currentParticle.GetDefinition() == aPiMinus ) &&
( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
( G4UniformRand() <= atomicWeight/300.0 ) )
{
if (eAvailable > nucleonMass - piMass) {
if( G4UniformRand() > atomicNumber/atomicWeight )
currentParticle.SetDefinitionAndUpdateE( aNeutron );
else
currentParticle.SetDefinitionAndUpdateE( aProton );
incidentHasChanged = true;
}
}
if( antiTest && (
targetParticle.GetDefinition() == aPiPlus ||
targetParticle.GetDefinition() == aPiZero ||
targetParticle.GetDefinition() == aPiMinus ) &&
( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
( G4UniformRand() <= atomicWeight/300.0 ) )
{
if (eAvailable > nucleonMass - piMass) {
if( G4UniformRand() > atomicNumber/atomicWeight )
targetParticle.SetDefinitionAndUpdateE( aNeutron );
else
targetParticle.SetDefinitionAndUpdateE( aProton );
targetHasChanged = true;
}
}
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
for( G4int i=0; i<vecLen; ++i )
{
if( antiTest && (
vec[i]->GetDefinition() == aPiPlus ||
vec[i]->GetDefinition() == aPiZero ||
vec[i]->GetDefinition() == aPiMinus ) &&
( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
( G4UniformRand() <= atomicWeight/300.0 ) )
{
if (eAvailable > nucleonMass - piMass) {
if( G4UniformRand() > atomicNumber/atomicWeight )
vec[i]->SetDefinitionAndUpdateE( aNeutron );
else
vec[i]->SetDefinitionAndUpdateE( aProton );
}
}
}
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
G4bool G4ReactionDynamics::TwoCluster(
G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
G4int &vecLen,
G4ReactionProduct &modifiedOriginal, // Fermi motion & evap. effects included
const G4HadProjectile *originalIncident, // the original incident particle
G4ReactionProduct &currentParticle,
G4ReactionProduct &targetParticle,
const G4DynamicParticle* originalTarget,
const G4Nucleus &targetNucleus,
G4bool &incidentHasChanged,
G4bool &targetHasChanged,
G4bool leadFlag,
G4ReactionProduct &leadingStrangeParticle )
{
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
// derived from original FORTRAN code TWOCLU by H. Fesefeldt (11-Oct-1987)
//
// Generation of X- and PT- values for incident, target, and all secondary particles
//
// A simple two cluster model is used.
// This should be sufficient for low energy interactions.
//
// + debugging
// raise(SIGSEGV);
// - debugging
G4int i;
G4ParticleDefinition *aProton = G4Proton::Proton();
G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
G4bool veryForward = false;
const G4double protonMass = aProton->GetPDGMass()/MeV;
const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
targetMass*targetMass +
2.0*targetMass*etOriginal ); // GeV
G4double currentMass = currentParticle.GetMass()/GeV;
targetMass = targetParticle.GetMass()/GeV;
if( currentMass == 0.0 && targetMass == 0.0 )
{
G4double ek = currentParticle.GetKineticEnergy();
G4ThreeVector mom = currentParticle.GetMomentum();
currentParticle = *vec[0];
targetParticle = *vec[1];
for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
if(vecLen<2)
{
for(G4int j=0; j<vecLen; j++) delete vec[j];
vecLen = 0;
throw G4HadReentrentException(__FILE__, __LINE__,
"G4ReactionDynamics::TwoCluster: Negative number of particles");
}
delete vec[vecLen-1];
delete vec[vecLen-2];
vecLen -= 2;
currentMass = currentParticle.GetMass()/GeV;
targetMass = targetParticle.GetMass()/GeV;
incidentHasChanged = true;
targetHasChanged = true;
currentParticle.SetKineticEnergy( ek );
currentParticle.SetMomentum( mom );
veryForward = true;
}
const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
//
// particles have been distributed in forward and backward hemispheres
// in center of mass system of the hadron nucleon interaction
//
// incident is always in forward hemisphere
G4int forwardCount = 1; // number of particles in forward hemisphere
currentParticle.SetSide( 1 );
G4double forwardMass = currentParticle.GetMass()/GeV;
G4double cMass = forwardMass;
// target is always in backward hemisphere
G4int backwardCount = 1; // number of particles in backward hemisphere
G4int backwardNucleonCount = 1; // number of nucleons in backward hemisphere
targetParticle.SetSide( -1 );
G4double backwardMass = targetParticle.GetMass()/GeV;
G4double bMass = backwardMass;
for( i=0; i<vecLen; ++i )
{
if( vec[i]->GetSide() < 0 )vec[i]->SetSide( -1 ); // added by JLC, 2Jul97
// to take care of the case where vec has been preprocessed by GenerateXandPt
// and some of them have been set to -2 or -3
if( vec[i]->GetSide() == -1 )
{
++backwardCount;
backwardMass += vec[i]->GetMass()/GeV;
}
else
{
++forwardCount;
forwardMass += vec[i]->GetMass()/GeV;
}
}
//
// nucleons and some pions from intranuclear cascade
//
G4double term1 = std::log(centerofmassEnergy*centerofmassEnergy);
if(term1 < 0) term1 = 0.0001; // making sure xtarg<0;
const G4double afc = 0.312 + 0.2 * std::log(term1);
G4double xtarg;
if( centerofmassEnergy < 2.0+G4UniformRand() ) // added +2 below, JLC 4Jul97
xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount+vecLen+2)/2.0;
else
xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount);
if( xtarg <= 0.0 )xtarg = 0.01;
G4int nuclearExcitationCount = Poisson( xtarg );
if(atomicWeight<1.0001) nuclearExcitationCount = 0;
G4int extraNucleonCount = 0;
G4double extraMass = 0.0;
G4double extraNucleonMass = 0.0;
if( nuclearExcitationCount > 0 )
{
G4int momentumBin = std::min( 4, G4int(pOriginal/3.0) );
const G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4 };
//
// NOTE: in TWOCLU, these new particles were given negative codes
// here we use NewlyAdded = true instead
//
// G4ReactionProduct *pVec = new G4ReactionProduct [nuclearExcitationCount];
for( i=0; i<nuclearExcitationCount; ++i )
{
G4ReactionProduct* pVec = new G4ReactionProduct();
if( G4UniformRand() < nucsup[momentumBin] ) // add proton or neutron
{
if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
pVec->SetDefinition( aProton );
else
pVec->SetDefinition( aNeutron );
++backwardNucleonCount;
++extraNucleonCount;
extraNucleonMass += pVec->GetMass()/GeV;
}
else
{ // add a pion
G4double ran = G4UniformRand();
if( ran < 0.3181 )
pVec->SetDefinition( aPiPlus );
else if( ran < 0.6819 )
pVec->SetDefinition( aPiZero );
else
pVec->SetDefinition( aPiMinus );
}
pVec->SetSide( -2 ); // backside particle
extraMass += pVec->GetMass()/GeV;
pVec->SetNewlyAdded( true );
vec.SetElement( vecLen++, pVec );
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
}
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
G4double forwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +cMass - forwardMass;
G4double backwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +bMass - backwardMass;
G4double eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
G4bool secondaryDeleted;
G4double pMass;
while( eAvailable <= 0.0 ) // must eliminate a particle
{
secondaryDeleted = false;
for( i=(vecLen-1); i>=0; --i )
{
if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
{
pMass = vec[i]->GetMass()/GeV;
for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
--forwardCount;
forwardEnergy += pMass;
forwardMass -= pMass;
secondaryDeleted = true;
break;
}
else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
{
pMass = vec[i]->GetMass()/GeV;
for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
--backwardCount;
backwardEnergy += pMass;
backwardMass -= pMass;
secondaryDeleted = true;
break;
}
} // breaks go down to here
if( secondaryDeleted )
{
delete vec[vecLen-1];
--vecLen;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
else
{
if( vecLen == 0 )
{
return false; // all secondaries have been eliminated
}
if( targetParticle.GetSide() == -1 )
{
pMass = targetParticle.GetMass()/GeV;
targetParticle = *vec[0];
for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
--backwardCount;
backwardEnergy += pMass;
backwardMass -= pMass;
secondaryDeleted = true;
}
else if( targetParticle.GetSide() == 1 )
{
pMass = targetParticle.GetMass()/GeV;
targetParticle = *vec[0];
for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
--forwardCount;
forwardEnergy += pMass;
forwardMass -= pMass;
secondaryDeleted = true;
}
if( secondaryDeleted )
{
delete vec[vecLen-1];
--vecLen;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
else
{
if( currentParticle.GetSide() == -1 )
{
pMass = currentParticle.GetMass()/GeV;
currentParticle = *vec[0];
for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
--backwardCount;
backwardEnergy += pMass;
backwardMass -= pMass;
secondaryDeleted = true;
}
else if( currentParticle.GetSide() == 1 )
{
pMass = currentParticle.GetMass()/GeV;
currentParticle = *vec[0];
for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
--forwardCount;
forwardEnergy += pMass;
forwardMass -= pMass;
secondaryDeleted = true;
}
if( secondaryDeleted )
{
delete vec[vecLen-1];
--vecLen;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
else break;
}
}
eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
}
//
// This is the start of the TwoCluster function
// Choose masses for the 3 clusters:
// forward cluster
// backward meson cluster
// backward nucleon cluster
//
G4double rmc = 0.0, rmd = 0.0;
const G4double cpar[] = { 0.6, 0.6, 0.35, 0.15, 0.10 };
const G4double gpar[] = { 2.6, 2.6, 1.8, 1.30, 1.20 };
if (forwardCount <= 0 || backwardCount <= 0) return false; // array bounds protection
if (forwardCount == 1) rmc = forwardMass;
else
{
G4int ntc = std::max(1, std::min(5,forwardCount))-1; // check if offset by 1 @@
rmc = forwardMass + std::pow(-std::log(1.0-G4UniformRand()),cpar[ntc-1])/gpar[ntc-1];
}
if (backwardCount == 1) rmd = backwardMass;
else
{
G4int ntc = std::max(1, std::min(5,backwardCount)); // check, if offfset by 1 @@
rmd = backwardMass + std::pow(-std::log(1.0-G4UniformRand()),cpar[ntc-1])/gpar[ntc-1];
}
while( rmc+rmd > centerofmassEnergy )
{
if( (rmc <= forwardMass) && (rmd <= backwardMass) )
{
G4double temp = 0.999*centerofmassEnergy/(rmc+rmd);
rmc *= temp;
rmd *= temp;
}
else
{
rmc = 0.1*forwardMass + 0.9*rmc;
rmd = 0.1*backwardMass + 0.9*rmd;
}
}
G4ReactionProduct pseudoParticle[8];
for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
pseudoParticle[1].SetMass( mOriginal*GeV );
pseudoParticle[1].SetTotalEnergy( etOriginal*GeV );
pseudoParticle[1].SetMomentum( 0.0, 0.0, pOriginal*GeV );
pseudoParticle[2].SetMass( protonMass*MeV );
pseudoParticle[2].SetTotalEnergy( protonMass*MeV );
pseudoParticle[2].SetMomentum( 0.0, 0.0, 0.0 );
//
// transform into centre of mass system
//
pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[0] );
pseudoParticle[2].Lorentz( pseudoParticle[2], pseudoParticle[0] );
const G4double pfMin = 0.0001;
G4double pf = (centerofmassEnergy*centerofmassEnergy+rmd*rmd-rmc*rmc);
pf *= pf;
pf -= 4*centerofmassEnergy*centerofmassEnergy*rmd*rmd;
pf = std::sqrt( std::max(pf,pfMin) )/(2.0*centerofmassEnergy);
//
// set final state masses and energies in centre of mass system
//
pseudoParticle[3].SetMass( rmc*GeV );
pseudoParticle[3].SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*GeV );
pseudoParticle[4].SetMass( rmd*GeV );
pseudoParticle[4].SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*GeV );
//
// set |T| and |TMIN|
//
const G4double bMin = 0.01;
const G4double b1 = 4.0;
const G4double b2 = 1.6;
G4double pin = pseudoParticle[1].GetMomentum().mag()/GeV;
G4double dtb = 4.0*pin*pf*std::max( bMin, b1+b2*std::log(pOriginal) );
G4double factor = 1.0 - std::exp(-dtb);
G4double costheta = 1.0 + 2.0*std::log(1.0 - G4UniformRand()*factor) / dtb;
costheta = std::max(-1.0, std::min(1.0, costheta));
G4double sintheta = std::sqrt(1.0-costheta*costheta);
G4double phi = G4UniformRand() * twopi;
//
// calculate final state momenta in centre of mass system
//
pseudoParticle[3].SetMomentum( pf*sintheta*std::cos(phi)*GeV,
pf*sintheta*std::sin(phi)*GeV,
pf*costheta*GeV );
pseudoParticle[4].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
//
// simulate backward nucleon cluster in lab. system and transform in cms
//
G4double pp, pp1;
if( nuclearExcitationCount > 0 )
{
const G4double ga = 1.2;
G4double ekit1 = 0.04;
G4double ekit2 = 0.6;
if( ekOriginal <= 5.0 )
{
ekit1 *= ekOriginal*ekOriginal/25.0;
ekit2 *= ekOriginal*ekOriginal/25.0;
}
const G4double a = (1.0-ga)/(std::pow(ekit2,(1.0-ga)) - std::pow(ekit1,(1.0-ga)));
for( i=0; i<vecLen; ++i )
{
if( vec[i]->GetSide() == -2 )
{
G4double kineticE =
std::pow( (G4UniformRand()*(1.0-ga)/a+std::pow(ekit1,(1.0-ga))), (1.0/(1.0-ga)) );
vec[i]->SetKineticEnergy( kineticE*GeV );
G4double vMass = vec[i]->GetMass()/MeV;
G4double totalE = kineticE*GeV + vMass;
pp = std::sqrt( std::abs(totalE*totalE-vMass*vMass) );
G4double cost = std::min( 1.0, std::max( -1.0, std::log(2.23*G4UniformRand()+0.383)/0.96 ) );
G4double sint = std::sqrt( std::max( 0.0, (1.0-cost*cost) ) );
phi = twopi*G4UniformRand();
vec[i]->SetMomentum( pp*sint*std::sin(phi)*MeV,
pp*sint*std::cos(phi)*MeV,
pp*cost*MeV );
vec[i]->Lorentz( *vec[i], pseudoParticle[0] );
}
}
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
//
// fragmentation of forward cluster and backward meson cluster
//
currentParticle.SetMomentum( pseudoParticle[3].GetMomentum() );
currentParticle.SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
targetParticle.SetMomentum( pseudoParticle[4].GetMomentum() );
targetParticle.SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
pseudoParticle[5].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
pseudoParticle[5].SetMass( pseudoParticle[3].GetMass() );
pseudoParticle[5].SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
pseudoParticle[6].SetMomentum( pseudoParticle[4].GetMomentum() * (-1.0) );
pseudoParticle[6].SetMass( pseudoParticle[4].GetMass() );
pseudoParticle[6].SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
G4double wgt;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
if( forwardCount > 1 ) // tempV will contain the forward particles
{
G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
tempV.Initialize( forwardCount );
G4bool constantCrossSection = true;
G4int tempLen = 0;
if( currentParticle.GetSide() == 1 )
tempV.SetElement( tempLen++, &currentParticle );
if( targetParticle.GetSide() == 1 )
tempV.SetElement( tempLen++, &targetParticle );
for( i=0; i<vecLen; ++i )
{
if( vec[i]->GetSide() == 1 )
{
if( tempLen < 18 )
tempV.SetElement( tempLen++, vec[i] );
else
{
vec[i]->SetSide( -1 );
continue;
}
}
}
if( tempLen >= 2 )
{
wgt = GenerateNBodyEvent( pseudoParticle[3].GetMass()/MeV,
constantCrossSection, tempV, tempLen );
if( currentParticle.GetSide() == 1 )
currentParticle.Lorentz( currentParticle, pseudoParticle[5] );
if( targetParticle.GetSide() == 1 )
targetParticle.Lorentz( targetParticle, pseudoParticle[5] );
for( i=0; i<vecLen; ++i )
{
if( vec[i]->GetSide() == 1 )vec[i]->Lorentz( *vec[i], pseudoParticle[5] );
}
}
}
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
if( backwardCount > 1 ) // tempV will contain the backward particles,
{ // but not those created from the intranuclear cascade
G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
tempV.Initialize( backwardCount );
G4bool constantCrossSection = true;
G4int tempLen = 0;
if( currentParticle.GetSide() == -1 )
tempV.SetElement( tempLen++, &currentParticle );
if( targetParticle.GetSide() == -1 )
tempV.SetElement( tempLen++, &targetParticle );
for( i=0; i<vecLen; ++i )
{
if( vec[i]->GetSide() == -1 )
{
if( tempLen < 18 )
tempV.SetElement( tempLen++, vec[i] );
else
{
vec[i]->SetSide( -2 );
vec[i]->SetKineticEnergy( 0.0 );
vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
continue;
}
}
}
if( tempLen >= 2 )
{
wgt = GenerateNBodyEvent( pseudoParticle[4].GetMass()/MeV,
constantCrossSection, tempV, tempLen );
if( currentParticle.GetSide() == -1 )
currentParticle.Lorentz( currentParticle, pseudoParticle[6] );
if( targetParticle.GetSide() == -1 )
targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
for( i=0; i<vecLen; ++i )
{
if( vec[i]->GetSide() == -1 )vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
}
}
}
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
//
// Lorentz transformation in lab system
//
currentParticle.Lorentz( currentParticle, pseudoParticle[2] );
targetParticle.Lorentz( targetParticle, pseudoParticle[2] );
for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[2] );
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
//
// sometimes the leading strange particle is lost, set it back
//
G4bool dum = true;
if( leadFlag )
{
// leadFlag will be true
// iff original particle is at least as heavy as K+ and not a proton or
// neutron AND if incident particle is at least as heavy as K+ and it is
// not a proton or neutron leadFlag is set to the incident particle
// or
// target particle is at least as heavy as K+ and it is not a proton or
// neutron leadFlag is set to the target particle
//
if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
dum = false;
else if( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
dum = false;
else
{
for( i=0; i<vecLen; ++i )
{
if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
{
dum = false;
break;
}
}
}
if( dum )
{
G4double leadMass = leadingStrangeParticle.GetMass()/MeV;
G4double ekin;
if( ((leadMass < protonMass) && (targetParticle.GetMass()/MeV < protonMass)) ||
((leadMass >= protonMass) && (targetParticle.GetMass()/MeV >= protonMass)) )
{
ekin = targetParticle.GetKineticEnergy()/GeV;
pp1 = targetParticle.GetMomentum().mag()/MeV; // old momentum
targetParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
targetParticle.SetKineticEnergy( ekin*GeV );
pp = targetParticle.GetTotalMomentum()/MeV; // new momentum
if( pp1 < 1.0e-3 )
{
G4double costheta2 = 2.*G4UniformRand() - 1.;
G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
G4double phi2 = twopi*G4UniformRand();
targetParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
pp*sintheta2*std::sin(phi2)*MeV,
pp*costheta2*MeV ) ;
}
else
targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
targetHasChanged = true;
}
else
{
ekin = currentParticle.GetKineticEnergy()/GeV;
pp1 = currentParticle.GetMomentum().mag()/MeV;
currentParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
currentParticle.SetKineticEnergy( ekin*GeV );
pp = currentParticle.GetTotalMomentum()/MeV;
if( pp1 < 1.0e-3 )
{
G4double costheta2 = 2.*G4UniformRand() - 1.;
G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
G4double phi2 = twopi*G4UniformRand();
currentParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
pp*sintheta2*std::sin(phi2)*MeV,
pp*costheta2*MeV ) ;
}
else
currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
incidentHasChanged = true;
}
}
} // end of if( leadFlag )
// Get number of final state nucleons and nucleons remaining in
// target nucleus
std::pair<G4int, G4int> finalStateNucleons =
GetFinalStateNucleons(originalTarget, vec, vecLen);
G4int protonsInFinalState = finalStateNucleons.first;
G4int neutronsInFinalState = finalStateNucleons.second;
G4int numberofFinalStateNucleons =
protonsInFinalState + neutronsInFinalState;
if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
originalIncident->GetDefinition()->GetPDGMass() <
G4Lambda::Lambda()->GetPDGMass())
numberofFinalStateNucleons++;
numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
G4int PinNucleus = std::max(0,
targetNucleus.GetZ_asInt() - protonsInFinalState);
G4int NinNucleus = std::max(0,
targetNucleus.GetN_asInt() - neutronsInFinalState);
//
// for various reasons, the energy balance is not sufficient,
// check that, energy balance, angle of final system, etc.
//
pseudoParticle[4].SetMass( mOriginal*GeV );
pseudoParticle[4].SetTotalEnergy( etOriginal*GeV );
pseudoParticle[4].SetMomentum( 0.0, 0.0, pOriginal*GeV );
G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
G4int diff = 0;
if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
if(numberofFinalStateNucleons == 1) diff = 0;
pseudoParticle[5].SetMomentum( 0.0, 0.0, 0.0 );
pseudoParticle[5].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV);
pseudoParticle[5].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV);
G4double theoreticalKinetic =
pseudoParticle[4].GetTotalEnergy()/GeV + pseudoParticle[5].GetTotalEnergy()/GeV;
pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[6] );
pseudoParticle[5].Lorentz( pseudoParticle[5], pseudoParticle[6] );
if( vecLen < 16 )
{
G4ReactionProduct tempR[130];
tempR[0] = currentParticle;
tempR[1] = targetParticle;
for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
tempV.Initialize( vecLen+2 );
G4bool constantCrossSection = true;
G4int tempLen = 0;
for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
if( tempLen >= 2 )
{
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
wgt = GenerateNBodyEvent( pseudoParticle[4].GetTotalEnergy()/MeV +
pseudoParticle[5].GetTotalEnergy()/MeV,
constantCrossSection, tempV, tempLen );
if (wgt == -1) {
G4double Qvalue = 0;
for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
wgt = GenerateNBodyEvent( Qvalue/MeV,
constantCrossSection, tempV, tempLen );
}
theoreticalKinetic = 0.0;
for( i=0; i<vecLen+2; ++i )
{
pseudoParticle[7].SetMomentum( tempV[i]->GetMomentum() );
pseudoParticle[7].SetMass( tempV[i]->GetMass() );
pseudoParticle[7].SetTotalEnergy( tempV[i]->GetTotalEnergy() );
pseudoParticle[7].Lorentz( pseudoParticle[7], pseudoParticle[5] );
theoreticalKinetic += pseudoParticle[7].GetKineticEnergy()/GeV;
}
}
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
else
{
theoreticalKinetic -=
( currentParticle.GetMass()/GeV + targetParticle.GetMass()/GeV );
for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/GeV;
}
G4double simulatedKinetic =
currentParticle.GetKineticEnergy()/GeV + targetParticle.GetKineticEnergy()/GeV;
for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/GeV;
//
// make sure that kinetic energies are correct
// the backward nucleon cluster is not produced within proper kinematics!!!
//
if( simulatedKinetic != 0.0 )
{
wgt = (theoreticalKinetic)/simulatedKinetic;
currentParticle.SetKineticEnergy( wgt*currentParticle.GetKineticEnergy() );
pp = currentParticle.GetTotalMomentum()/MeV;
pp1 = currentParticle.GetMomentum().mag()/MeV;
if( pp1 < 0.001*MeV )
{
G4double costheta2 = 2.*G4UniformRand() - 1.;
G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
G4double phi2 = twopi*G4UniformRand();
currentParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
pp*sintheta2*std::sin(phi2)*MeV,
pp*costheta2*MeV ) ;
}
else
currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
targetParticle.SetKineticEnergy( wgt*targetParticle.GetKineticEnergy() );
pp = targetParticle.GetTotalMomentum()/MeV;
pp1 = targetParticle.GetMomentum().mag()/MeV;
if( pp1 < 0.001*MeV )
{
G4double costheta2 = 2.*G4UniformRand() - 1.;
G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
G4double phi2 = twopi*G4UniformRand();
targetParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
pp*sintheta2*std::sin(phi2)*MeV,
pp*costheta2*MeV ) ;
}
else
targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
for( i=0; i<vecLen; ++i )
{
vec[i]->SetKineticEnergy( wgt*vec[i]->GetKineticEnergy() );
pp = vec[i]->GetTotalMomentum()/MeV;
pp1 = vec[i]->GetMomentum().mag()/MeV;
if( pp1 < 0.001 )
{
G4double costheta2 = 2.*G4UniformRand() - 1.;
G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
G4double phi2 = twopi*G4UniformRand();
vec[i]->SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
pp*sintheta2*std::sin(phi2)*MeV,
pp*costheta2*MeV ) ;
}
else
vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
}
}
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
Rotate( numberofFinalStateNucleons, pseudoParticle[4].GetMomentum(),
modifiedOriginal, originalIncident, targetNucleus,
currentParticle, targetParticle, vec, vecLen );
//
// add black track particles
// the total number of particles produced is restricted to 198
// this may have influence on very high energies
//
if( atomicWeight >= 1.5 )
{
// npnb is number of proton/neutron black track particles
// ndta is the number of deuterons, tritons, and alphas produced
// epnb is the kinetic energy available for proton/neutron black track
// particles
// edta is the kinetic energy available for deuteron/triton/alpha
// particles
G4int npnb = 0;
G4int ndta = 0;
G4double epnb, edta;
if (veryForward) {
epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
} else {
epnb = targetNucleus.GetPNBlackTrackEnergy();
edta = targetNucleus.GetDTABlackTrackEnergy();
}
const G4double pnCutOff = 0.001; // GeV
const G4double dtaCutOff = 0.001; // GeV
const G4double kineticMinimum = 1.e-6;
const G4double kineticFactor = -0.005;
G4double sprob = 0.0; // sprob = probability of self-absorption in
// heavy molecules
const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
if( epnb >= pnCutOff )
{
npnb = Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
if( numberofFinalStateNucleons + npnb > atomicWeight )
npnb = G4int(atomicWeight - numberofFinalStateNucleons);
npnb = std::min( npnb, 127-vecLen );
}
if( edta >= dtaCutOff )
{
ndta = Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
ndta = std::min( ndta, 127-vecLen );
}
if (npnb == 0 && ndta == 0) npnb = 1;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum,
kineticFactor, modifiedOriginal,
PinNucleus, NinNucleus, targetNucleus,
vec, vecLen );
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
//if( centerofmassEnergy <= (4.0+G4UniformRand()) )
// MomentumCheck( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
//
// calculate time delay for nuclear reactions
//
if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
else
currentParticle.SetTOF( 1.0 );
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
return true;
}
void G4ReactionDynamics::TwoBody(
G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
G4int &vecLen,
G4ReactionProduct &modifiedOriginal,
const G4DynamicParticle* originalTarget,
G4ReactionProduct &currentParticle,
G4ReactionProduct &targetParticle,
const G4Nucleus &targetNucleus,
G4bool &/* targetHasChanged*/ )
{
//
// derived from original FORTRAN code TWOB by H. Fesefeldt (15-Sep-1987)
//
// Generation of momenta for elastic and quasi-elastic 2 body reactions
//
// The simple formula ds/d|t| = s0* std::exp(-b*|t|) is used.
// The b values are parametrizations from experimental data.
// Not available values are taken from those of similar reactions.
//
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
static const G4double expxu = 82.; // upper bound for arg. of exp
static const G4double expxl = -expxu; // lower bound for arg. of exp
const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
G4double currentMass = currentParticle.GetMass()/GeV;
G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
targetMass = targetParticle.GetMass()/GeV;
const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
G4double etCurrent = currentParticle.GetTotalEnergy()/GeV;
G4double pCurrent = currentParticle.GetTotalMomentum()/GeV;
G4double cmEnergy = std::sqrt( currentMass*currentMass +
targetMass*targetMass +
2.0*targetMass*etCurrent ); // in GeV
//if( (pOriginal < 0.1) ||
// (centerofmassEnergy < 0.01) ) // 2-body scattering not possible
// Continue with original particle, but spend the nuclear evaporation energy
// targetParticle.SetMass( 0.0 ); // flag that the target doesn't exist
//else // Two-body scattering is possible
if( (pCurrent < 0.1) || (cmEnergy < 0.01) ) // 2-body scattering not possible
{
targetParticle.SetMass( 0.0 ); // flag that the target particle doesn't exist
}
else
{
// moved this if-block to a later stage, i.e. to the assignment of the scattering angle
// @@@@@ double-check.
// if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
// targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
// if( G4UniformRand() < 0.5 )
// targetParticle.SetDefinitionAndUpdateE( aNeutron );
// else
// targetParticle.SetDefinitionAndUpdateE( aProton );
// targetHasChanged = true;
// targetMass = targetParticle.GetMass()/GeV;
// }
//
// Set masses and momenta for final state particles
//
G4double pf = cmEnergy*cmEnergy + targetMass*targetMass - currentMass*currentMass;
pf = pf*pf - 4*cmEnergy*cmEnergy*targetMass*targetMass;
if( pf < 0.001 )
{
for(G4int i=0; i<vecLen; i++) delete vec[i];
vecLen = 0;
throw G4HadronicException(__FILE__, __LINE__, "G4ReactionDynamics::TwoBody: pf is too small ");
}
pf = std::sqrt( pf ) / ( 2.0*cmEnergy );
//
// Set beam and target in centre of mass system
//
G4ReactionProduct pseudoParticle[3];
if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
pseudoParticle[0].SetMass( targetMass*GeV );
pseudoParticle[0].SetTotalEnergy( etOriginal*GeV );
pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
pseudoParticle[1].SetMass( mOriginal*GeV );
pseudoParticle[1].SetKineticEnergy( 0.0 );
} else {
pseudoParticle[0].SetMass( currentMass*GeV );
pseudoParticle[0].SetTotalEnergy( etCurrent*GeV );
pseudoParticle[0].SetMomentum( 0.0, 0.0, pCurrent*GeV );
pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
pseudoParticle[1].SetMass( targetMass*GeV );
pseudoParticle[1].SetKineticEnergy( 0.0 );
}
//
// Transform into centre of mass system
//
pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
//
// Set final state masses and energies in centre of mass system
//
currentParticle.SetTotalEnergy( std::sqrt(pf*pf+currentMass*currentMass)*GeV );
targetParticle.SetTotalEnergy( std::sqrt(pf*pf+targetMass*targetMass)*GeV );
//
// Set |t| and |tmin|
//
const G4double cb = 0.01;
const G4double b1 = 4.225;
const G4double b2 = 1.795;
//
// Calculate slope b for elastic scattering on proton/neutron
//
G4double b = std::max( cb, b1+b2*std::log(pOriginal) );
G4double btrang = b * 4.0 * pf * pseudoParticle[0].GetMomentum().mag()/GeV;
G4double exindt = -1.0;
exindt += std::exp(std::max(-btrang,expxl));
//
// Calculate sqr(std::sin(teta/2.) and std::cos(teta), set azimuth angle phi
//
G4double ctet = 1.0 + 2*std::log( 1.0+G4UniformRand()*exindt ) / btrang;
if( std::fabs(ctet) > 1.0 )ctet > 0.0 ? ctet = 1.0 : ctet = -1.0;
G4double stet = std::sqrt( (1.0-ctet)*(1.0+ctet) );
G4double phi = twopi * G4UniformRand();
//
// Calculate final state momenta in centre of mass system
//
if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
currentParticle.SetMomentum( -pf*stet*std::sin(phi)*GeV,
-pf*stet*std::cos(phi)*GeV,
-pf*ctet*GeV );
} else {
currentParticle.SetMomentum( pf*stet*std::sin(phi)*GeV,
pf*stet*std::cos(phi)*GeV,
pf*ctet*GeV );
}
targetParticle.SetMomentum( currentParticle.GetMomentum() * (-1.0) );
//
// Transform into lab system
//
currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
G4double pp, pp1, ekin;
if( atomicWeight >= 1.5 )
{
const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
pp1 = currentParticle.GetMomentum().mag()/MeV;
if( pp1 >= 1.0 )
{
ekin = currentParticle.GetKineticEnergy()/MeV - cfa*(1.0+0.5*normal())*GeV;
ekin = std::max( 0.0001*GeV, ekin );
currentParticle.SetKineticEnergy( ekin*MeV );
pp = currentParticle.GetTotalMomentum()/MeV;
currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
}
pp1 = targetParticle.GetMomentum().mag()/MeV;
if( pp1 >= 1.0 )
{
ekin = targetParticle.GetKineticEnergy()/MeV - cfa*(1.0+normal()/2.)*GeV;
ekin = std::max( 0.0001*GeV, ekin );
targetParticle.SetKineticEnergy( ekin*MeV );
pp = targetParticle.GetTotalMomentum()/MeV;
targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
}
}
}
// Get number of final state nucleons and nucleons remaining in
// target nucleus
std::pair<G4int, G4int> finalStateNucleons =
GetFinalStateNucleons(originalTarget, vec, vecLen);
G4int protonsInFinalState = finalStateNucleons.first;
G4int neutronsInFinalState = finalStateNucleons.second;
G4int PinNucleus = std::max(0,
targetNucleus.GetZ_asInt() - protonsInFinalState);
G4int NinNucleus = std::max(0,
targetNucleus.GetN_asInt() - neutronsInFinalState);
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
if( atomicWeight >= 1.5 )
{
// Add black track particles
// npnb is number of proton/neutron black track particles
// ndta is the number of deuterons, tritons, and alphas produced
// epnb is the kinetic energy available for proton/neutron black track particles
// edta is the kinetic energy available for deuteron/triton/alpha particles
//
G4double epnb, edta;
G4int npnb=0, ndta=0;
epnb = targetNucleus.GetPNBlackTrackEnergy(); // was enp1 in fortran code
edta = targetNucleus.GetDTABlackTrackEnergy(); // was enp3 in fortran code
const G4double pnCutOff = 0.0001; // GeV
const G4double dtaCutOff = 0.0001; // GeV
const G4double kineticMinimum = 0.0001;
const G4double kineticFactor = -0.010;
G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules
if( epnb >= pnCutOff )
{
npnb = Poisson( epnb/0.02 );
if( npnb > atomicWeight )npnb = G4int(atomicWeight);
if( (epnb > pnCutOff) && (npnb <= 0) )npnb = 1;
npnb = std::min( npnb, 127-vecLen );
}
if( edta >= dtaCutOff )
{
ndta = G4int(2.0 * std::log(atomicWeight));
ndta = std::min( ndta, 127-vecLen );
}
if (npnb == 0 && ndta == 0) npnb = 1;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum,
kineticFactor, modifiedOriginal,
PinNucleus, NinNucleus, targetNucleus,
vec, vecLen);
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
//
// calculate time delay for nuclear reactions
//
if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
else
currentParticle.SetTOF( 1.0 );
return;
}
G4double G4ReactionDynamics::GenerateNBodyEvent(
const G4double totalEnergy, // MeV
const G4bool constantCrossSection,
G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
G4int &vecLen )
{
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
// derived from original FORTRAN code PHASP by H. Fesefeldt (02-Dec-1986)
// Returns the weight of the event
//
G4int i;
const G4double expxu = 82.; // upper bound for arg. of exp
const G4double expxl = -expxu; // lower bound for arg. of exp
if( vecLen < 2 )
{
G4cerr << "*** Error in G4ReactionDynamics::GenerateNBodyEvent" << G4endl;
G4cerr << " number of particles < 2" << G4endl;
G4cerr << "totalEnergy = " << totalEnergy << "MeV, vecLen = " << vecLen << G4endl;
return -1.0;
}
G4double mass[18]; // mass of each particle
G4double energy[18]; // total energy of each particle
G4double pcm[3][18]; // pcm is an array with 3 rows and vecLen columns
G4double totalMass = 0.0;
G4double extraMass = 0;
G4double sm[18];
for( i=0; i<vecLen; ++i )
{
mass[i] = vec[i]->GetMass()/GeV;
if(vec[i]->GetSide() == -2) extraMass+=vec[i]->GetMass()/GeV;
vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
pcm[0][i] = 0.0; // x-momentum of i-th particle
pcm[1][i] = 0.0; // y-momentum of i-th particle
pcm[2][i] = 0.0; // z-momentum of i-th particle
energy[i] = mass[i]; // total energy of i-th particle
totalMass += mass[i];
sm[i] = totalMass;
}
G4double totalE = totalEnergy/GeV;
if( totalMass > totalE )
{
//G4cerr << "*** Error in G4ReactionDynamics::GenerateNBodyEvent" << G4endl;
//G4cerr << " total mass (" << totalMass*GeV << "MeV) > total energy ("
// << totalEnergy << "MeV)" << G4endl;
totalE = totalMass;
return -1.0;
}
G4double kineticEnergy = totalE - totalMass;
G4double emm[18];
//G4double *emm = new G4double [vecLen];
emm[0] = mass[0];
emm[vecLen-1] = totalE;
if( vecLen > 2 ) // the random numbers are sorted
{
G4double ran[18];
for( i=0; i<vecLen; ++i )ran[i] = G4UniformRand();
for( i=0; i<vecLen-2; ++i )
{
for( G4int j=vecLen-2; j>i; --j )
{
if( ran[i] > ran[j] )
{
G4double temp = ran[i];
ran[i] = ran[j];
ran[j] = temp;
}
}
}
for( i=1; i<vecLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
}
// Weight is the sum of logarithms of terms instead of the product of terms
G4bool lzero = true;
G4double wtmax = 0.0;
if( constantCrossSection ) // this is KGENEV=1 in PHASP
{
G4double emmax = kineticEnergy + mass[0];
G4double emmin = 0.0;
for( i=1; i<vecLen; ++i )
{
emmin += mass[i-1];
emmax += mass[i];
G4double wtfc = 0.0;
if( emmax*emmax > 0.0 )
{
G4double arg = emmax*emmax
+ (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
- 2.0*(emmin*emmin+mass[i]*mass[i]);
if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
}
if( wtfc == 0.0 )
{
lzero = false;
break;
}
wtmax += std::log( wtfc );
}
if( lzero )
wtmax = -wtmax;
else
wtmax = expxu;
}
else
{
// ffq(n) = pi*(2*pi)^(n-2)/(n-2)!
const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
256.3704, 268.4705, 240.9780, 189.2637,
132.1308, 83.0202, 47.4210, 24.8295,
12.0006, 5.3858, 2.2560, 0.8859 };
wtmax = std::log( std::pow( kineticEnergy, vecLen-2 ) * ffq[vecLen-1] / totalE );
}
lzero = true;
G4double pd[50];
//G4double *pd = new G4double [vecLen-1];
for( i=0; i<vecLen-1; ++i )
{
pd[i] = 0.0;
if( emm[i+1]*emm[i+1] > 0.0 )
{
G4double arg = emm[i+1]*emm[i+1]
+ (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
/(emm[i+1]*emm[i+1])
- 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
}
if( pd[i] <= 0.0 ) // changed from == on 02 April 98
lzero = false;
else
wtmax += std::log( pd[i] );
}
G4double weight = 0.0; // weight is returned by GenerateNBodyEvent
if( lzero )weight = std::exp( std::max(std::min(wtmax,expxu),expxl) );
G4double bang, cb, sb, s0, s1, s2, c, ss, esys, a, b, gama, beta;
pcm[0][0] = 0.0;
pcm[1][0] = pd[0];
pcm[2][0] = 0.0;
for( i=1; i<vecLen; ++i )
{
pcm[0][i] = 0.0;
pcm[1][i] = -pd[i-1];
pcm[2][i] = 0.0;
bang = twopi*G4UniformRand();
cb = std::cos(bang);
sb = std::sin(bang);
c = 2.0*G4UniformRand() - 1.0;
ss = std::sqrt( std::fabs( 1.0-c*c ) );
if( i < vecLen-1 )
{
esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
beta = pd[i]/esys;
gama = esys/emm[i];
for( G4int j=0; j<=i; ++j )
{
s0 = pcm[0][j];
s1 = pcm[1][j];
s2 = pcm[2][j];
energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
a = s0*c - s1*ss; // rotation
pcm[1][j] = s0*ss + s1*c;
b = pcm[2][j];
pcm[0][j] = a*cb - b*sb;
pcm[2][j] = a*sb + b*cb;
pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
}
}
else
{
for( G4int j=0; j<=i; ++j )
{
s0 = pcm[0][j];
s1 = pcm[1][j];
s2 = pcm[2][j];
energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
a = s0*c - s1*ss; // rotation
pcm[1][j] = s0*ss + s1*c;
b = pcm[2][j];
pcm[0][j] = a*cb - b*sb;
pcm[2][j] = a*sb + b*cb;
}
}
}
for( i=0; i<vecLen; ++i )
{
vec[i]->SetMomentum( pcm[0][i]*GeV, pcm[1][i]*GeV, pcm[2][i]*GeV );
vec[i]->SetTotalEnergy( energy[i]*GeV );
}
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
return weight;
}
G4double
G4ReactionDynamics::normal()
{
G4double ran = -6.0;
for( G4int i=0; i<12; ++i )ran += G4UniformRand();
return ran;
}
G4int
G4ReactionDynamics::Poisson( G4double x ) // generation of poisson distribution
{
G4int iran;
G4double ran;
if( x > 9.9 ) // use normal distribution with sigma^2 = <x>
iran = static_cast<G4int>(std::max( 0.0, x+normal()*std::sqrt(x) ) );
else {
G4int mn = G4int(5.0*x);
if( mn <= 0 ) // for very small x try iran=1,2,3
{
G4double p1 = x*std::exp(-x);
G4double p2 = x*p1/2.0;
G4double p3 = x*p2/3.0;
ran = G4UniformRand();
if( ran < p3 )
iran = 3;
else if( ran < p2 ) // this is original Geisha, it should be ran < p2+p3
iran = 2;
else if( ran < p1 ) // should be ran < p1+p2+p3
iran = 1;
else
iran = 0;
}
else
{
iran = 0;
G4double r = std::exp(-x);
ran = G4UniformRand();
if( ran > r )
{
G4double rrr;
G4double rr = r;
for( G4int i=1; i<=mn; ++i )
{
iran++;
if( i > 5 ) // Stirling's formula for large numbers
rrr = std::exp(i*std::log(x)-(i+0.5)*std::log((G4double)i)+i-0.9189385);
else
rrr = std::pow(x,i)/Factorial(i);
rr += r*rrr;
if( ran <= rr )break;
}
}
}
}
return iran;
}
G4int
G4ReactionDynamics::Factorial( G4int n )
{ // calculates factorial( n ) = n*(n-1)*(n-2)*...*1
G4int mn = std::min(n,10);
G4int result = 1;
if( mn <= 1 )return result;
for( G4int i=2; i<=mn; ++i )result *= i;
return result;
}
void G4ReactionDynamics::Defs1(
const G4ReactionProduct &modifiedOriginal,
G4ReactionProduct &currentParticle,
G4ReactionProduct &targetParticle,
G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
G4int &vecLen )
{
const G4double pjx = modifiedOriginal.GetMomentum().x()/MeV;
const G4double pjy = modifiedOriginal.GetMomentum().y()/MeV;
const G4double pjz = modifiedOriginal.GetMomentum().z()/MeV;
const G4double p = modifiedOriginal.GetMomentum().mag()/MeV;
if( pjx*pjx+pjy*pjy > 0.0 )
{
G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
cost = pjz/p;
sint = 0.5 * ( std::sqrt(std::abs((1.0-cost)*(1.0+cost))) + std::sqrt(pjx*pjx+pjy*pjy)/p );
if( pjy < 0.0 )
ph = 3*halfpi;
else
ph = halfpi;
if( std::abs( pjx ) > 0.001*MeV )ph = std::atan2(pjy,pjx);
cosp = std::cos(ph);
sinp = std::sin(ph);
pix = currentParticle.GetMomentum().x()/MeV;
piy = currentParticle.GetMomentum().y()/MeV;
piz = currentParticle.GetMomentum().z()/MeV;
currentParticle.SetMomentum( cost*cosp*pix*MeV - sinp*piy+sint*cosp*piz*MeV,
cost*sinp*pix*MeV + cosp*piy+sint*sinp*piz*MeV,
-sint*pix*MeV + cost*piz*MeV );
pix = targetParticle.GetMomentum().x()/MeV;
piy = targetParticle.GetMomentum().y()/MeV;
piz = targetParticle.GetMomentum().z()/MeV;
targetParticle.SetMomentum( cost*cosp*pix*MeV - sinp*piy+sint*cosp*piz*MeV,
cost*sinp*pix*MeV + cosp*piy+sint*sinp*piz*MeV,
-sint*pix*MeV + cost*piz*MeV );
for( G4int i=0; i<vecLen; ++i )
{
pix = vec[i]->GetMomentum().x()/MeV;
piy = vec[i]->GetMomentum().y()/MeV;
piz = vec[i]->GetMomentum().z()/MeV;
vec[i]->SetMomentum( cost*cosp*pix*MeV - sinp*piy+sint*cosp*piz*MeV,
cost*sinp*pix*MeV + cosp*piy+sint*sinp*piz*MeV,
-sint*pix*MeV + cost*piz*MeV );
}
}
else
{
if( pjz < 0.0 )
{
currentParticle.SetMomentum( -currentParticle.GetMomentum().z() );
targetParticle.SetMomentum( -targetParticle.GetMomentum().z() );
for( G4int i=0; i<vecLen; ++i )
vec[i]->SetMomentum( -vec[i]->GetMomentum().z() );
}
}
}
void G4ReactionDynamics::Rotate(
const G4double numberofFinalStateNucleons,
const G4ThreeVector &temp,
const G4ReactionProduct &modifiedOriginal, // Fermi motion & evap. effect included
const G4HadProjectile *originalIncident, // original incident particle
const G4Nucleus &targetNucleus,
G4ReactionProduct &currentParticle,
G4ReactionProduct &targetParticle,
G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
G4int &vecLen )
{
// derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
//
// Rotate in direction of z-axis, this does disturb in some way our
// inclusive distributions, but it is necessary for momentum conservation
//
const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
const G4double logWeight = std::log(atomicWeight);
G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
G4int i;
G4ThreeVector pseudoParticle[4];
for( i=0; i<4; ++i )pseudoParticle[i].set(0,0,0);
pseudoParticle[0] = currentParticle.GetMomentum()
+ targetParticle.GetMomentum();
for( i=0; i<vecLen; ++i )
pseudoParticle[0] = pseudoParticle[0] + (vec[i]->GetMomentum());
//
// Some smearing in transverse direction from Fermi motion
//
G4float pp, pp1;
G4double alekw, p;
G4double r1, r2, a1, ran1, ran2, xxh, exh, pxTemp, pyTemp, pzTemp;
r1 = twopi*G4UniformRand();
r2 = G4UniformRand();
a1 = std::sqrt(-2.0*std::log(r2));
ran1 = a1*std::sin(r1)*0.020*numberofFinalStateNucleons*GeV;
ran2 = a1*std::cos(r1)*0.020*numberofFinalStateNucleons*GeV;
G4ThreeVector frm(ran1, ran2, 0);
pseudoParticle[0] = pseudoParticle[0]+frm; // all particles + fermi
pseudoParticle[2] = temp; // original in cms system
pseudoParticle[3] = pseudoParticle[0];
pseudoParticle[1] = pseudoParticle[2].cross(pseudoParticle[3]);
G4double rotation = 2.*pi*G4UniformRand();
pseudoParticle[1] = pseudoParticle[1].rotate(rotation, pseudoParticle[3]);
pseudoParticle[2] = pseudoParticle[3].cross(pseudoParticle[1]);
for(G4int ii=1; ii<=3; ii++)
{
p = pseudoParticle[ii].mag();
if( p == 0.0 )
pseudoParticle[ii]= G4ThreeVector( 0.0, 0.0, 0.0 );
else
pseudoParticle[ii]= pseudoParticle[ii] * (1./p);
}
pxTemp = pseudoParticle[1].dot(currentParticle.GetMomentum());
pyTemp = pseudoParticle[2].dot(currentParticle.GetMomentum());
pzTemp = pseudoParticle[3].dot(currentParticle.GetMomentum());
currentParticle.SetMomentum( pxTemp, pyTemp, pzTemp );
pxTemp = pseudoParticle[1].dot(targetParticle.GetMomentum());
pyTemp = pseudoParticle[2].dot(targetParticle.GetMomentum());
pzTemp = pseudoParticle[3].dot(targetParticle.GetMomentum());
targetParticle.SetMomentum( pxTemp, pyTemp, pzTemp );
for( i=0; i<vecLen; ++i )
{
pxTemp = pseudoParticle[1].dot(vec[i]->GetMomentum());
pyTemp = pseudoParticle[2].dot(vec[i]->GetMomentum());
pzTemp = pseudoParticle[3].dot(vec[i]->GetMomentum());
vec[i]->SetMomentum( pxTemp, pyTemp, pzTemp );
}
//
// Rotate in direction of primary particle, subtract binding energies
// and make some further corrections if required
//
Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
G4double ekin;
G4double dekin = 0.0;
G4double ek1 = 0.0;
G4int npions = 0;
if( atomicWeight >= 1.5 ) // self-absorption in heavy molecules
{
// corrections for single particle spectra (shower particles)
//
const G4double alem[] = { 1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00 };
const G4double val0[] = { 0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65 };
alekw = std::log( originalIncident->GetKineticEnergy()/GeV );
exh = 1.0;
if( alekw > alem[0] ) // get energy bin
{
exh = val0[6];
for( G4int j=1; j<7; ++j )
{
if( alekw < alem[j] ) // use linear interpolation/extrapolation
{
G4double rcnve = (val0[j] - val0[j-1]) / (alem[j] - alem[j-1]);
exh = rcnve * alekw + val0[j-1] - rcnve * alem[j-1];
break;
}
}
exh = 1.0 - exh;
}
const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
ekin = currentParticle.GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
ekin = std::max( 1.0e-6, ekin );
xxh = 1.0;
if( (modifiedOriginal.GetDefinition() == aPiPlus ||
modifiedOriginal.GetDefinition() == aPiMinus) &&
currentParticle.GetDefinition() == aPiZero &&
G4UniformRand() <= logWeight) xxh = exh;
dekin += ekin*(1.0-xxh);
ekin *= xxh;
if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
++npions;
ek1 += ekin;
}
currentParticle.SetKineticEnergy( ekin*GeV );
pp = currentParticle.GetTotalMomentum()/MeV;
pp1 = currentParticle.GetMomentum().mag()/MeV;
if( pp1 < 0.001*MeV )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi = twopi*G4UniformRand();
currentParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
pp*sintheta*std::sin(phi)*MeV,
pp*costheta*MeV ) ;
}
else
currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
ekin = targetParticle.GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
ekin = std::max( 1.0e-6, ekin );
xxh = 1.0;
if( (modifiedOriginal.GetDefinition() == aPiPlus ||
modifiedOriginal.GetDefinition() == aPiMinus) &&
targetParticle.GetDefinition() == aPiZero &&
G4UniformRand() < logWeight) xxh = exh;
dekin += ekin*(1.0-xxh);
ekin *= xxh;
if (targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
++npions;
ek1 += ekin;
}
targetParticle.SetKineticEnergy( ekin*GeV );
pp = targetParticle.GetTotalMomentum()/MeV;
pp1 = targetParticle.GetMomentum().mag()/MeV;
if( pp1 < 0.001*MeV )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi = twopi*G4UniformRand();
targetParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
pp*sintheta*std::sin(phi)*MeV,
pp*costheta*MeV ) ;
}
else
targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
for( i=0; i<vecLen; ++i )
{
ekin = vec[i]->GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
ekin = std::max( 1.0e-6, ekin );
xxh = 1.0;
if( (modifiedOriginal.GetDefinition() == aPiPlus ||
modifiedOriginal.GetDefinition() == aPiMinus) &&
vec[i]->GetDefinition() == aPiZero &&
G4UniformRand() < logWeight) xxh = exh;
dekin += ekin*(1.0-xxh);
ekin *= xxh;
if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
++npions;
ek1 += ekin;
}
vec[i]->SetKineticEnergy( ekin*GeV );
pp = vec[i]->GetTotalMomentum()/MeV;
pp1 = vec[i]->GetMomentum().mag()/MeV;
if( pp1 < 0.001*MeV )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi = twopi*G4UniformRand();
vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*MeV,
pp*sintheta*std::sin(phi)*MeV,
pp*costheta*MeV ) ;
}
else
vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
}
}
if( (ek1 != 0.0) && (npions > 0) )
{
dekin = 1.0 + dekin/ek1;
//
// first do the incident particle
//
if (currentParticle.GetDefinition()->GetParticleSubType() == "pi")
{
currentParticle.SetKineticEnergy(
std::max( 0.001*MeV, dekin*currentParticle.GetKineticEnergy() ) );
pp = currentParticle.GetTotalMomentum()/MeV;
pp1 = currentParticle.GetMomentum().mag()/MeV;
if( pp1 < 0.001 )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi = twopi*G4UniformRand();
currentParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
pp*sintheta*std::sin(phi)*MeV,
pp*costheta*MeV ) ;
} else {
currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
}
}
if (targetParticle.GetDefinition()->GetParticleSubType() == "pi")
{
targetParticle.SetKineticEnergy(
std::max( 0.001*MeV, dekin*targetParticle.GetKineticEnergy() ) );
pp = targetParticle.GetTotalMomentum()/MeV;
pp1 = targetParticle.GetMomentum().mag()/MeV;
if( pp1 < 0.001 )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi = twopi*G4UniformRand();
targetParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
pp*sintheta*std::sin(phi)*MeV,
pp*costheta*MeV ) ;
} else {
targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
}
}
for( i=0; i<vecLen; ++i )
{
if (vec[i]->GetDefinition()->GetParticleSubType() == "pi")
{
vec[i]->SetKineticEnergy( std::max( 0.001*MeV, dekin*vec[i]->GetKineticEnergy() ) );
pp = vec[i]->GetTotalMomentum()/MeV;
pp1 = vec[i]->GetMomentum().mag()/MeV;
if( pp1 < 0.001 )
{
G4double costheta = 2.*G4UniformRand() - 1.;
G4double sintheta = std::sqrt(1. - costheta*costheta);
G4double phi = twopi*G4UniformRand();
vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*MeV,
pp*sintheta*std::sin(phi)*MeV,
pp*costheta*MeV ) ;
} else {
vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
}
}
} // for i
} // if (ek1 != 0)
}
void G4ReactionDynamics::AddBlackTrackParticles(
const G4double epnb, // GeV
const G4int npnb,
const G4double edta, // GeV
const G4int ndta,
const G4double sprob,
const G4double kineticMinimum, // GeV
const G4double kineticFactor, // GeV
const G4ReactionProduct &modifiedOriginal,
G4int PinNucleus,
G4int NinNucleus,
const G4Nucleus &targetNucleus,
G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
G4int &vecLen )
{
// derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
//
// npnb is number of proton/neutron black track particles
// ndta is the number of deuterons, tritons, and alphas produced
// epnb is the kinetic energy available for proton/neutron black track particles
// edta is the kinetic energy available for deuteron/triton/alpha particles
//
G4ParticleDefinition *aProton = G4Proton::Proton();
G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
G4ParticleDefinition *aTriton = G4Triton::Triton();
G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/MeV;
const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
const G4double ika1 = 3.6;
const G4double ika2 = 35.56;
const G4double ika3 = 6.45;
G4int i;
G4double pp;
G4double kinCreated = 0;
G4double cfa = 0.025*((atomicWeight-1.0)/120.0) * std::exp(-(atomicWeight-1.0)/120.0);
// First add protons and neutrons to final state
if (npnb > 0)
{
G4double backwardKinetic = 0.0;
G4int local_npnb = npnb;
for( i=0; i<npnb; ++i ) if( G4UniformRand() < sprob ) local_npnb--;
G4double local_epnb = epnb;
if (ndta == 0) local_epnb += edta; // Retrieve unused kinetic energy
G4double ekin = local_epnb/std::max(1,local_npnb);
for( i=0; i<local_npnb; ++i )
{
G4ReactionProduct * p1 = new G4ReactionProduct();
if( backwardKinetic > local_epnb )
{
delete p1;
break;
}
G4double ran = G4UniformRand();
G4double kinetic = -ekin*std::log(ran) - cfa*(1.0+0.5*normal());
if( kinetic < 0.0 )kinetic = -0.010*std::log(ran);
backwardKinetic += kinetic;
if( backwardKinetic > local_epnb )
kinetic = std::max( kineticMinimum, local_epnb-(backwardKinetic-kinetic) );
if (G4UniformRand() > (1.0-atomicNumber/atomicWeight)) {
// Boil off a proton if there are any left, otherwise a neutron
if (PinNucleus > 0) {
p1->SetDefinition( aProton );
PinNucleus--;
} else if (NinNucleus > 0) {
p1->SetDefinition( aNeutron );
NinNucleus--;
} else {
delete p1;
break; // no nucleons left in nucleus
}
} else {
// Boil off a neutron if there are any left, otherwise a proton
if (NinNucleus > 0) {
p1->SetDefinition( aNeutron );
NinNucleus--;
} else if (PinNucleus > 0) {
p1->SetDefinition( aProton );
PinNucleus--;
} else {
delete p1;
break; // no nucleons left in nucleus
}
}
vec.SetElement( vecLen, p1 );
G4double cost = G4UniformRand() * 2.0 - 1.0;
G4double sint = std::sqrt(std::fabs(1.0-cost*cost));
G4double phi = twopi * G4UniformRand();
vec[vecLen]->SetNewlyAdded( true );
vec[vecLen]->SetKineticEnergy( kinetic*GeV );
kinCreated+=kinetic;
pp = vec[vecLen]->GetTotalMomentum()/MeV;
vec[vecLen]->SetMomentum( pp*sint*std::sin(phi)*MeV,
pp*sint*std::cos(phi)*MeV,
pp*cost*MeV );
vecLen++;
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
if (NinNucleus > 0) {
if( (atomicWeight >= 10.0) && (ekOriginal <= 2.0*GeV) )
{
G4double ekw = ekOriginal/GeV;
G4int ika, kk = 0;
if( ekw > 1.0 )ekw *= ekw;
ekw = std::max( 0.1, ekw );
ika = G4int(ika1*std::exp((atomicNumber*atomicNumber/
atomicWeight-ika2)/ika3)/ekw);
if( ika > 0 )
{
for( i=(vecLen-1); i>=0; --i )
{
if( (vec[i]->GetDefinition() == aProton) && vec[i]->GetNewlyAdded() )
{
vec[i]->SetDefinitionAndUpdateE( aNeutron ); // modified 22-Oct-97
PinNucleus++;
NinNucleus--;
if( ++kk > ika )break;
}
}
}
}
} // if (NinNucleus >0)
} // if (npnb > 0)
// Next try to add deuterons, tritons and alphas to final state
if (ndta > 0)
{
G4double backwardKinetic = 0.0;
G4int local_ndta=ndta;
for( i=0; i<ndta; ++i )if( G4UniformRand() < sprob )local_ndta--;
G4double local_edta = edta;
if (npnb == 0) local_edta += epnb; // Retrieve unused kinetic energy
G4double ekin = local_edta/std::max(1,local_ndta);
for( i=0; i<local_ndta; ++i )
{
G4ReactionProduct *p2 = new G4ReactionProduct();
if( backwardKinetic > local_edta )
{
delete p2;
break;
}
G4double ran = G4UniformRand();
G4double kinetic = -ekin*std::log(ran)-cfa*(1.+0.5*normal());
if( kinetic < 0.0 )kinetic = kineticFactor*std::log(ran);
backwardKinetic += kinetic;
if( backwardKinetic > local_edta )kinetic = local_edta-(backwardKinetic-kinetic);
if( kinetic < 0.0 )kinetic = kineticMinimum;
G4double cost = 2.0*G4UniformRand() - 1.0;
G4double sint = std::sqrt(std::max(0.0,(1.0-cost*cost)));
G4double phi = twopi*G4UniformRand();
ran = G4UniformRand();
if (ran < 0.60) {
if (PinNucleus > 0 && NinNucleus > 0) {
p2->SetDefinition( aDeuteron );
PinNucleus--;
NinNucleus--;
} else if (NinNucleus > 0) {
p2->SetDefinition( aNeutron );
NinNucleus--;
} else if (PinNucleus > 0) {
p2->SetDefinition( aProton );
PinNucleus--;
} else {
delete p2;
break;
}
} else if (ran < 0.90) {
if (PinNucleus > 0 && NinNucleus > 1) {
p2->SetDefinition( aTriton );
PinNucleus--;
NinNucleus -= 2;
} else if (PinNucleus > 0 && NinNucleus > 0) {
p2->SetDefinition( aDeuteron );
PinNucleus--;
NinNucleus--;
} else if (NinNucleus > 0) {
p2->SetDefinition( aNeutron );
NinNucleus--;
} else if (PinNucleus > 0) {
p2->SetDefinition( aProton );
PinNucleus--;
} else {
delete p2;
break;
}
} else {
if (PinNucleus > 1 && NinNucleus > 1) {
p2->SetDefinition( anAlpha );
PinNucleus -= 2;
NinNucleus -= 2;
} else if (PinNucleus > 0 && NinNucleus > 1) {
p2->SetDefinition( aTriton );
PinNucleus--;
NinNucleus -= 2;
} else if (PinNucleus > 0 && NinNucleus > 0) {
p2->SetDefinition( aDeuteron );
PinNucleus--;
NinNucleus--;
} else if (NinNucleus > 0) {
p2->SetDefinition( aNeutron );
NinNucleus--;
} else if (PinNucleus > 0) {
p2->SetDefinition( aProton );
PinNucleus--;
} else {
delete p2;
break;
}
}
vec.SetElement( vecLen, p2 );
vec[vecLen]->SetNewlyAdded( true );
vec[vecLen]->SetKineticEnergy( kinetic*GeV );
kinCreated+=kinetic;
pp = vec[vecLen]->GetTotalMomentum()/MeV;
vec[vecLen++]->SetMomentum( pp*sint*std::sin(phi)*MeV,
pp*sint*std::cos(phi)*MeV,
pp*cost*MeV );
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
} // if (ndta > 0)
// G4double delta = epnb+edta - kinCreated;
}
std::pair<G4int, G4int> G4ReactionDynamics::GetFinalStateNucleons(
const G4DynamicParticle* originalTarget,
const G4FastVector<G4ReactionProduct,GHADLISTSIZE>& vec,
const G4int& vecLen)
{
// Get number of protons and neutrons removed from the target nucleus
G4int protonsRemoved = 0;
G4int neutronsRemoved = 0;
if (originalTarget->GetDefinition()->GetParticleName() == "proton")
protonsRemoved++;
else
neutronsRemoved++;
G4String secName;
for (G4int i = 0; i < vecLen; i++) {
secName = vec[i]->GetDefinition()->GetParticleName();
if (secName == "proton") {
protonsRemoved++;
} else if (secName == "neutron") {
neutronsRemoved++;
} else if (secName == "anti_proton") {
protonsRemoved--;
} else if (secName == "anti_neutron") {
neutronsRemoved--;
}
}
return std::pair<G4int, G4int>(protonsRemoved, neutronsRemoved);
}
void G4ReactionDynamics::MomentumCheck(
const G4ReactionProduct &modifiedOriginal,
G4ReactionProduct &currentParticle,
G4ReactionProduct &targetParticle,
G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
G4int &vecLen )
{
const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/MeV;
G4double testMomentum = currentParticle.GetMomentum().mag()/MeV;
G4double pMass;
if( testMomentum >= pOriginal )
{
pMass = currentParticle.GetMass()/MeV;
currentParticle.SetTotalEnergy(
std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
currentParticle.SetMomentum(
currentParticle.GetMomentum() * (pOriginal/testMomentum) );
}
testMomentum = targetParticle.GetMomentum().mag()/MeV;
if( testMomentum >= pOriginal )
{
pMass = targetParticle.GetMass()/MeV;
targetParticle.SetTotalEnergy(
std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
targetParticle.SetMomentum(
targetParticle.GetMomentum() * (pOriginal/testMomentum) );
}
for( G4int i=0; i<vecLen; ++i )
{
testMomentum = vec[i]->GetMomentum().mag()/MeV;
if( testMomentum >= pOriginal )
{
pMass = vec[i]->GetMass()/MeV;
vec[i]->SetTotalEnergy(
std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
vec[i]->SetMomentum( vec[i]->GetMomentum() * (pOriginal/testMomentum) );
}
}
}
void G4ReactionDynamics::ProduceStrangeParticlePairs(
G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
G4int &vecLen,
const G4ReactionProduct &modifiedOriginal,
const G4DynamicParticle *originalTarget,
G4ReactionProduct &currentParticle,
G4ReactionProduct &targetParticle,
G4bool &incidentHasChanged,
G4bool &targetHasChanged )
{
// derived from original FORTRAN code STPAIR by H. Fesefeldt (16-Dec-1987)
//
// Choose charge combinations K+ K-, K+ K0B, K0 K0B, K0 K-,
// K+ Y0, K0 Y+, K0 Y-
// For antibaryon induced reactions half of the cross sections KB YB
// pairs are produced. Charge is not conserved, no experimental data available
// for exclusive reactions, therefore some average behaviour assumed.
// The ratio L/SIGMA is taken as 3:1 (from experimental low energy)
//
if( vecLen == 0 )return;
//
// the following protects against annihilation processes
//
if( currentParticle.GetMass() == 0.0 || targetParticle.GetMass() == 0.0 )return;
const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass()/GeV;
G4double targetMass = originalTarget->GetDefinition()->GetPDGMass()/GeV;
G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
targetMass*targetMass +
2.0*targetMass*etOriginal ); // GeV
G4double currentMass = currentParticle.GetMass()/GeV;
G4double availableEnergy = centerofmassEnergy-(targetMass+currentMass);
if( availableEnergy <= 1.0 )return;
G4ParticleDefinition *aProton = G4Proton::Proton();
G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
G4ParticleDefinition *aLambda = G4Lambda::Lambda();
G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
const G4double protonMass = aProton->GetPDGMass()/GeV;
const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass()/GeV;
//
// determine the center of mass energy bin
//
const G4double avrs[] = {3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.};
G4int ibin, i3, i4;
G4double avk, avy, avn, ran;
G4int i = 1;
while( (i<12) && (centerofmassEnergy>avrs[i]) )++i;
if( i == 12 )
ibin = 11;
else
ibin = i;
//
// the fortran code chooses a random replacement of produced kaons
// but does not take into account charge conservation
//
if( vecLen == 1 ) // we know that vecLen > 0
{
i3 = 0;
i4 = 1; // note that we will be adding a new secondary particle in this case only
}
else // otherwise 0 <= i3,i4 < vecLen
{
ran = G4UniformRand();
while( ran == 1.0 )ran = G4UniformRand();
i4 = i3 = G4int( vecLen*ran );
while( i3 == i4 )
{
ran = G4UniformRand();
while( ran == 1.0 )ran = G4UniformRand();
i4 = G4int( vecLen*ran );
}
}
//
// use linear interpolation or extrapolation by y=centerofmassEnergy*x+b
//
const G4double avkkb[] = { 0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075,
0.0975, 0.123, 0.28, 0.398, 0.495, 0.573 };
const G4double avky[] = { 0.005, 0.03, 0.064, 0.095, 0.115, 0.13,
0.145, 0.155, 0.20, 0.205, 0.210, 0.212 };
const G4double avnnb[] = { 0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02,
0.04, 0.05, 0.12, 0.15, 0.18, 0.20 };
avk = (std::log(avkkb[ibin])-std::log(avkkb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
/(avrs[ibin]-avrs[ibin-1]) + std::log(avkkb[ibin-1]);
avk = std::exp(avk);
avy = (std::log(avky[ibin])-std::log(avky[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
/(avrs[ibin]-avrs[ibin-1]) + std::log(avky[ibin-1]);
avy = std::exp(avy);
avn = (std::log(avnnb[ibin])-std::log(avnnb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
/(avrs[ibin]-avrs[ibin-1]) + std::log(avnnb[ibin-1]);
avn = std::exp(avn);
if( avk+avy+avn <= 0.0 )return;
if( currentMass < protonMass )avy /= 2.0;
if( targetMass < protonMass )avy = 0.0;
avy += avk+avn;
avk += avn;
ran = G4UniformRand();
if( ran < avn )
{
if( availableEnergy < 2.0 )return;
if( vecLen == 1 ) // add a new secondary
{
G4ReactionProduct *p1 = new G4ReactionProduct;
if( G4UniformRand() < 0.5 )
{
vec[0]->SetDefinition( aNeutron );
p1->SetDefinition( anAntiNeutron );
(G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
vec[0]->SetMayBeKilled(false);
p1->SetMayBeKilled(false);
}
else
{
vec[0]->SetDefinition( aProton );
p1->SetDefinition( anAntiProton );
(G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
vec[0]->SetMayBeKilled(false);
p1->SetMayBeKilled(false);
}
vec.SetElement( vecLen++, p1 );
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
else
{ // replace two secondaries
if( G4UniformRand() < 0.5 )
{
vec[i3]->SetDefinition( aNeutron );
vec[i4]->SetDefinition( anAntiNeutron );
vec[i3]->SetMayBeKilled(false);
vec[i4]->SetMayBeKilled(false);
}
else
{
vec[i3]->SetDefinition( aProton );
vec[i4]->SetDefinition( anAntiProton );
vec[i3]->SetMayBeKilled(false);
vec[i4]->SetMayBeKilled(false);
}
}
}
else if( ran < avk )
{
if( availableEnergy < 1.0 )return;
const G4double kkb[] = { 0.2500, 0.3750, 0.5000, 0.5625, 0.6250,
0.6875, 0.7500, 0.8750, 1.000 };
const G4int ipakkb1[] = { 10, 10, 10, 11, 11, 12, 12, 11, 12 };
const G4int ipakkb2[] = { 13, 11, 12, 11, 12, 11, 12, 13, 13 };
ran = G4UniformRand();
i = 0;
while( (i<9) && (ran>=kkb[i]) )++i;
if( i == 9 )return;
//
// ipakkb[] = { 10,13, 10,11, 10,12, 11,11, 11,12, 12,11, 12,12, 11,13, 12,13 };
// charge + - + 0 + 0 0 0 0 0 0 0 0 0 0 - 0 -
//
switch( ipakkb1[i] )
{
case 10:
vec[i3]->SetDefinition( aKaonPlus );
vec[i3]->SetMayBeKilled(false);
break;
case 11:
vec[i3]->SetDefinition( aKaonZS );
vec[i3]->SetMayBeKilled(false);
break;
case 12:
vec[i3]->SetDefinition( aKaonZL );
vec[i3]->SetMayBeKilled(false);
break;
}
if( vecLen == 1 ) // add a secondary
{
G4ReactionProduct *p1 = new G4ReactionProduct;
switch( ipakkb2[i] )
{
case 11:
p1->SetDefinition( aKaonZS );
p1->SetMayBeKilled(false);
break;
case 12:
p1->SetDefinition( aKaonZL );
p1->SetMayBeKilled(false);
break;
case 13:
p1->SetDefinition( aKaonMinus );
p1->SetMayBeKilled(false);
break;
}
(G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
vec.SetElement( vecLen++, p1 );
// DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
}
else // replace
{
switch( ipakkb2[i] )
{
case 11:
vec[i4]->SetDefinition( aKaonZS );
vec[i4]->SetMayBeKilled(false);
break;
case 12:
vec[i4]->SetDefinition( aKaonZL );
vec[i4]->SetMayBeKilled(false);
break;
case 13:
vec[i4]->SetDefinition( aKaonMinus );
vec[i4]->SetMayBeKilled(false);
break;
}
}
}
else if( ran < avy )
{
if( availableEnergy < 1.6 )return;
const G4double ky[] = { 0.200, 0.300, 0.400, 0.550, 0.625, 0.700,
0.800, 0.850, 0.900, 0.950, 0.975, 1.000 };
const G4int ipaky1[] = { 18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22 };
const G4int ipaky2[] = { 10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12 };
const G4int ipakyb1[] = { 19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25 };
const G4int ipakyb2[] = { 13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11 };
ran = G4UniformRand();
i = 0;
while( (i<12) && (ran>ky[i]) )++i;
if( i == 12 )return;
if( (currentMass<protonMass) || (G4UniformRand()<0.5) )
{
// ipaky[] = { 18,10, 18,11, 18,12, 20,10, 20,11, 20,12,
// 0 + 0 0 0 0 + + + 0 + 0
//
// 21,10, 21,11, 21,12, 22,10, 22,11, 22,12 }
// 0 + 0 0 0 0 - + - 0 - 0
switch( ipaky1[i] )
{
case 18:
targetParticle.SetDefinition( aLambda );
break;
case 20:
targetParticle.SetDefinition( aSigmaPlus );
break;
case 21:
targetParticle.SetDefinition( aSigmaZero );
break;
case 22:
targetParticle.SetDefinition( aSigmaMinus );
break;
}
targetHasChanged = true;
switch( ipaky2[i] )
{
case 10:
vec[i3]->SetDefinition( aKaonPlus );
vec[i3]->SetMayBeKilled(false);
break;
case 11:
vec[i3]->SetDefinition( aKaonZS );
vec[i3]->SetMayBeKilled(false);
break;
case 12:
vec[i3]->SetDefinition( aKaonZL );
vec[i3]->SetMayBeKilled(false);
break;
}
}
else // (currentMass >= protonMass) && (G4UniformRand() >= 0.5)
{
// ipakyb[] = { 19,13, 19,12, 19,11, 23,13, 23,12, 23,11,
// 24,13, 24,12, 24,11, 25,13, 25,12, 25,11 };
if( (currentParticle.GetDefinition() == anAntiProton) ||
(currentParticle.GetDefinition() == anAntiNeutron) ||
(currentParticle.GetDefinition() == anAntiLambda) ||
(currentMass > sigmaMinusMass) )
{
switch( ipakyb1[i] )
{
case 19:
currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
break;
case 23:
currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
break;
case 24:
currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
break;
case 25:
currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
break;
}
incidentHasChanged = true;
switch( ipakyb2[i] )
{
case 11:
vec[i3]->SetDefinition( aKaonZS );
vec[i3]->SetMayBeKilled(false);
break;
case 12:
vec[i3]->SetDefinition( aKaonZL );
vec[i3]->SetMayBeKilled(false);
break;
case 13:
vec[i3]->SetDefinition( aKaonMinus );
vec[i3]->SetMayBeKilled(false);
break;
}
}
else
{
switch( ipaky1[i] )
{
case 18:
currentParticle.SetDefinitionAndUpdateE( aLambda );
break;
case 20:
currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
break;
case 21:
currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
break;
case 22:
currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
break;
}
incidentHasChanged = true;
switch( ipaky2[i] )
{
case 10:
vec[i3]->SetDefinition( aKaonPlus );
vec[i3]->SetMayBeKilled(false);
break;
case 11:
vec[i3]->SetDefinition( aKaonZS );
vec[i3]->SetMayBeKilled(false);
break;
case 12:
vec[i3]->SetDefinition( aKaonZL );
vec[i3]->SetMayBeKilled(false);
break;
}
}
}
}
else return;
//
// check the available energy
// if there is not enough energy for kkb/ky pair production
// then reduce the number of secondary particles
// NOTE:
// the number of secondaries may have been changed
// the incident and/or target particles may have changed
// charge conservation is ignored (as well as strangness conservation)
//
currentMass = currentParticle.GetMass()/GeV;
targetMass = targetParticle.GetMass()/GeV;
G4double energyCheck = centerofmassEnergy-(currentMass+targetMass);
for( i=0; i<vecLen; ++i )
{
energyCheck -= vec[i]->GetMass()/GeV;
if( energyCheck < 0.0 ) // chop off the secondary List
{
vecLen = std::max( 0, --i ); // looks like a memory leak @@@@@@@@@@@@
G4int j;
for(j=i; j<vecLen; j++) delete vec[j];
break;
}
}
return;
}
void
G4ReactionDynamics::NuclearReaction(
G4FastVector<G4ReactionProduct,4> &vec,
G4int &vecLen,
const G4HadProjectile *originalIncident,
const G4Nucleus &targetNucleus,
const G4double theAtomicMass,
const G4double *mass )
{
// derived from original FORTRAN code NUCREC by H. Fesefeldt (12-Feb-1987)
//
// Nuclear reaction kinematics at low energies
//
G4ParticleDefinition *aGamma = G4Gamma::Gamma();
G4ParticleDefinition *aProton = G4Proton::Proton();
G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
G4ParticleDefinition *aTriton = G4Triton::Triton();
G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
const G4double aProtonMass = aProton->GetPDGMass()/MeV;
const G4double aNeutronMass = aNeutron->GetPDGMass()/MeV;
const G4double aDeuteronMass = aDeuteron->GetPDGMass()/MeV;
const G4double aTritonMass = aTriton->GetPDGMass()/MeV;
const G4double anAlphaMass = anAlpha->GetPDGMass()/MeV;
G4ReactionProduct currentParticle;
currentParticle = *originalIncident;
//
// Set beam particle, take kinetic energy of current particle as the
// fundamental quantity. Due to the difficult kinematic, all masses have to
// be assigned the best measured values
//
G4double p = currentParticle.GetTotalMomentum();
G4double pp = currentParticle.GetMomentum().mag();
if( pp <= 0.001*MeV )
{
G4double phinve = twopi*G4UniformRand();
G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
p*std::sin(rthnve)*std::sin(phinve),
p*std::cos(rthnve) );
}
else
currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
//
// calculate Q-value of reactions
//
G4double currentKinetic = currentParticle.GetKineticEnergy()/MeV;
G4double currentMass = currentParticle.GetDefinition()->GetPDGMass()/MeV;
G4double qv = currentKinetic + theAtomicMass + currentMass;
G4double qval[9];
qval[0] = qv - mass[0];
qval[1] = qv - mass[1] - aNeutronMass;
qval[2] = qv - mass[2] - aProtonMass;
qval[3] = qv - mass[3] - aDeuteronMass;
qval[4] = qv - mass[4] - aTritonMass;
qval[5] = qv - mass[5] - anAlphaMass;
qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
qval[8] = qv - mass[8] - aProtonMass - aProtonMass;
if( currentParticle.GetDefinition() == aNeutron )
{
const G4double A = G4double(targetNucleus.GetA_asInt()); // atomic weight
if( G4UniformRand() > ((A-1.0)/230.0)*((A-1.0)/230.0) )
qval[0] = 0.0;
if( G4UniformRand() >= currentKinetic/7.9254*A )
qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
}
else
qval[0] = 0.0;
G4int i;
qv = 0.0;
for( i=0; i<9; ++i )
{
if( mass[i] < 500.0*MeV )qval[i] = 0.0;
if( qval[i] < 0.0 )qval[i] = 0.0;
qv += qval[i];
}
G4double qv1 = 0.0;
G4double ran = G4UniformRand();
G4int index;
for( index=0; index<9; ++index )
{
if( qval[index] > 0.0 )
{
qv1 += qval[index]/qv;
if( ran <= qv1 )break;
}
}
if( index == 9 ) // loop continued to the end
{
throw G4HadronicException(__FILE__, __LINE__,
"G4ReactionDynamics::NuclearReaction: inelastic reaction kinematically not possible");
}
G4double ke = currentParticle.GetKineticEnergy()/GeV;
G4int nt = 2;
if( (index>=6) || (G4UniformRand()<std::min(0.5,ke*10.0)) )nt = 3;
G4ReactionProduct **v = new G4ReactionProduct * [3];
v[0] = new G4ReactionProduct;
v[1] = new G4ReactionProduct;
v[2] = new G4ReactionProduct;
v[0]->SetMass( mass[index]*MeV );
switch( index )
{
case 0:
v[1]->SetDefinition( aGamma );
v[2]->SetDefinition( aGamma );
break;
case 1:
v[1]->SetDefinition( aNeutron );
v[2]->SetDefinition( aGamma );
break;
case 2:
v[1]->SetDefinition( aProton );
v[2]->SetDefinition( aGamma );
break;
case 3:
v[1]->SetDefinition( aDeuteron );
v[2]->SetDefinition( aGamma );
break;
case 4:
v[1]->SetDefinition( aTriton );
v[2]->SetDefinition( aGamma );
break;
case 5:
v[1]->SetDefinition( anAlpha );
v[2]->SetDefinition( aGamma );
break;
case 6:
v[1]->SetDefinition( aNeutron );
v[2]->SetDefinition( aNeutron );
break;
case 7:
v[1]->SetDefinition( aNeutron );
v[2]->SetDefinition( aProton );
break;
case 8:
v[1]->SetDefinition( aProton );
v[2]->SetDefinition( aProton );
break;
}
//
// calculate centre of mass energy
//
G4ReactionProduct pseudo1;
pseudo1.SetMass( theAtomicMass*MeV );
pseudo1.SetTotalEnergy( theAtomicMass*MeV );
G4ReactionProduct pseudo2 = currentParticle + pseudo1;
pseudo2.SetMomentum( pseudo2.GetMomentum() * (-1.0) );
//
// use phase space routine in centre of mass system
//
G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
tempV.Initialize( nt );
G4int tempLen = 0;
tempV.SetElement( tempLen++, v[0] );
tempV.SetElement( tempLen++, v[1] );
if( nt == 3 )tempV.SetElement( tempLen++, v[2] );
G4bool constantCrossSection = true;
GenerateNBodyEvent( pseudo2.GetMass()/MeV, constantCrossSection, tempV, tempLen );
v[0]->Lorentz( *v[0], pseudo2 );
v[1]->Lorentz( *v[1], pseudo2 );
if( nt == 3 )v[2]->Lorentz( *v[2], pseudo2 );
G4bool particleIsDefined = false;
if( v[0]->GetMass()/MeV - aProtonMass < 0.1 )
{
v[0]->SetDefinition( aProton );
particleIsDefined = true;
}
else if( v[0]->GetMass()/MeV - aNeutronMass < 0.1 )
{
v[0]->SetDefinition( aNeutron );
particleIsDefined = true;
}
else if( v[0]->GetMass()/MeV - aDeuteronMass < 0.1 )
{
v[0]->SetDefinition( aDeuteron );
particleIsDefined = true;
}
else if( v[0]->GetMass()/MeV - aTritonMass < 0.1 )
{
v[0]->SetDefinition( aTriton );
particleIsDefined = true;
}
else if( v[0]->GetMass()/MeV - anAlphaMass < 0.1 )
{
v[0]->SetDefinition( anAlpha );
particleIsDefined = true;
}
currentParticle.SetKineticEnergy(
std::max( 0.001, currentParticle.GetKineticEnergy()/MeV ) );
p = currentParticle.GetTotalMomentum();
pp = currentParticle.GetMomentum().mag();
if( pp <= 0.001*MeV )
{
G4double phinve = twopi*G4UniformRand();
G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
p*std::sin(rthnve)*std::sin(phinve),
p*std::cos(rthnve) );
}
else
currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
if( particleIsDefined )
{
v[0]->SetKineticEnergy(
std::max( 0.001, 0.5*G4UniformRand()*v[0]->GetKineticEnergy()/MeV ) );
p = v[0]->GetTotalMomentum();
pp = v[0]->GetMomentum().mag();
if( pp <= 0.001*MeV )
{
G4double phinve = twopi*G4UniformRand();
G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
v[0]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
p*std::sin(rthnve)*std::sin(phinve),
p*std::cos(rthnve) );
}
else
v[0]->SetMomentum( v[0]->GetMomentum() * (p/pp) );
}
if( (v[1]->GetDefinition() == aDeuteron) ||
(v[1]->GetDefinition() == aTriton) ||
(v[1]->GetDefinition() == anAlpha) )
v[1]->SetKineticEnergy(
std::max( 0.001, 0.5*G4UniformRand()*v[1]->GetKineticEnergy()/MeV ) );
else
v[1]->SetKineticEnergy( std::max( 0.001, v[1]->GetKineticEnergy()/MeV ) );
p = v[1]->GetTotalMomentum();
pp = v[1]->GetMomentum().mag();
if( pp <= 0.001*MeV )
{
G4double phinve = twopi*G4UniformRand();
G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
v[1]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
p*std::sin(rthnve)*std::sin(phinve),
p*std::cos(rthnve) );
}
else
v[1]->SetMomentum( v[1]->GetMomentum() * (p/pp) );
if( nt == 3 )
{
if( (v[2]->GetDefinition() == aDeuteron) ||
(v[2]->GetDefinition() == aTriton) ||
(v[2]->GetDefinition() == anAlpha) )
v[2]->SetKineticEnergy(
std::max( 0.001, 0.5*G4UniformRand()*v[2]->GetKineticEnergy()/MeV ) );
else
v[2]->SetKineticEnergy( std::max( 0.001, v[2]->GetKineticEnergy()/MeV ) );
p = v[2]->GetTotalMomentum();
pp = v[2]->GetMomentum().mag();
if( pp <= 0.001*MeV )
{
G4double phinve = twopi*G4UniformRand();
G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
v[2]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
p*std::sin(rthnve)*std::sin(phinve),
p*std::cos(rthnve) );
}
else
v[2]->SetMomentum( v[2]->GetMomentum() * (p/pp) );
}
G4int del;
for(del=0; del<vecLen; del++) delete vec[del];
vecLen = 0;
if( particleIsDefined )
{
vec.SetElement( vecLen++, v[0] );
}
else
{
delete v[0];
}
vec.SetElement( vecLen++, v[1] );
if( nt == 3 )
{
vec.SetElement( vecLen++, v[2] );
}
else
{
delete v[2];
}
delete [] v;
return;
}
/* end of file */
Reference in New Issue View Git Blame Copy Permalink