Main method to apply the INCL physics model.
164{
171
172
173 if((isIonTrack && (trackZ<=0 || trackA<=trackZ)) || (nucleusA>1 && (nucleusZ<=0 || nucleusA<=nucleusZ))) {
178 return &theResult;
179 }
180
181
182 if(trackA<=1 && nucleusA<=1) {
183 return theBackupModelNucleon->
ApplyYourself(aTrack, theNucleus);
184 }
185
186
187
189 if(trackA > theMaxProjMassINCL && nucleusA > theMaxProjMassINCL) {
190 if(!complainedAboutBackupModel) {
191 complainedAboutBackupModel = true;
192 std::stringstream ss;
193 ss << "INCL++ refuses to handle reactions between nuclei with A>"
194 << theMaxProjMassINCL
195 << ". A backup model ("
197 << ") will be used instead.";
199 }
201 }
202
203
207 && trackKinE < cascadeMinEnergyPerNucleon) {
208 if(!complainedAboutPreCompound) {
209 complainedAboutPreCompound = true;
210 std::stringstream ss;
211 ss << "INCL++ refuses to handle nucleon-induced reactions below "
212 << cascadeMinEnergyPerNucleon / MeV
213 << " MeV. A PreCoumpound model ("
215 << ") will be used instead.";
217 }
218 return thePreCompoundModel->
ApplyYourself(aTrack, theNucleus);
219 }
220
221
224 const G4double theTrackEnergy = trackKinE + theTrackMass;
225 const G4double theTrackMomentumAbs2 = theTrackEnergy*theTrackEnergy - theTrackMass*theTrackMass;
226 const G4double theTrackMomentumAbs = ((theTrackMomentumAbs2>0.0) ? std::sqrt(theTrackMomentumAbs2) : 0.0);
228
231 fourMomentumIn.
setE(theTrackEnergy + theNucleusMass);
232 fourMomentumIn.
setVect(theTrackMomentum);
233
234
235 const G4bool inverseKinematics = AccurateProjectile(aTrack, theNucleus);
236
237
241 G4Nucleus *theTargetNucleus = &theNucleus;
242 if(inverseKinematics) {
245
246 if(oldProjectileDef != 0 && oldTargetDef != 0) {
249
250 if(newTargetA > 0 && newTargetZ > 0) {
251
252 theTargetNucleus =
new G4Nucleus(newTargetA, newTargetZ);
255 G4DynamicParticle swappedProjectileParticle(oldTargetDef, (*toInverseKinematics) * theProjectile4Momentum);
257 } else {
258 G4String message =
"badly defined target after swapping. Falling back to normal (non-swapped) mode.";
261 }
262 } else {
263 G4String message =
"oldProjectileDef or oldTargetDef was null";
266 }
267 }
268
269
270
271
272
274
275
276
277
278
279
280
281
282
288
289
290
291
292
293
294
295
298
299 std::list<G4Fragment> remnants;
300
301 const G4int maxTries = 200;
303
304
305
307 do {
309 const G4double kineticEnergy = toINCLKineticEnergy(*aProjectileTrack);
310
311
313
315
317 if(eventIsOK) {
318
319
320
321 if(inverseKinematics) {
323 }
324
326
331
337 if(p != 0) {
339
340
341 momentum *= toLabFrame;
342
343 if(inverseKinematics) {
344 momentum *= *toDirectKinematics;
346 }
347
348
350 fourMomentumOut += momentum;
352
353 } else {
354 G4String message =
"the model produced a particle that couldn't be converted to Geant4 particle.";
356 }
357 }
358
362
368 eventInfo.
jxRem[i]*hbar_Planck,
369 eventInfo.
jyRem[i]*hbar_Planck,
370 eventInfo.
jzRem[i]*hbar_Planck
371 );
374 const G4double scaling = remnant4MomentumScaling(nuclearMass,
375 kinE,
376 px, py, pz);
378 nuclearMass + kinE);
379 if(std::abs(scaling - 1.0) > 0.01) {
380 std::stringstream ss;
381 ss << "momentum scaling = " << scaling
382 << "\n Lorentz vector = " << fourMomentum
383 <<
")\n A = " <<
A <<
", Z = " << Z
384 << "\n E* = " << excitationE << ", nuclearMass = " << nuclearMass
385 <<
"\n remnant i=" << i <<
", nRemnants=" << eventInfo.
nRemnants
389 << ", in " << (inverseKinematics ? "inverse" : "direct") << " kinematics.";
391 }
392
393
394 fourMomentum *= toLabFrame;
395 spin *= toLabFrame3;
396
397 if(inverseKinematics) {
398 fourMomentum *= *toDirectKinematics;
399 fourMomentum.setVect(-fourMomentum.vect());
400 }
401
402 fourMomentumOut += fourMomentum;
404 remnant.SetAngularMomentum(spin);
405 if(dumpRemnantInfo) {
406 G4cerr <<
"G4INCLXX_DUMP_REMNANT: " << remnant <<
" spin: " << spin <<
G4endl;
407 }
408 remnants.push_back(remnant);
409 }
410
411
412 const G4LorentzVector violation4Momentum = fourMomentumOut - fourMomentumIn;
413 const G4double energyViolation = std::abs(violation4Momentum.
e());
414 const G4double momentumViolation = violation4Momentum.
rho();
416 std::stringstream ss;
417 ss << "energy conservation violated by " << energyViolation/MeV << " MeV in "
420 << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics. Will resample.";
422 eventIsOK = false;
424 for(
G4int j=0; j<nSecondaries; ++j)
428 remnants.clear();
430 std::stringstream ss;
431 ss << "momentum conservation violated by " << momentumViolation/MeV << " MeV in "
434 << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics. Will resample.";
436 eventIsOK = false;
438 for(
G4int j=0; j<nSecondaries; ++j)
442 remnants.clear();
443 }
444 }
445 nTries++;
446 } while(!eventIsOK && nTries < maxTries);
447
448
449 if(inverseKinematics) {
450 delete aProjectileTrack;
451 delete theTargetNucleus;
452 delete toInverseKinematics;
453 delete toDirectKinematics;
454 }
455
456 if(!eventIsOK) {
457 std::stringstream ss;
458 ss << "maximum number of tries exceeded for the proposed "
460 <<
" + " << theIonTable->
GetIonName(nucleusZ, nucleusA, 0)
461 << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics.";
466 return &theResult;
467 }
468
469
470
472 for(std::list<G4Fragment>::iterator i = remnants.begin();
473 i != remnants.end(); ++i) {
475
476 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
477 fragment != deExcitationResult->end(); ++fragment) {
479 if(def != 0) {
482 }
483 }
484
485 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
486 fragment != deExcitationResult->end(); ++fragment) {
487 delete (*fragment);
488 }
489 deExcitationResult->clear();
490 delete deExcitationResult;
491 }
492 }
493
494 remnants.clear();
495
496 if((theTally = theInterfaceStore->
GetTally()))
497 theTally->
Tally(aTrack, theNucleus, theResult);
498
499 return &theResult;
500}
double A(double temperature)
CLHEP::HepLorentzRotation G4LorentzRotation
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4GLOB_DLL std::ostream G4cerr
HepLorentzRotation inverse() const
void setVect(const Hep3Vector &)
HepRotation inverse() const
HepRotation & rotateZ(double delta)
HepRotation & rotateY(double delta)
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
static G4GenericIon * GenericIon()
void SetStatusChange(G4HadFinalStateStatus aS)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
std::size_t GetNumberOfSecondaries() const
void SetEnergyChange(G4double anEnergy)
G4HadSecondary * GetSecondary(size_t i)
void SetMomentumChange(const G4ThreeVector &aV)
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4LorentzVector & Get4Momentum() const
G4DynamicParticle * GetParticle()
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
const G4String & GetModelName() const
G4int GetMaxProjMassINCL() const
Getter for theMaxProjMassINCL.
G4INCL::INCL * GetINCLModel()
Get the cached INCL model engine.
G4double GetConservationTolerance() const
Getter for conservationTolerance.
G4double GetCascadeMinEnergyPerNucleon() const
Getter for cascadeMinEnergyPerNucleon.
G4INCLXXVInterfaceTally * GetTally() const
Getter for the interface tally.
virtual void Tally(G4HadProjectile const &aTrack, G4Nucleus const &theNucleus, G4HadFinalState const &result)=0
const EventInfo & processEvent()
const G4String & GetIonName(G4int Z, G4int A, G4int lvl=0) const
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
static G4Neutron * NeutronDefinition()
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4int GetAtomicNumber() const
const G4String & GetParticleType() const
G4double GetPDGMass() const
G4int GetAtomicMass() const
G4double GetPDGCharge() const
const G4String & GetParticleName() const
static G4Proton * ProtonDefinition()
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)=0
Short_t nRemnants
Number of remnants.
Bool_t transparent
True if the event is transparent.
Short_t A[maxSizeParticles]
Particle mass number.
Float_t EKinRem[maxSizeRemnants]
Remnant kinetic energy [MeV].
Float_t jzRem[maxSizeRemnants]
Remnant angular momentum, z component [ ].
Float_t EStarRem[maxSizeRemnants]
Remnant excitation energy [MeV].
Short_t Z[maxSizeParticles]
Particle charge number.
Float_t EKin[maxSizeParticles]
Particle kinetic energy [MeV].
Short_t nParticles
Number of particles in the final state.
Float_t jxRem[maxSizeRemnants]
Remnant angular momentum, x component [ ].
Float_t px[maxSizeParticles]
Particle momentum, x component [MeV/c].
Float_t pxRem[maxSizeRemnants]
Remnant momentum, x component [MeV/c].
Float_t pyRem[maxSizeRemnants]
Remnant momentum, y component [MeV/c].
Float_t pzRem[maxSizeRemnants]
Remnant momentum, z component [MeV/c].
Float_t pz[maxSizeParticles]
Particle momentum, z component [MeV/c].
Int_t PDGCode[maxSizeParticles]
PDG numbering of the particles.
Float_t jyRem[maxSizeRemnants]
Remnant angular momentum, y component [ ].
Float_t py[maxSizeParticles]
Particle momentum, y component [MeV/c].
Short_t ARem[maxSizeRemnants]
Remnant mass number.
Short_t ZRem[maxSizeRemnants]
Remnant charge number.