73{
74
75 theParticleChange->Clear();
78
79
80
81
82
83
84
85
86
87 const G4double energyThresholdForCharmAndBottomHadrons = 100.0*CLHEP::MeV;
88 if ( thePrimary.
GetKineticEnergy() < energyThresholdForCharmAndBottomHadrons &&
93 theParticleChange->SetStatusChange(
isAlive );
96 return theParticleChange;
97 }
98
99
100
101
102
103
104
105
106
107
108
109
110 const G4double energyThresholdForHyperNuclei = 100.0*CLHEP::MeV;
113 theParticleChange->SetStatusChange(
isAlive );
116 return theParticleChange;
117 }
118
119
120
122
123 if ( theQuasielastic )
124 {
125 if ( theQuasielastic->GetFraction(theNucleus, aPart) >
G4UniformRand() )
126 {
127
128 G4KineticTrackVector *result= theQuasielastic->Scatter(theNucleus, aPart);
129
130 if (result)
131 {
132 for(auto & ptr : *result)
133 {
134 G4DynamicParticle * aNew =
135 new G4DynamicParticle(ptr->GetDefinition(),
136 ptr->Get4Momentum().e(),
137 ptr->Get4Momentum().vect());
138 theParticleChange->AddSecondary(aNew, ptr->GetCreatorModelID());
139 delete ptr;
140 }
141 delete result;
142 }
143 else
144 {
145 theParticleChange->SetStatusChange(
isAlive);
148 }
149 return theParticleChange;
150 }
151 }
152
153
154 G4KineticTrackVector * theInitialResult =
155 theHighEnergyGenerator->Scatter(theNucleus, aPart);
156
157
158 for ( auto & ptr : *theInitialResult ) {
159 ptr->SetCreatorModelID( theStringModelID );
160 }
161
162
163 #ifdef DEBUG_initial_result
166 for(auto & ptr : *theInitialResult)
167 {
168
169
170 E_out += ptr->Get4Momentum().e();
171 }
173 G4double init_E=aPart.Get4Momentum().e();
174
175
176 const std::vector<G4Nucleon> & thy = theHighEnergyGenerator->GetWoundedNucleus()->GetNucleons();
177
178 G4int resZ(0),resA(0);
180 for(auto & nuc : thy)
181 {
182 if(
nuc.AreYouHit()) {
183 ++resA;
185 ++resZ;
186 delta_m += CLHEP::proton_mass_c2;
187 } else {
188 delta_m += CLHEP::neutron_mass_c2;
189 }
190 }
191 }
192
196 }
197 G4double E_excit=init_mass + init_E - final_mass - E_out;
198 G4cout <<
" Corrected delta mass " << init_mass - final_mass - delta_m <<
G4endl;
199 G4cout <<
"initial E, mass = " << init_E <<
", " << init_mass <<
G4endl;
200 G4cout <<
" final E, mass = " << E_out <<
", " << final_mass <<
" excitation_E " << E_excit <<
G4endl;
201 #endif
202
204
205 G4V3DNucleus* theProjectileNucleus = theHighEnergyGenerator->GetProjectileNucleus();
206 if(theProjectileNucleus == nullptr)
207 {
209 const std::vector<G4Nucleon>& they = theHighEnergyGenerator->GetWoundedNucleus()->GetNucleons();
210 for(auto & nuc : they)
211 {
212 if(
nuc.AreYouHit()) ++hitCount;
213 }
214 if(hitCount != theHighEnergyGenerator->GetWoundedNucleus()->GetMassNumber() )
215 {
216 theTransport->SetPrimaryProjectile(thePrimary);
217 theTransportResult =
218 theTransport->Propagate(theInitialResult,
219 theHighEnergyGenerator->GetWoundedNucleus());
220 if ( !theTransportResult ) {
221 G4cout <<
"G4TheoFSGenerator: null ptr from transport propagate " <<
G4endl;
222 throw G4HadronicException(__FILE__, __LINE__, "Null ptr from transport propagate");
223 }
224 }
225 else
226 {
227 theTransportResult = theDecay.Propagate(theInitialResult,
228 theHighEnergyGenerator->GetWoundedNucleus());
229 if ( theTransportResult == nullptr ) {
230 G4cout <<
"G4TheoFSGenerator: null ptr from decay propagate " <<
G4endl;
231 throw G4HadronicException(__FILE__, __LINE__, "Null ptr from decay propagate");
232 }
233 }
234 }
235 else
236 {
237 theTransport->SetPrimaryProjectile(thePrimary);
238 theTransportResult = theTransport->PropagateNuclNucl(theInitialResult,
239 theHighEnergyGenerator->GetWoundedNucleus(),
240 theProjectileNucleus);
241 if ( !theTransportResult ) {
242 G4cout <<
"G4TheoFSGenerator: null ptr from transport propagate " <<
G4endl;
243 throw G4HadronicException(__FILE__, __LINE__, "Null ptr from transport propagate");
244 }
245 }
246
247
248
249
250
252 if(nullptr == theCosmicCoalescence) {
253 theCosmicCoalescence = (G4CRCoalescence*)
255 if(nullptr == theCosmicCoalescence) {
256 theCosmicCoalescence = new G4CRCoalescence();
257 }
258 }
259 theCosmicCoalescence->SetP0Coalescence( thePrimary, theHighEnergyGenerator->GetModelName() );
260 theCosmicCoalescence->GenerateDeuterons( theTransportResult );
261 }
262
263
264 for(auto & ptr : *theTransportResult)
265 {
266 G4DynamicParticle * aNewDP =
267 new G4DynamicParticle(ptr->GetDefinition(),
268 ptr->GetTotalEnergy(),
269 ptr->GetMomentum());
270 G4HadSecondary aNew = G4HadSecondary(aNewDP);
271 G4double time = std::max(ptr->GetFormationTime(), 0.0);
272 aNew.
SetTime(timePrimary + time);
276 theParticleChange->AddSecondary(aNew);
277 delete ptr;
278 }
279
280
281 delete theTransportResult;
282
283
284 return theParticleChange;
285}
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4LorentzVector & Get4Momentum() const
G4double GetGlobalTime() const
void SetParentResonanceDef(const G4ParticleDefinition *parentDef)
void SetTime(G4double aT)
void SetCreatorModelID(G4int id)
void SetParentResonanceID(const G4int parentID)
G4HadronicInteraction * FindModel(const G4String &name)
static G4HadronicInteractionRegistry * Instance()
static G4HadronicParameters * Instance()
G4double GetIonMass(G4int Z, G4int A, G4int nL=0, G4int lvl=0) const
G4int GetQuarkContent(G4int flavor) const
G4bool IsHypernucleus() const
G4int GetAntiQuarkContent(G4int flavor) const
G4IonTable * GetIonTable() const
static G4ParticleTable * GetParticleTable()
static G4Proton * Proton()