55 {
56
57 G4bool isNN = theNucleus->isNucleusNucleusCollision();
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
87 if(isNN) {
88
89 ProjectileRemnant * const projectileRemnant = theNucleus->getProjectileRemnant();
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122 G4double theProjectileExcitationEnergy = 0;
123 G4double theProjectileEffectiveMass =0;
124 if (theParticle->isNucleonorLambda()){
125 theProjectileExcitationEnergy = (projectileRemnant->getA()-theParticle->getA()>1) ? (projectileRemnant->computeExcitationEnergyExcept(theParticle->getID())) : 0.;
126 theProjectileEffectiveMass =
127 ParticleTable::getTableMass(projectileRemnant->getA() - theParticle->getA(), projectileRemnant->getZ() - theParticle->getZ(), projectileRemnant->getS() - theParticle->getS())
128 + theProjectileExcitationEnergy;
129 }
130 else if (theParticle->isAntiNucleon()){
131 theProjectileExcitationEnergy = (projectileRemnant->getA() -theParticle->getA()<-1) ? (projectileRemnant->computeExcitationEnergyExcept(theParticle->getID())) : 0;
132 theProjectileEffectiveMass =
133 ParticleTable::getTableMass(-(projectileRemnant->getA() - theParticle->getA()), -(projectileRemnant->getZ() - theParticle->getZ()), projectileRemnant->getS() - theParticle->getS())
134 + theProjectileExcitationEnergy;
135 }
136
137
138
139
140 const ThreeVector &theProjectileMomentum = projectileRemnant->getMomentum() - theParticle->getMomentum();
141 const G4double theProjectileEnergy = std::sqrt(theProjectileMomentum.mag2() + theProjectileEffectiveMass*theProjectileEffectiveMass);
142 const G4double theProjectileCorrection = theProjectileEnergy - (projectileRemnant->getEnergy() - theParticle->getEnergy());
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165 theCorrection = theParticle->getEmissionQValueCorrection(
166 theNucleus->getA() + theParticle->getA(),
167 theNucleus->getZ() + theParticle->getZ(),
168 theNucleus->getS() + theParticle->getS())
169 + theParticle->getTableMass() - theParticle->getINCLMass()
170 + theProjectileCorrection;
171
172
173
174 projectileRemnant->removeParticle(theParticle, theProjectileCorrection);
175 } else {
176 const G4int ACN = theNucleus->getA() + theParticle->getA();
177 const G4int ZCN = theNucleus->getZ() + theParticle->getZ();
178 const G4int SCN = theNucleus->getS() + theParticle->getS();
179
180 theCorrection = theParticle->getEmissionQValueCorrection(ACN,ZCN,SCN);
181 INCL_DEBUG(
"The following Particle enters with correction " << theCorrection <<
'\n'
182 << theParticle->print() << '\n');
183 }
184
185 const G4double energyBefore = theParticle->getEnergy() - theCorrection;
187 if(isNN && theParticle->isAntiNucleon() && (theParticle->getEnergy() - theCorrection <=theParticle->
getINCLMass()) ){
188 success =true;
189 G4double energyInside = theParticle->getEnergy() + theNucleus->getPotential()->computePotentialEnergy(theParticle) - theCorrection;
190 theParticle->setEnergy(energyInside);
191 theParticle->setPotentialEnergy(theNucleus->getPotential()->computePotentialEnergy(theParticle));
192 theParticle->setMomentum(theParticle->getMomentum());
193 theParticle->adjustMomentumFromEnergy();
194 }
195 else{
196 success = particleEnters(theCorrection);
197 }
198 fs->addEnteringParticle(theParticle);
199
200 if(!success) {
201 fs->makeParticleBelowZero();
202 } else if(theParticle->isNucleonorLambda() &&
203 theParticle->getKineticEnergy()<theNucleus->getPotential()->getFermiEnergy(theParticle)) {
204
205
206 fs->makeParticleBelowFermi();
207 } else if(theParticle->isKaon()) theNucleus->setNumberOfKaon(theNucleus->getNumberOfKaon()+1);
208
209 fs->setTotalEnergyBeforeInteraction(energyBefore);
210 }
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
G4double getINCLMass(const G4int A, const G4int Z, const G4int S)
Get INCL nuclear mass (in MeV/c^2).