86{
88
89
91
92
97 G4int numberOfHoles = 0;
103
104
105#ifdef exactExcitationEnergy
107#endif
108 G4KineticTrackVector::iterator iter;
109 for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
110 {
112 G4double e = (*iter)->Get4Momentum().e();
113 G4double mass = (*iter)->Get4Momentum().mag();
115 if((part != proton && part != neutron) ||
116 (e > mass + CaptureThreshold) ||
117 ((*iter)->GetPosition().mag() > R)) {
121 theTotalResult->push_back(theNew);
122#ifdef exactExcitationEnergy
123 secondary4Momemtum += (*iter)->Get4Momentum();
124#endif
125 } else {
126
127
128 ++anA;
129 ++numberOfEx;
131 aZ += Z;
132 numberOfCh += Z;
133 captured3Momentum += mom;
134 exEnergy += (e - mass);
135 }
136 delete (*iter);
137 }
138 delete theSecondaries;
139
140
143 while(theCurrentNucleon) {
145 ++numberOfHoles;
146 ++numberOfEx;
147 --anA;
150
152 }
154 }
155 exciton3Momentum = captured3Momentum - wounded3Momentum;
156
157 if(anA>0 && aZ>0) {
159#ifdef exactExcitationEnergy
160
164 G4double Efinal = fMass + secondary4Momemtum.e();
165 if ( (Einitial - Efinal) > 0 ) {
166
167
168 exEnergy=Einitial - Efinal;
169 }
170 else {
171
172
173 exEnergy=0.;
174 }
175#endif
176 fMass += exEnergy;
178 #ifdef G4GPI_debug_excitation
179 G4cout <<
"momentum balance init/final " << balance <<
" value " << balance.
mag() <<
G4endl
180 <<
"primary / secondaries "<< primary->
Get4Momentum() <<
" / "
181 << secondary4Momemtum << " captured/wounded: " << captured3Momentum << " / " << wounded3Momentum
182 <<
" exciton " << exciton3Momentum <<
G4endl
183 << secondary4Momemtum.vect() + exciton3Momentum <<
G4endl;
184 #endif
185#ifdef exactExcitationEnergy
187#else
189 std::sqrt(exciton3Momentum.mag2() + fMass*fMass));
190#endif
191 if ( exEnergy > 0.0 ) {
192 G4Fragment anInitialState(anA, aZ, exciton4Momentum);
193 anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
194 anInitialState.SetNumberOfCharged(numberOfCh);
195 anInitialState.SetNumberOfHoles(numberOfHoles);
196
198
199 theTotalResult->insert(theTotalResult->end(),aPrecoResult->begin(),aPrecoResult->end() );
200 delete aPrecoResult;
201
202 } else {
203
205 if (anA == 1 && aZ == 0) {
207 } else if (anA == 1 && aZ == 1) {
209 } else if (anA == 2 && aZ == 1) {
211 } else if (anA == 3 && aZ == 1) {
213 } else if (anA == 3 && aZ == 2) {
215 } else if (anA == 4 && aZ == 2) {
217 } else {
218 theKindOfFragment =
220 }
221 if (theKindOfFragment != 0) {
225
226 theTotalResult->push_back(theNew);
227 }
228 }
229 }
230
231 return theTotalResult;
232}
std::vector< G4ReactionProduct * > G4ReactionProductVector
CLHEP::Hep3Vector G4ThreeVector
static G4Alpha * AlphaDefinition()
static G4Deuteron * DeuteronDefinition()
const G4LorentzVector & Get4Momentum() const
static G4He3 * He3Definition()
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int J=0)
static G4Neutron * NeutronDefinition()
static G4double GetNuclearMass(const G4double A, const G4double Z)
virtual const G4LorentzVector & Get4Momentum() const
G4double GetBindingEnergy() const
virtual G4ParticleDefinition * GetDefinition() const
G4double GetPDGCharge() const
static G4ParticleTable * GetParticleTable()
G4IonTable * GetIonTable()
static G4Proton * ProtonDefinition()
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetTotalEnergy(const G4double en)
static G4Triton * TritonDefinition()
virtual G4Nucleon * GetNextNucleon()=0
virtual G4int GetCharge()=0
virtual G4bool StartLoop()=0
virtual G4double GetNuclearRadius()=0
virtual G4int GetMassNumber()=0
G4VPreCompoundModel * theDeExcitation
const G4HadProjectile * GetPrimaryProjectile() const
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)=0