| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35
36 #include "globals.hh"
37
38 #include <cmath>
39
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4INCLXXInterface.hh"
43 #include "G4GenericIon.hh"
44 #include "G4INCLCascade.hh"
45 #include "G4ReactionProductVector.hh"
46 #include "G4ReactionProduct.hh"
47 #include "G4HadSecondary.hh"
48 #include "G4ParticleTable.hh"
49 #include "G4INCLXXInterfaceStore.hh"
50 #include "G4INCLXXVInterfaceTally.hh"
51 #include "G4String.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4HadronicInteractionRegistry.hh"
55 #include "G4INCLVersion.hh"
56 #include "G4VEvaporation.hh"
57 #include "G4VEvaporationChannel.hh"
58 #include "G4CompetitiveFission.hh"
59 #include "G4FissionLevelDensityParameterINCLXX
60 #include "G4PhysicsModelCatalog.hh"
61
62 #include "G4HyperNucleiProperties.hh"
63 #include "G4HyperTriton.hh"
64 #include "G4HyperH4.hh"
65 #include "G4HyperAlpha.hh"
66 #include "G4DoubleHyperH4.hh"
67 #include "G4DoubleHyperDoubleNeutron.hh"
68 #include "G4HyperHe5.hh"
69
70 G4INCLXXInterface::G4INCLXXInterface(G4VPreCom
71 G4VIntraNuclearTransportModel(G4INCLXXInterf
72 theINCLModel(NULL),
73 thePreCompoundModel(aPreCompound),
74 theInterfaceStore(G4INCLXXInterfaceStore::Ge
75 theTally(NULL),
76 complainedAboutBackupModel(false),
77 complainedAboutPreCompound(false),
78 theIonTable(G4IonTable::GetIonTable()),
79 theINCLXXLevelDensity(NULL),
80 theINCLXXFissionProbability(NULL),
81 secID(-1)
82 {
83 if(!thePreCompoundModel) {
84 G4HadronicInteraction* p =
85 G4HadronicInteractionRegistry::Instance(
86 thePreCompoundModel = static_cast<G4VPreCo
87 if(!thePreCompoundModel) { thePreCompoundM
88 }
89
90 // Use the environment variable G4INCLXX_NO_
91 if(std::getenv("G4INCLXX_NO_DE_EXCITATION"))
92 G4String message = "de-excitation is compl
93 theInterfaceStore->EmitWarning(message);
94 theDeExcitation = 0;
95 } else {
96 G4HadronicInteraction* p =
97 G4HadronicInteractionRegistry::Instance(
98 theDeExcitation = static_cast<G4VPreCompou
99 if(!theDeExcitation) { theDeExcitation = n
100
101 // set the fission parameters for G4Excita
102 G4VEvaporationChannel * const theFissionCh
103 theDeExcitation->GetExcitationHandler()-
104 G4CompetitiveFission * const theFissionCha
105 if(theFissionChannelCast) {
106 theINCLXXLevelDensity = new G4FissionLev
107 theFissionChannelCast->SetLevelDensityPa
108 theINCLXXFissionProbability = new G4Fiss
109 theINCLXXFissionProbability->SetFissionL
110 theFissionChannelCast->SetEmissionStrate
111 theInterfaceStore->EmitBigWarning("INCL+
112 } else {
113 theInterfaceStore->EmitBigWarning("INCL+
114 }
115 }
116
117 // use the envvar G4INCLXX_DUMP_REMNANT to d
118 // remnants on stdout
119 if(std::getenv("G4INCLXX_DUMP_REMNANT"))
120 dumpRemnantInfo = true;
121 else
122 dumpRemnantInfo = false;
123
124 theBackupModel = new G4BinaryLightIonReactio
125 theBackupModelNucleon = new G4BinaryCascade;
126 secID = G4PhysicsModelCatalog::GetModelID( "
127 }
128
129 G4INCLXXInterface::~G4INCLXXInterface()
130 {
131 delete theINCLXXLevelDensity;
132 delete theINCLXXFissionProbability;
133 }
134
135 G4bool G4INCLXXInterface::AccurateProjectile(c
136 // Use direct kinematics if the projectile i
137 const G4ParticleDefinition *projectileDef =
138 if(std::abs(projectileDef->GetBaryonNumber()
139 return false;
140
141 // Here all projectiles should be nuclei
142 const G4int pA = projectileDef->GetAtomicMas
143 if(pA<=0) {
144 std::stringstream ss;
145 ss << "the model does not know how to hand
146 << projectileDef->GetParticleName() << "
147 << theNucleus.GetZ_asInt() << ", A=" <<
148 theInterfaceStore->EmitBigWarning(ss.str()
149 return true;
150 }
151
152 // If either nucleus is a LCP (A<=4), run th
153 const G4int tA = theNucleus.GetA_asInt();
154 if(tA<=4 || pA<=4) {
155 if(pA<tA)
156 return false;
157 else
158 return true;
159 }
160
161 // If one of the nuclei is heavier than theM
162 // as light on heavy.
163 // Note that here we are sure that either th
164 // smaller than theMaxProjMass; otherwise th
165 // called.
166 const G4int theMaxProjMassINCL = theInterfac
167 if(pA > theMaxProjMassINCL)
168 return true;
169 else if(tA > theMaxProjMassINCL)
170 return false;
171 else
172 // In all other cases, use the global sett
173 return theInterfaceStore->GetAccurateProje
174 }
175
176 G4HadFinalState* G4INCLXXInterface::ApplyYours
177 {
178 G4ParticleDefinition const * const trackDefi
179 const G4bool isIonTrack = trackDefinition->G
180 const G4int trackA = trackDefinition->GetAto
181 const G4int trackZ = (G4int) trackDefinition
182 const G4int trackL = trackDefinition->GetNum
183 const G4int nucleusA = theNucleus.GetA_asInt
184 const G4int nucleusZ = theNucleus.GetZ_asInt
185
186 // For reactions induced by weird projectile
187 if((isIonTrack && ((trackZ<=0 && trackL==0)
188 (nucleusA>1 && (nucleusZ<=
189 theResult.Clear();
190 theResult.SetStatusChange(isAlive);
191 theResult.SetEnergyChange(aTrack.GetKineti
192 theResult.SetMomentumChange(aTrack.Get4Mom
193 return &theResult;
194 }
195
196 // For reactions on nucleons, use the backup
197 // except for anti_proton projectile (in thi
198 if(trackA<=1 && nucleusA<=1 && (trackZ>=0 ||
199 return theBackupModelNucleon->ApplyYoursel
200 }
201
202 // For systems heavier than theMaxProjMassIN
203 // BIC)
204 const G4int theMaxProjMassINCL = theInterfac
205 if(trackA > theMaxProjMassINCL && nucleusA >
206 if(!complainedAboutBackupModel) {
207 complainedAboutBackupModel = true;
208 std::stringstream ss;
209 ss << "INCL++ refuses to handle reaction
210 << theMaxProjMassINCL
211 << ". A backup model ("
212 << theBackupModel->GetModelName()
213 << ") will be used instead.";
214 theInterfaceStore->EmitBigWarning(ss.str
215 }
216 return theBackupModel->ApplyYourself(aTrac
217 }
218
219 // For energies lower than cascadeMinEnergyP
220 const G4double cascadeMinEnergyPerNucleon =
221 const G4double trackKinE = aTrack.GetKinetic
222 if((trackDefinition==G4Neutron::NeutronDefin
223 && trackKinE < cascadeMinEnergyPerNucleo
224 if(!complainedAboutPreCompound) {
225 complainedAboutPreCompound = true;
226 std::stringstream ss;
227 ss << "INCL++ refuses to handle nucleon-
228 << cascadeMinEnergyPerNucleon / MeV
229 << " MeV. A PreCoumpound model ("
230 << thePreCompoundModel->GetModelName()
231 << ") will be used instead.";
232 theInterfaceStore->EmitBigWarning(ss.str
233 }
234 return thePreCompoundModel->ApplyYourself(
235 }
236
237 // Calculate the total four-momentum in the
238 const G4double theNucleusMass = theIonTable-
239 const G4double theTrackMass = trackDefinitio
240 const G4double theTrackEnergy = trackKinE +
241 const G4double theTrackMomentumAbs2 = theTra
242 const G4double theTrackMomentumAbs = ((theTr
243 const G4ThreeVector theTrackMomentum = aTrac
244
245 G4LorentzVector goodTrack4Momentum(theTrackM
246 G4LorentzVector fourMomentumIn;
247 fourMomentumIn.setE(theTrackEnergy + theNucl
248 fourMomentumIn.setVect(theTrackMomentum);
249
250 // Check if inverse kinematics should be use
251 const G4bool inverseKinematics = AccuratePro
252
253 // If we are running in inverse kinematics,
254 G4LorentzRotation *toInverseKinematics = NUL
255 G4LorentzRotation *toDirectKinematics = NULL
256 G4HadProjectile const *aProjectileTrack = &a
257 G4Nucleus *theTargetNucleus = &theNucleus;
258 if(inverseKinematics) {
259 G4ParticleDefinition *oldTargetDef = theIo
260 const G4ParticleDefinition *oldProjectileD
261
262 if(oldProjectileDef != 0 && oldTargetDef !
263 const G4int newTargetA = oldProjectileDe
264 const G4int newTargetZ = oldProjectileDe
265 const G4int newTargetL = oldProjectileDe
266 if(newTargetA > 0 && newTargetZ > 0) {
267 // This should give us the same energy
268 theTargetNucleus = new G4Nucleus(newTa
269 toInverseKinematics = new G4LorentzRot
270 G4LorentzVector theProjectile4Momentum
271 G4DynamicParticle swappedProjectilePar
272 aProjectileTrack = new G4HadProjectile
273 } else {
274 G4String message = "badly defined targ
275 theInterfaceStore->EmitWarning(message
276 toInverseKinematics = new G4LorentzRot
277 }
278 } else {
279 G4String message = "oldProjectileDef or
280 theInterfaceStore->EmitWarning(message);
281 toInverseKinematics = new G4LorentzRotat
282 }
283 }
284
285 // INCL assumes the projectile particle is g
286 // the Z-axis. Here we construct proper rota
287 // momentum vectors of the outcoming particl
288 // coordinate system.
289 G4LorentzVector projectileMomentum = aProjec
290
291 // INCL++ assumes that the projectile is goi
292 // the z-axis. In principle, if the coordina
293 // hadronic framework is defined differently
294 // transform the INCL++ reaction products to
295 // frame. Please note that it isn't necessar
296 // transform to the projectile because when
297 // projectile particle; \see{toINCLKineticEn
298 // projectile energy (direction is simply as
299 G4RotationMatrix toZ;
300 toZ.rotateZ(-projectileMomentum.phi());
301 toZ.rotateY(-projectileMomentum.theta());
302 G4RotationMatrix toLabFrame3 = toZ.inverse()
303 G4LorentzRotation toLabFrame(toLabFrame3);
304 // However, it turns out that the projectile
305 // hadronic framework is already going in th
306 // z-axis so this rotation is actually unnec
307 // toLabFrame turn out to be unit matrices a
308 // uncommenting the folowing two lines:
309 // G4cout <<"toZ = " << toZ << G4endl;
310 // G4cout <<"toLabFrame = " << toLabFrame <<
311
312 theResult.Clear(); // Make sure the output d
313 theResult.SetStatusChange(stopAndKill);
314
315 std::list<G4Fragment> remnants;
316
317 const G4int maxTries = 200;
318 G4int nTries = 0;
319 // INCL can generate transparent events. How
320 // only in the standalone code. In Geant4 we
321 // produce a valid cascade.
322 G4bool eventIsOK = false;
323 do {
324 const G4INCL::ParticleSpecies theSpecies =
325 const G4double kineticEnergy = toINCLKinet
326
327 // The INCL model will be created at the f
328 theINCLModel = G4INCLXXInterfaceStore::Get
329
330 const G4INCL::EventInfo eventInfo = theINC
331 theTargetNucleus->GetA_asIn
332 theTargetNucleus->GetZ_asIn
333 -theTargetNucleus->GetL());
334 // eventIsOK = !eventInfo.transparent &
335 eventIsOK = !eventInfo.transparent;
336 if(eventIsOK) {
337
338 // If the collision was run in inverse k
339 // all the reaction products
340 if(inverseKinematics) {
341 toDirectKinematics = new G4LorentzRota
342 }
343
344 G4LorentzVector fourMomentumOut;
345
346 for(G4int i = 0; i < eventInfo.nParticle
347 G4int A = eventInfo.A[i];
348 G4int Z = eventInfo.Z[i];
349 G4int S = eventInfo.S[i]; // Strangen
350 G4int PDGCode = eventInfo.PDGCode[i];
351 G4double kinE = eventInfo.EKin[i];
352 G4double px = eventInfo.px[i];
353 G4double py = eventInfo.py[i];
354 G4double pz = eventInfo.pz[i];
355 G4DynamicParticle *p = toG4Particle(A,
356 if(p != 0) {
357 G4LorentzVector momentum = p->Get4Mo
358
359 // Apply the toLabFrame rotation
360 momentum *= toLabFrame;
361 // Apply the inverse-kinematics boos
362 if(inverseKinematics) {
363 momentum *= *toDirectKinematics;
364 momentum.setVect(-momentum.vect())
365 }
366
367 // Set the four-momentum of the reaction p
368 p->Set4Momentum(momentum);
369 fourMomentumOut += momentum;
370
371 // Propagate the particle's parent resonan
372 G4HadSecondary secondary(p, 1.0, sec
373 G4ParticleDefinition* parentResonanc
374 if ( eventInfo.parentResonancePDGCod
375 parentResonanceDef = G4ParticleTab
376 }
377 secondary.SetParentResonanceDef(pare
378 secondary.SetParentResonanceID(event
379
380 theResult.AddSecondary(secondary);
381
382 } else {
383 G4String message = "the model produc
384 theInterfaceStore->EmitWarning(messa
385 }
386 }
387
388 for(G4int i = 0; i < eventInfo.nRemnants
389 const G4int A = eventInfo.ARem[i];
390 const G4int Z = eventInfo.ZRem[i];
391 const G4int S = eventInfo.SRem[i];
392 // Check that the remnant is a physica
393 if(( Z == 0 && S == 0 && A > 1 ) |
394 ( Z == 0 && S != 0 && A < 4 ) |
395 ( Z != 0 && S != 0 && A == Z +
396 std::stringstream ss;
397 ss << "unphysical residual fragment
398 << " skipping it and resampling
399 theInterfaceStore->EmitWarning(ss.st
400 eventIsOK = false;
401 continue;
402 }
403 const G4double kinE = eventInfo.EKinRe
404 const G4double px = eventInfo.pxRem[i]
405 const G4double py = eventInfo.pyRem[i]
406 const G4double pz = eventInfo.pzRem[i]
407 G4ThreeVector spin(
408 eventInfo.jxRem[i]*hbar_Planck,
409 eventInfo.jyRem[i]*hbar_Planck,
410 eventInfo.jzRem[i]*hbar_Planck
411 );
412 const G4double excitationE = eventInfo
413 G4double nuclearMass = excitationE;
414 if ( S == 0 ) {
415 nuclearMass += G4NucleiProperties::G
416 } else {
417 // Assumed that the opposite of the strang
418 nuclearMass += G4HyperNucleiProperti
419 }
420 const G4double scaling = remnant4Momen
421 G4LorentzVector fourMomentum(scaling *
422 nuclearMass + kinE);
423 if(std::abs(scaling - 1.0) > 0.01) {
424 std::stringstream ss;
425 ss << "momentum scaling = " << scali
426 << "\n Lorentz vec
427 << ")\n A = " << A
428 << "\n E* = " << e
429 << "\n remnant i="
430 << "\n Reaction wa
431 << "-MeV " << trackDefinition->Ge
432 << theIonTable->GetIonName(theNuc
433 << ", in " << (inverseKinematics
434 theInterfaceStore->EmitWarning(ss.st
435 }
436
437 // Apply the toLabFrame rotation
438 fourMomentum *= toLabFrame;
439 spin *= toLabFrame3;
440 // Apply the inverse-kinematics boost
441 if(inverseKinematics) {
442 fourMomentum *= *toDirectKinematics;
443 fourMomentum.setVect(-fourMomentum.v
444 }
445
446 fourMomentumOut += fourMomentum;
447 G4Fragment remnant(A, Z, std::abs(S),
448 remnant.SetAngularMomentum(spin);
449 remnant.SetCreatorModelID(secID);
450 if(dumpRemnantInfo) {
451 G4cerr << "G4INCLXX_DUMP_REMNANT: "
452 }
453 remnants.push_back(remnant);
454 }
455
456 // Give up is the event is not ok (e.g.
457 if(!eventIsOK) {
458 const G4int nSecondaries = (G4int)theR
459 for(G4int j=0; j<nSecondaries; ++j) de
460 theResult.Clear();
461 theResult.SetStatusChange(stopAndKill)
462 remnants.clear();
463 } else {
464 // Check four-momentum conservation
465 const G4LorentzVector violation4Moment
466 const G4double energyViolation = std::
467 const G4double momentumViolation = vio
468 if(energyViolation > G4INCLXXInterface
469 std::stringstream ss;
470 ss << "energy conservation violated
471 << aTrack.GetKineticEnergy()/MeV
472 << " + " << theIonTable->GetIonNa
473 << " inelastic reaction, in " <<
474 theInterfaceStore->EmitWarning(ss.st
475 eventIsOK = false;
476 const G4int nSecondaries = (G4int)th
477 for(G4int j=0; j<nSecondaries; ++j)
478 theResult.Clear();
479 theResult.SetStatusChange(stopAndKil
480 remnants.clear();
481 } else if(momentumViolation > G4INCLXX
482 std::stringstream ss;
483 ss << "momentum conservation violate
484 << aTrack.GetKineticEnergy()/MeV
485 << " + " << theIonTable->GetIonNa
486 << " inelastic reaction, in " <<
487 theInterfaceStore->EmitWarning(ss.st
488 eventIsOK = false;
489 const G4int nSecondaries = (G4int)th
490 for(G4int j=0; j<nSecondaries; ++j)
491 theResult.Clear();
492 theResult.SetStatusChange(stopAndKil
493 remnants.clear();
494 }
495 }
496 }
497 nTries++;
498 } while(!eventIsOK && nTries < maxTries); /*
499
500 // Clean up the objects that we created for
501 if(inverseKinematics) {
502 delete aProjectileTrack;
503 delete theTargetNucleus;
504 delete toInverseKinematics;
505 delete toDirectKinematics;
506 }
507
508 if(!eventIsOK) {
509 std::stringstream ss;
510 ss << "maximum number of tries exceeded fo
511 << aTrack.GetKineticEnergy()/MeV << "-MeV
512 << " + " << theIonTable->GetIonName(nucleu
513 << " inelastic reaction, in " << (inverseK
514 theInterfaceStore->EmitWarning(ss.str());
515 theResult.SetStatusChange(isAlive);
516 theResult.SetEnergyChange(aTrack.GetKineti
517 theResult.SetMomentumChange(aTrack.Get4Mom
518 return &theResult;
519 }
520
521 // De-excitation:
522
523 if(theDeExcitation != 0) {
524 for(std::list<G4Fragment>::iterator i = re
525 i != remnants.end(); ++i) {
526 G4ReactionProductVector *deExcitationRes
527
528 for(G4ReactionProductVector::iterator fr
529 fragment != deExcitationResult->end(); ++f
530 const G4ParticleDefinition *def = (*fragment
531 if(def != 0) {
532 G4DynamicParticle *theFragment = new G4Dyn
533 theResult.AddSecondary(theFragment, (*frag
534 }
535 }
536
537 for(G4ReactionProductVector::iterator fr
538 fragment != deExcitationResult->end(); ++f
539 delete (*fragment);
540 }
541 deExcitationResult->clear();
542 delete deExcitationResult;
543 }
544 }
545
546 remnants.clear();
547
548 if((theTally = theInterfaceStore->GetTally()
549 theTally->Tally(aTrack, theNucleus, theRes
550
551 return &theResult;
552 }
553
554 G4ReactionProductVector* G4INCLXXInterface::Pr
555 return 0;
556 }
557
558 G4INCL::ParticleType G4INCLXXInterface::toINCL
559 if( pdef == G4Proton::Proton())
560 else if(pdef == G4Neutron::Neutron())
561 else if(pdef == G4PionPlus::PionPlus())
562 else if(pdef == G4PionMinus::PionMinus())
563 else if(pdef == G4PionZero::PionZero())
564 else if(pdef == G4KaonPlus::KaonPlus())
565 else if(pdef == G4KaonZero::KaonZero())
566 else if(pdef == G4KaonMinus::KaonMinus())
567 else if(pdef == G4AntiKaonZero::AntiKaonZero
568 // For K0L & K0S we do not take into account
569 else if(pdef == G4KaonZeroLong::KaonZeroLong
570 else if(pdef == G4KaonZeroShort::KaonZeroSho
571 else if(pdef == G4Deuteron::Deuteron())
572 else if(pdef == G4Triton::Triton())
573 else if(pdef == G4He3::He3())
574 else if(pdef == G4Alpha::Alpha())
575 else if(pdef == G4AntiProton::AntiProton())
576 else if(pdef->GetParticleType() == G4Generic
577 else
578 }
579
580 G4INCL::ParticleSpecies G4INCLXXInterface::toI
581 const G4ParticleDefinition *pdef = aTrack.Ge
582 const G4INCL::ParticleType theType = toINCLP
583 if(theType!=G4INCL::Composite)
584 return G4INCL::ParticleSpecies(theType);
585 else {
586 G4INCL::ParticleSpecies theSpecies;
587 theSpecies.theType=theType;
588 theSpecies.theA=pdef->GetAtomicMass();
589 theSpecies.theZ=pdef->GetAtomicNumber();
590 return theSpecies;
591 }
592 }
593
594 G4double G4INCLXXInterface::toINCLKineticEnerg
595 return aTrack.GetKineticEnergy();
596 }
597
598 G4ParticleDefinition *G4INCLXXInterface::toG4P
599 if (PDGCode == 2212) { return G4Proton
600 } else if(PDGCode == 2112) { return G4Neutro
601 } else if(PDGCode == 211) { return G4PionPl
602 } else if(PDGCode == 111) { return G4PionZe
603 } else if(PDGCode == -211) { return G4PionMi
604
605 } else if(PDGCode == 221) { return G4Eta::E
606 } else if(PDGCode == 22) { return G4Gamma:
607
608 } else if(PDGCode == 3122) { return G4Lambda
609 } else if(PDGCode == 3222) { return G4SigmaP
610 } else if(PDGCode == 3212) { return G4SigmaZ
611 } else if(PDGCode == 3112) { return G4SigmaM
612 } else if(PDGCode == 321) { return G4KaonPl
613 } else if(PDGCode == -321) { return G4KaonMi
614 } else if(PDGCode == 130) { return G4KaonZe
615 } else if(PDGCode == 310) { return G4KaonZe
616
617 } else if(PDGCode == 1002) { return G4Deuter
618 } else if(PDGCode == 1003) { return G4Triton
619 } else if(PDGCode == 2003) { return G4He3::H
620 } else if(PDGCode == 2004) { return G4Alpha:
621
622 } else if(PDGCode == -2212) { return G4AntiP
623 } else if(S != 0) { // Assumed that -S give
624 if (A == 3 && Z == 1 && S == -1 ) return G
625 if (A == 4 && Z == 1 && S == -1 ) return G
626 if (A == 4 && Z == 2 && S == -1 ) return G
627 if (A == 4 && Z == 1 && S == -2 ) return G
628 if (A == 4 && Z == 0 && S == -2 ) return G
629 if (A == 5 && Z == 2 && S == -1 ) return G
630 } else if(A > 0 && Z > 0 && A > Z) { // Retu
631 return theIonTable->GetIon(Z, A, 0);
632 }
633 return 0; // Error, unrecognized particle
634 }
635
636 G4DynamicParticle *G4INCLXXInterface::toG4Part
637 G4double kinE, G4double px,
638
639 const G4ParticleDefinition *def = toG4Partic
640 if(def == 0) { // Check if we have a valid p
641 return 0;
642 }
643 const G4double energy = kinE * MeV;
644 const G4ThreeVector momentum(px, py, pz);
645 const G4ThreeVector momentumDirection = mome
646 G4DynamicParticle *p = new G4DynamicParticle
647 return p;
648 }
649
650 G4double G4INCLXXInterface::remnant4MomentumSc
651 G4double kineticE,
652 G4double px, G4double py,
653 G4double pz) const {
654 const G4double p2 = px*px + py*py + pz*pz;
655 if(p2 > 0.0) {
656 const G4double pnew2 = kineticE*kineticE +
657 return std::sqrt(pnew2)/std::sqrt(p2);
658 } else {
659 return 1.0;
660 }
661 }
662
663 void G4INCLXXInterface::ModelDescription(std::
664 outFile
665 << "The Liège Intranuclear Cascade (INCL
666 << "by nucleons, pions and light ion on a
667 << "described as an avalanche of binary n
668 << "lead to the emission of energetic par
669 << "excited thermalised nucleus (remnant)
670 << "outside the scope of INCL++ and is ty
671 << "INCL++ has been reasonably well teste
672 << "pion (idem) and light-ion projectiles
673 << "Most tests involved target nuclei clo
674 << "numbers between 4 and 250.\n\n"
675 << "Reference: D. Mancusi et al., Phys. R
676 }
677
678 G4String const &G4INCLXXInterface::GetDeExcita
679 return theDeExcitation->GetModelName();
680 }
681