43 #define MAX_SECONDARIES 100
47 G4AntiNeutronAnnihilationAtRest::G4AntiNeutronAnnihilationAtRest(
const G4String& processName,
105 return ( &particle == pdefAntiNeutron );
138 G4cout <<
"G4AntiNeutronAnnihilationAtRestProcess::AtRestGetPhysicalInteractionLength ";
175 for (
G4int i1=0; i1 < numberOfElements; i1++ )
177 normalization += theAtomicNumberDensity[i1] ;
182 for (
G4int i2=0; i2 < numberOfElements; i2++ )
184 runningSum += theAtomicNumberDensity[i2];
186 if (random<=runningSum)
188 targetCharge =
G4double( ((*theElementVector)[i2])->GetZ());
189 targetAtomicMass = (*theElementVector)[i2]->GetN();
192 if (random>runningSum)
194 targetCharge =
G4double( ((*theElementVector)[numberOfElements-1])->GetZ());
195 targetAtomicMass = (*theElementVector)[numberOfElements-1]->GetN();
199 G4cout <<
"G4AntiNeutronAnnihilationAtRest::AtRestDoIt is invoked " <<
G4endl;
207 GenerateSecondaries();
211 for (
G4int isec = 0; isec < ngkine; isec++ ) {
216 localtime = globalTime + gkin[isec].
GetTOF();
218 G4Track* aNewTrack =
new G4Track( aNewParticle, localtime*
s, position );
237 void G4AntiNeutronAnnihilationAtRest::GenerateSecondaries()
253 result.
SetMass( massAntiNeutron );
260 AntiNeutronAnnihilation(&nopt);
279 for (l = 1; l <= ntot; ++l) {
285 gkin[ngkine] = eve[
index];
286 gkin[ngkine].
SetTOF( eve[index].GetTOF() * 5e-11 );
305 void G4AntiNeutronAnnihilationAtRest::Poisso(
G4float xav,
G4int *iran)
320 ran1 = xav + ran1 * std::sqrt(xav);
334 for (i = 1; i <= fivex; ++i) {
337 rrr = std::pow(xav,
G4float(i)) / NFac(i);
341 rrr = std::exp(i * std::log(xav) -
354 p1 = xav * std::exp(-
G4double(xav));
380 G4int G4AntiNeutronAnnihilationAtRest::NFac(
G4int n)
395 for (i = 2; i <= j; ++i) {
404 void G4AntiNeutronAnnihilationAtRest::Normal(
G4float *ran)
412 for (i = 1; i <= 12; ++i) {
419 void G4AntiNeutronAnnihilationAtRest::AntiNeutronAnnihilation(
G4int *nopt)
451 pv[1].
SetMass( massAntiNeutron );
470 rmnve1 = massPionPlus;
471 rmnve2 = massPionMinus;
473 rmnve1 = massPionZero;
476 rmnve2 = massPionZero;
479 rmnve2 = massPionZero;
487 ek = massNeutron + massAntiNeutron - rmnve1 - rmnve2;
494 en = ek + (rmnve1 + rmnve2) /
G4float(2.);
495 r__1 = en * en - rmnve1 * rmnve2;
496 pcm = r__1 > 0 ? std::sqrt(r__1) : 0;
507 pv[3].
SetEnergy( std::sqrt(pv[3].GetMass()*pv[3].GetMass()+pcm*pcm) );
531 if (targetAtomicMass >=
G4float(1.5)) {
538 black = (targ *
G4float(1.25) +
539 G4float(1.5)) * evapEnergy1 / (evapEnergy1 + evapEnergy3);
542 nbl =
G4int(targetAtomicMass - targ);
550 for (i = 1; i <= nbl; ++i) {
559 ekin1 = -
G4double(ekin) * std::log(ran1) -
562 ekin1 = std::log(ran1) *
G4float(-.01);
567 ekin1 = tex - (ekin2 - ekin1);
573 pnrat =
G4float(1.) - targetCharge / targetAtomicMass;
595 for (i = 1; i <=
nt; ++i) {
597 if (pv[ii].GetParticleDef() != pdefProton) {
617 (evapEnergy1 + evapEnergy3);
625 for (i = 1; i <= nbl; ++i) {
634 ekin1 = -
G4double(ekin) * std::log(ran1) -
637 ekin1 = std::log(ran1) *
G4float(-.01);
642 ekin1 = tex - (ekin2 - ekin1);
669 for (i = 3; i <=
nt; ++i) {
693 if (targetAtomicMass >=
G4float(1.5)) {
707 cfa =
G4float(.13043478260869565);
708 cfa = cfa * std::log(ekin1) +
G4float(.35);
713 atno3 = targetAtomicMass;
725 gfa = (targetAtomicMass -
G4float(1.)) /
728 evapEnergy1 = ret_val * fpdiv;
729 evapEnergy3 = ret_val - evapEnergy1;
736 evapEnergy1 *= ran1 * gfa +
G4float(1.);
737 if (evapEnergy1 <
G4float(0.)) {
740 evapEnergy3 *= ran2 * gfa +
G4float(1.);
741 if (evapEnergy3 <
G4float(0.)) {
744 while ((ret_val = evapEnergy1 + evapEnergy3) >= ek1) {
G4double condition(const G4ErrorSymMatrix &m)
void DeRegisterExtraProcess(G4VProcess *)
G4double GetMeanLifeTime(const G4Track &, G4ForceCondition *)
void SetParticleDef(G4ParticleDefinition *c)
void SetMomentum(const G4ThreeVector &momentum)
std::vector< G4Element * > G4ElementVector
static G4HadronicProcessStore * Instance()
const G4DynamicParticle * GetDynamicParticle() const
~G4AntiNeutronAnnihilationAtRest()
G4ParticleDefinition * GetParticleDef()
const G4String & GetName() const
const G4ThreeVector & GetPosition() const
void DumpInfo(G4int mode=0) const
G4GHEKinematicsVector * GetSecondaryKinematics()
G4double theNumberOfInteractionLengthLeft
void SetTouchableHandle(const G4TouchableHandle &apValue)
void BuildPhysicsTable(const G4ParticleDefinition &)
virtual void ResetNumberOfInteractionLengthLeft()
#define G4HadronicDeprecate(name)
const G4ElementVector * GetElementVector() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetKineticEnergyAndUpdate(G4double ekin)
G4VParticleChange * AtRestDoIt(const G4Track &, const G4Step &)
void RegisterParticleForExtraProcess(G4VProcess *, const G4ParticleDefinition *)
G4GLOB_DLL std::ostream G4cout
G4int GetNumberOfSecondaries()
G4double AtRestGetPhysicalInteractionLength(const G4Track &, G4ForceCondition *)
G4double currentInteractionLength
void SetProcessSubType(G4int)
G4double GetGlobalTime() const
const G4String & GetProcessName() const
const G4double * GetAtomicNumDensityVector() const
const G4TouchableHandle & GetTouchableHandle() const
G4Material * GetMaterial() const
void SetEnergyAndUpdate(G4double e)
void RegisterExtraProcess(G4VProcess *)
G4bool IsApplicable(const G4ParticleDefinition &)
void PreparePhysicsTable(const G4ParticleDefinition &)
virtual void Initialize(const G4Track &)
void SetEnergy(G4double e)
void SetMomentumAndUpdate(G4ParticleMomentum mom)
void SetNumberOfSecondaries(G4int totSecondaries)
G4ParticleChange aParticleChange
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void SetMass(G4double mas)
void AddSecondary(G4Track *aSecondary)
size_t GetNumberOfElements() const
void ProposeTrackStatus(G4TrackStatus status)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void PrintInfo(const G4ParticleDefinition *)