| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
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 Helsinki Institute of Physics, Finland
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.hh"
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(G4VPreCompoundModel * const aPreCompound) :
71 G4VIntraNuclearTransportModel(G4INCLXXInterfaceStore::GetInstance()->getINCLXXVersionName()),
72 theINCLModel(NULL),
73 thePreCompoundModel(aPreCompound),
74 theInterfaceStore(G4INCLXXInterfaceStore::GetInstance()),
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()->FindModel("PRECO");
86 thePreCompoundModel = static_cast<G4VPreCompoundModel*>(p);
87 if(!thePreCompoundModel) { thePreCompoundModel = new G4PreCompoundModel; }
88 }
89
90 // Use the environment variable G4INCLXX_NO_DE_EXCITATION to disable de-excitation
91 if(std::getenv("G4INCLXX_NO_DE_EXCITATION")) {
92 G4String message = "de-excitation is completely disabled!";
93 theInterfaceStore->EmitWarning(message);
94 theDeExcitation = 0;
95 } else {
96 G4HadronicInteraction* p =
97 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
98 theDeExcitation = static_cast<G4VPreCompoundModel*>(p);
99 if(!theDeExcitation) { theDeExcitation = new G4PreCompoundModel; }
100
101 // set the fission parameters for G4ExcitationHandler
102 G4VEvaporationChannel * const theFissionChannel =
103 theDeExcitation->GetExcitationHandler()->GetEvaporation()->GetFissionChannel();
104 G4CompetitiveFission * const theFissionChannelCast = dynamic_cast<G4CompetitiveFission *>(theFissionChannel);
105 if(theFissionChannelCast) {
106 theINCLXXLevelDensity = new G4FissionLevelDensityParameterINCLXX;
107 theFissionChannelCast->SetLevelDensityParameter(theINCLXXLevelDensity);
108 theINCLXXFissionProbability = new G4FissionProbability;
109 theINCLXXFissionProbability->SetFissionLevelDensityParameter(theINCLXXLevelDensity);
110 theFissionChannelCast->SetEmissionStrategy(theINCLXXFissionProbability);
111 theInterfaceStore->EmitBigWarning("INCL++/G4ExcitationHandler uses its own level-density parameter for fission");
112 } else {
113 theInterfaceStore->EmitBigWarning("INCL++/G4ExcitationHandler could not use its own level-density parameter for fission");
114 }
115 }
116
117 // use the envvar G4INCLXX_DUMP_REMNANT to dump information about the
118 // remnants on stdout
119 if(std::getenv("G4INCLXX_DUMP_REMNANT"))
120 dumpRemnantInfo = true;
121 else
122 dumpRemnantInfo = false;
123
124 theBackupModel = new G4BinaryLightIonReaction;
125 theBackupModelNucleon = new G4BinaryCascade;
126 secID = G4PhysicsModelCatalog::GetModelID( "model_INCLXXCascade" );
127 }
128
129 G4INCLXXInterface::~G4INCLXXInterface()
130 {
131 delete theINCLXXLevelDensity;
132 delete theINCLXXFissionProbability;
133 }
134
135 G4bool G4INCLXXInterface::AccurateProjectile(const G4HadProjectile &aTrack, const G4Nucleus &theNucleus) const {
136 // Use direct kinematics if the projectile is a nucleon or a pion
137 const G4ParticleDefinition *projectileDef = aTrack.GetDefinition();
138 if(std::abs(projectileDef->GetBaryonNumber()) < 2) // Every non-composite particle. abs()-> in case of anti-nucleus projectile
139 return false;
140
141 // Here all projectiles should be nuclei
142 const G4int pA = projectileDef->GetAtomicMass();
143 if(pA<=0) {
144 std::stringstream ss;
145 ss << "the model does not know how to handle a collision between a "
146 << projectileDef->GetParticleName() << " projectile and a Z="
147 << theNucleus.GetZ_asInt() << ", A=" << theNucleus.GetA_asInt();
148 theInterfaceStore->EmitBigWarning(ss.str());
149 return true;
150 }
151
152 // If either nucleus is a LCP (A<=4), run the collision as light on heavy
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 theMaxProjMassINCL, run the collision
162 // as light on heavy.
163 // Note that here we are sure that either the projectile or the target is
164 // smaller than theMaxProjMass; otherwise theBackupModel would have been
165 // called.
166 const G4int theMaxProjMassINCL = theInterfaceStore->GetMaxProjMassINCL();
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 setting
173 return theInterfaceStore->GetAccurateProjectile();
174 }
175
176 G4HadFinalState* G4INCLXXInterface::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
177 {
178 G4ParticleDefinition const * const trackDefinition = aTrack.GetDefinition();
179 const G4bool isIonTrack = trackDefinition->GetParticleType()==G4GenericIon::GenericIon()->GetParticleType();
180 const G4int trackA = trackDefinition->GetAtomicMass();
181 const G4int trackZ = (G4int) trackDefinition->GetPDGCharge();
182 const G4int trackL = trackDefinition->GetNumberOfLambdasInHypernucleus();
183 const G4int nucleusA = theNucleus.GetA_asInt();
184 const G4int nucleusZ = theNucleus.GetZ_asInt();
185
186 // For reactions induced by weird projectiles (e.g. He2), bail out
187 if((isIonTrack && ((trackZ<=0 && trackL==0) || trackA<=trackZ)) ||
188 (nucleusA>1 && (nucleusZ<=0 || nucleusA<=nucleusZ))) {
189 theResult.Clear();
190 theResult.SetStatusChange(isAlive);
191 theResult.SetEnergyChange(aTrack.GetKineticEnergy());
192 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
193 return &theResult;
194 }
195
196 // For reactions on nucleons, use the backup model (without complaining),
197 // except for anti_proton projectile (in this case, INCLXX is used).
198 if(trackA<=1 && nucleusA<=1 && (trackZ>=0 || trackA==0)) {
199 return theBackupModelNucleon->ApplyYourself(aTrack, theNucleus);
200 }
201
202 // For systems heavier than theMaxProjMassINCL, use another model (typically
203 // BIC)
204 const G4int theMaxProjMassINCL = theInterfaceStore->GetMaxProjMassINCL();
205 if(trackA > theMaxProjMassINCL && nucleusA > theMaxProjMassINCL) {
206 if(!complainedAboutBackupModel) {
207 complainedAboutBackupModel = true;
208 std::stringstream ss;
209 ss << "INCL++ refuses to handle reactions between nuclei with A>"
210 << theMaxProjMassINCL
211 << ". A backup model ("
212 << theBackupModel->GetModelName()
213 << ") will be used instead.";
214 theInterfaceStore->EmitBigWarning(ss.str());
215 }
216 return theBackupModel->ApplyYourself(aTrack, theNucleus);
217 }
218
219 // For energies lower than cascadeMinEnergyPerNucleon, use PreCompound
220 const G4double cascadeMinEnergyPerNucleon = theInterfaceStore->GetCascadeMinEnergyPerNucleon();
221 const G4double trackKinE = aTrack.GetKineticEnergy();
222 if((trackDefinition==G4Neutron::NeutronDefinition() || trackDefinition==G4Proton::ProtonDefinition())
223 && trackKinE < cascadeMinEnergyPerNucleon) {
224 if(!complainedAboutPreCompound) {
225 complainedAboutPreCompound = true;
226 std::stringstream ss;
227 ss << "INCL++ refuses to handle nucleon-induced reactions below "
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(aTrack, theNucleus);
235 }
236
237 // Calculate the total four-momentum in the entrance channel
238 const G4double theNucleusMass = theIonTable->GetIonMass(nucleusZ, nucleusA);
239 const G4double theTrackMass = trackDefinition->GetPDGMass();
240 const G4double theTrackEnergy = trackKinE + theTrackMass;
241 const G4double theTrackMomentumAbs2 = theTrackEnergy*theTrackEnergy - theTrackMass*theTrackMass;
242 const G4double theTrackMomentumAbs = ((theTrackMomentumAbs2>0.0) ? std::sqrt(theTrackMomentumAbs2) : 0.0);
243 const G4ThreeVector theTrackMomentum = aTrack.Get4Momentum().getV().unit() * theTrackMomentumAbs;
244
245 G4LorentzVector goodTrack4Momentum(theTrackMomentum, theTrackEnergy);
246 G4LorentzVector fourMomentumIn;
247 fourMomentumIn.setE(theTrackEnergy + theNucleusMass);
248 fourMomentumIn.setVect(theTrackMomentum);
249
250 // Check if inverse kinematics should be used
251 const G4bool inverseKinematics = AccurateProjectile(aTrack, theNucleus);
252
253 // If we are running in inverse kinematics, redefine aTrack and theNucleus
254 G4LorentzRotation *toInverseKinematics = NULL;
255 G4LorentzRotation *toDirectKinematics = NULL;
256 G4HadProjectile const *aProjectileTrack = &aTrack;
257 G4Nucleus *theTargetNucleus = &theNucleus;
258 if(inverseKinematics) {
259 G4ParticleDefinition *oldTargetDef = theIonTable->GetIon(nucleusZ, nucleusA, 0);
260 const G4ParticleDefinition *oldProjectileDef = trackDefinition;
261
262 if(oldProjectileDef != 0 && oldTargetDef != 0) {
263 const G4int newTargetA = oldProjectileDef->GetAtomicMass();
264 const G4int newTargetZ = oldProjectileDef->GetAtomicNumber();
265 const G4int newTargetL = oldProjectileDef->GetNumberOfLambdasInHypernucleus();
266 if(newTargetA > 0 && newTargetZ > 0) {
267 // This should give us the same energy per nucleon
268 theTargetNucleus = new G4Nucleus(newTargetA, newTargetZ, newTargetL);
269 toInverseKinematics = new G4LorentzRotation(goodTrack4Momentum.boostVector());
270 G4LorentzVector theProjectile4Momentum(0.0, 0.0, 0.0, theNucleusMass);
271 G4DynamicParticle swappedProjectileParticle(oldTargetDef, (*toInverseKinematics) * theProjectile4Momentum);
272 aProjectileTrack = new G4HadProjectile(swappedProjectileParticle);
273 } else {
274 G4String message = "badly defined target after swapping. Falling back to normal (non-swapped) mode.";
275 theInterfaceStore->EmitWarning(message);
276 toInverseKinematics = new G4LorentzRotation;
277 }
278 } else {
279 G4String message = "oldProjectileDef or oldTargetDef was null";
280 theInterfaceStore->EmitWarning(message);
281 toInverseKinematics = new G4LorentzRotation;
282 }
283 }
284
285 // INCL assumes the projectile particle is going in the direction of
286 // the Z-axis. Here we construct proper rotation to convert the
287 // momentum vectors of the outcoming particles to the original
288 // coordinate system.
289 G4LorentzVector projectileMomentum = aProjectileTrack->Get4Momentum();
290
291 // INCL++ assumes that the projectile is going in the direction of
292 // the z-axis. In principle, if the coordinate system used by G4
293 // hadronic framework is defined differently we need a rotation to
294 // transform the INCL++ reaction products to the appropriate
295 // frame. Please note that it isn't necessary to apply this
296 // transform to the projectile because when creating the INCL++
297 // projectile particle; \see{toINCLKineticEnergy} needs to use only the
298 // projectile energy (direction is simply assumed to be along z-axis).
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 given to us by G4
305 // hadronic framework is already going in the direction of the
306 // z-axis so this rotation is actually unnecessary. Both toZ and
307 // toLabFrame turn out to be unit matrices as can be seen by
308 // uncommenting the folowing two lines:
309 // G4cout <<"toZ = " << toZ << G4endl;
310 // G4cout <<"toLabFrame = " << toLabFrame << G4endl;
311
312 theResult.Clear(); // Make sure the output data structure is clean.
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. However, this is meaningful
320 // only in the standalone code. In Geant4 we must "force" INCL to
321 // produce a valid cascade.
322 G4bool eventIsOK = false;
323 do {
324 const G4INCL::ParticleSpecies theSpecies = toINCLParticleSpecies(*aProjectileTrack);
325 const G4double kineticEnergy = toINCLKineticEnergy(*aProjectileTrack);
326
327 // The INCL model will be created at the first use
328 theINCLModel = G4INCLXXInterfaceStore::GetInstance()->GetINCLModel();
329
330 const G4INCL::EventInfo eventInfo = theINCLModel->processEvent(theSpecies, kineticEnergy,
331 theTargetNucleus->GetA_asInt(),
332 theTargetNucleus->GetZ_asInt(),
333 -theTargetNucleus->GetL()); // Strangeness has opposite sign
334 // eventIsOK = !eventInfo.transparent && nTries < maxTries; // of the number of Lambdas
335 eventIsOK = !eventInfo.transparent;
336 if(eventIsOK) {
337
338 // If the collision was run in inverse kinematics, prepare to boost back
339 // all the reaction products
340 if(inverseKinematics) {
341 toDirectKinematics = new G4LorentzRotation(toInverseKinematics->inverse());
342 }
343
344 G4LorentzVector fourMomentumOut;
345
346 for(G4int i = 0; i < eventInfo.nParticles; ++i) {
347 G4int A = eventInfo.A[i];
348 G4int Z = eventInfo.Z[i];
349 G4int S = eventInfo.S[i]; // Strangeness
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, Z, S, PDGCode, kinE, px, py, pz);
356 if(p != 0) {
357 G4LorentzVector momentum = p->Get4Momentum();
358
359 // Apply the toLabFrame rotation
360 momentum *= toLabFrame;
361 // Apply the inverse-kinematics boost
362 if(inverseKinematics) {
363 momentum *= *toDirectKinematics;
364 momentum.setVect(-momentum.vect());
365 }
366
367 // Set the four-momentum of the reaction products
368 p->Set4Momentum(momentum);
369 fourMomentumOut += momentum;
370
371 // Propagate the particle's parent resonance information
372 G4HadSecondary secondary(p, 1.0, secID);
373 G4ParticleDefinition* parentResonanceDef = nullptr;
374 if ( eventInfo.parentResonancePDGCode[i] != 0 ) {
375 parentResonanceDef = G4ParticleTable::GetParticleTable()->FindParticle(eventInfo.parentResonancePDGCode[i]);
376 }
377 secondary.SetParentResonanceDef(parentResonanceDef);
378 secondary.SetParentResonanceID(eventInfo.parentResonanceID[i]);
379
380 theResult.AddSecondary(secondary);
381
382 } else {
383 G4String message = "the model produced a particle that couldn't be converted to Geant4 particle.";
384 theInterfaceStore->EmitWarning(message);
385 }
386 }
387
388 for(G4int i = 0; i < eventInfo.nRemnants; ++i) {
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 physical bound state: if not, resample the collision.
393 if(( Z == 0 && S == 0 && A > 1 ) || // No bound states for nn, nnn, nnnn, ...
394 ( Z == 0 && S != 0 && A < 4 ) || // No bound states for nl, ll, nnl, nll, lll
395 ( Z != 0 && S != 0 && A == Z + std::abs(S) )) { // No bound states for pl, ppl, pll, ...
396 std::stringstream ss;
397 ss << "unphysical residual fragment : Z=" << Z << " S=" << S << " A=" << A
398 << " skipping it and resampling the collision";
399 theInterfaceStore->EmitWarning(ss.str());
400 eventIsOK = false;
401 continue;
402 }
403 const G4double kinE = eventInfo.EKinRem[i];
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.EStarRem[i];
413 G4double nuclearMass = excitationE;
414 if ( S == 0 ) {
415 nuclearMass += G4NucleiProperties::GetNuclearMass(A, Z);
416 } else {
417 // Assumed that the opposite of the strangeness of the remnant gives the number of Lambdas inside it
418 nuclearMass += G4HyperNucleiProperties::GetNuclearMass(A, Z, std::abs(S));
419 }
420 const G4double scaling = remnant4MomentumScaling(nuclearMass, kinE, px, py, pz);
421 G4LorentzVector fourMomentum(scaling * px, scaling * py, scaling * pz,
422 nuclearMass + kinE);
423 if(std::abs(scaling - 1.0) > 0.01) {
424 std::stringstream ss;
425 ss << "momentum scaling = " << scaling
426 << "\n Lorentz vector = " << fourMomentum
427 << ")\n A = " << A << ", Z = " << Z << ", S = " << S
428 << "\n E* = " << excitationE << ", nuclearMass = " << nuclearMass
429 << "\n remnant i=" << i << ", nRemnants=" << eventInfo.nRemnants
430 << "\n Reaction was: " << aTrack.GetKineticEnergy()/MeV
431 << "-MeV " << trackDefinition->GetParticleName() << " + "
432 << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
433 << ", in " << (inverseKinematics ? "inverse" : "direct") << " kinematics.";
434 theInterfaceStore->EmitWarning(ss.str());
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.vect());
444 }
445
446 fourMomentumOut += fourMomentum;
447 G4Fragment remnant(A, Z, std::abs(S), fourMomentum); // Assumed that -strangeness gives the number of Lambdas
448 remnant.SetAngularMomentum(spin);
449 remnant.SetCreatorModelID(secID);
450 if(dumpRemnantInfo) {
451 G4cerr << "G4INCLXX_DUMP_REMNANT: " << remnant << " spin: " << spin << G4endl;
452 }
453 remnants.push_back(remnant);
454 }
455
456 // Give up is the event is not ok (e.g. unphysical residual)
457 if(!eventIsOK) {
458 const G4int nSecondaries = (G4int)theResult.GetNumberOfSecondaries();
459 for(G4int j=0; j<nSecondaries; ++j) delete theResult.GetSecondary(j)->GetParticle();
460 theResult.Clear();
461 theResult.SetStatusChange(stopAndKill);
462 remnants.clear();
463 } else {
464 // Check four-momentum conservation
465 const G4LorentzVector violation4Momentum = fourMomentumOut - fourMomentumIn;
466 const G4double energyViolation = std::abs(violation4Momentum.e());
467 const G4double momentumViolation = violation4Momentum.rho();
468 if(energyViolation > G4INCLXXInterfaceStore::GetInstance()->GetConservationTolerance()) {
469 std::stringstream ss;
470 ss << "energy conservation violated by " << energyViolation/MeV << " MeV in "
471 << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
472 << " + " << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
473 << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics. Will resample.";
474 theInterfaceStore->EmitWarning(ss.str());
475 eventIsOK = false;
476 const G4int nSecondaries = (G4int)theResult.GetNumberOfSecondaries();
477 for(G4int j=0; j<nSecondaries; ++j) delete theResult.GetSecondary(j)->GetParticle();
478 theResult.Clear();
479 theResult.SetStatusChange(stopAndKill);
480 remnants.clear();
481 } else if(momentumViolation > G4INCLXXInterfaceStore::GetInstance()->GetConservationTolerance()) {
482 std::stringstream ss;
483 ss << "momentum conservation violated by " << momentumViolation/MeV << " MeV in "
484 << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
485 << " + " << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
486 << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics. Will resample.";
487 theInterfaceStore->EmitWarning(ss.str());
488 eventIsOK = false;
489 const G4int nSecondaries = (G4int)theResult.GetNumberOfSecondaries();
490 for(G4int j=0; j<nSecondaries; ++j) delete theResult.GetSecondary(j)->GetParticle();
491 theResult.Clear();
492 theResult.SetStatusChange(stopAndKill);
493 remnants.clear();
494 }
495 }
496 }
497 nTries++;
498 } while(!eventIsOK && nTries < maxTries); /* Loop checking, 10.07.2015, D.Mancusi */
499
500 // Clean up the objects that we created for the inverse kinematics
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 for the proposed "
511 << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
512 << " + " << theIonTable->GetIonName(nucleusZ, nucleusA, 0)
513 << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics.";
514 theInterfaceStore->EmitWarning(ss.str());
515 theResult.SetStatusChange(isAlive);
516 theResult.SetEnergyChange(aTrack.GetKineticEnergy());
517 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
518 return &theResult;
519 }
520
521 // De-excitation:
522
523 if(theDeExcitation != 0) {
524 for(std::list<G4Fragment>::iterator i = remnants.begin();
525 i != remnants.end(); ++i) {
526 G4ReactionProductVector *deExcitationResult = theDeExcitation->DeExcite((*i));
527
528 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
529 fragment != deExcitationResult->end(); ++fragment) {
530 const G4ParticleDefinition *def = (*fragment)->GetDefinition();
531 if(def != 0) {
532 G4DynamicParticle *theFragment = new G4DynamicParticle(def, (*fragment)->GetMomentum());
533 theResult.AddSecondary(theFragment, (*fragment)->GetCreatorModelID());
534 }
535 }
536
537 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
538 fragment != deExcitationResult->end(); ++fragment) {
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, theResult);
550
551 return &theResult;
552 }
553
554 G4ReactionProductVector* G4INCLXXInterface::Propagate(G4KineticTrackVector* , G4V3DNucleus* ) {
555 return 0;
556 }
557
558 G4INCL::ParticleType G4INCLXXInterface::toINCLParticleType(G4ParticleDefinition const * const pdef) const {
559 if( pdef == G4Proton::Proton()) return G4INCL::Proton;
560 else if(pdef == G4Neutron::Neutron()) return G4INCL::Neutron;
561 else if(pdef == G4PionPlus::PionPlus()) return G4INCL::PiPlus;
562 else if(pdef == G4PionMinus::PionMinus()) return G4INCL::PiMinus;
563 else if(pdef == G4PionZero::PionZero()) return G4INCL::PiZero;
564 else if(pdef == G4KaonPlus::KaonPlus()) return G4INCL::KPlus;
565 else if(pdef == G4KaonZero::KaonZero()) return G4INCL::KZero;
566 else if(pdef == G4KaonMinus::KaonMinus()) return G4INCL::KMinus;
567 else if(pdef == G4AntiKaonZero::AntiKaonZero()) return G4INCL::KZeroBar;
568 // For K0L & K0S we do not take into account K0/K0B oscillations
569 else if(pdef == G4KaonZeroLong::KaonZeroLong()) return G4UniformRand() < 0.5 ? G4INCL::KZeroBar : G4INCL::KZero;
570 else if(pdef == G4KaonZeroShort::KaonZeroShort()) return G4UniformRand() < 0.5 ? G4INCL::KZeroBar : G4INCL::KZero;
571 else if(pdef == G4Deuteron::Deuteron()) return G4INCL::Composite;
572 else if(pdef == G4Triton::Triton()) return G4INCL::Composite;
573 else if(pdef == G4He3::He3()) return G4INCL::Composite;
574 else if(pdef == G4Alpha::Alpha()) return G4INCL::Composite;
575 else if(pdef == G4AntiProton::AntiProton()) return G4INCL::antiProton;
576 else if(pdef->GetParticleType() == G4GenericIon::GenericIon()->GetParticleType()) return G4INCL::Composite;
577 else return G4INCL::UnknownParticle;
578 }
579
580 G4INCL::ParticleSpecies G4INCLXXInterface::toINCLParticleSpecies(G4HadProjectile const &aTrack) const {
581 const G4ParticleDefinition *pdef = aTrack.GetDefinition();
582 const G4INCL::ParticleType theType = toINCLParticleType(pdef);
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::toINCLKineticEnergy(G4HadProjectile const &aTrack) const {
595 return aTrack.GetKineticEnergy();
596 }
597
598 G4ParticleDefinition *G4INCLXXInterface::toG4ParticleDefinition(G4int A, G4int Z, G4int S, G4int PDGCode) const {
599 if (PDGCode == 2212) { return G4Proton::Proton();
600 } else if(PDGCode == 2112) { return G4Neutron::Neutron();
601 } else if(PDGCode == 211) { return G4PionPlus::PionPlus();
602 } else if(PDGCode == 111) { return G4PionZero::PionZero();
603 } else if(PDGCode == -211) { return G4PionMinus::PionMinus();
604
605 } else if(PDGCode == 221) { return G4Eta::Eta();
606 } else if(PDGCode == 22) { return G4Gamma::Gamma();
607
608 } else if(PDGCode == 3122) { return G4Lambda::Lambda();
609 } else if(PDGCode == 3222) { return G4SigmaPlus::SigmaPlus();
610 } else if(PDGCode == 3212) { return G4SigmaZero::SigmaZero();
611 } else if(PDGCode == 3112) { return G4SigmaMinus::SigmaMinus();
612 } else if(PDGCode == 321) { return G4KaonPlus::KaonPlus();
613 } else if(PDGCode == -321) { return G4KaonMinus::KaonMinus();
614 } else if(PDGCode == 130) { return G4KaonZeroLong::KaonZeroLong();
615 } else if(PDGCode == 310) { return G4KaonZeroShort::KaonZeroShort();
616
617 } else if(PDGCode == 1002) { return G4Deuteron::Deuteron();
618 } else if(PDGCode == 1003) { return G4Triton::Triton();
619 } else if(PDGCode == 2003) { return G4He3::He3();
620 } else if(PDGCode == 2004) { return G4Alpha::Alpha();
621
622 } else if(PDGCode == -2212) { return G4AntiProton::AntiProton();
623 } else if(S != 0) { // Assumed that -S gives the number of Lambdas
624 if (A == 3 && Z == 1 && S == -1 ) return G4HyperTriton::Definition();
625 if (A == 4 && Z == 1 && S == -1 ) return G4HyperH4::Definition();
626 if (A == 4 && Z == 2 && S == -1 ) return G4HyperAlpha::Definition();
627 if (A == 4 && Z == 1 && S == -2 ) return G4DoubleHyperH4::Definition();
628 if (A == 4 && Z == 0 && S == -2 ) return G4DoubleHyperDoubleNeutron::Definition();
629 if (A == 5 && Z == 2 && S == -1 ) return G4HyperHe5::Definition();
630 } else if(A > 0 && Z > 0 && A > Z) { // Returns ground state ion definition.
631 return theIonTable->GetIon(Z, A, 0);
632 }
633 return 0; // Error, unrecognized particle
634 }
635
636 G4DynamicParticle *G4INCLXXInterface::toG4Particle(G4int A, G4int Z, G4int S, G4int PDGCode,
637 G4double kinE, G4double px,
638 G4double py, G4double pz) const {
639 const G4ParticleDefinition *def = toG4ParticleDefinition(A, Z, S, PDGCode);
640 if(def == 0) { // Check if we have a valid particle definition
641 return 0;
642 }
643 const G4double energy = kinE * MeV;
644 const G4ThreeVector momentum(px, py, pz);
645 const G4ThreeVector momentumDirection = momentum.unit();
646 G4DynamicParticle *p = new G4DynamicParticle(def, momentumDirection, energy);
647 return p;
648 }
649
650 G4double G4INCLXXInterface::remnant4MomentumScaling(G4double mass,
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 + 2.0*kineticE*mass;
657 return std::sqrt(pnew2)/std::sqrt(p2);
658 } else {
659 return 1.0;
660 }
661 }
662
663 void G4INCLXXInterface::ModelDescription(std::ostream& outFile) const {
664 outFile
665 << "The Liège Intranuclear Cascade (INCL++) is a model for reactions induced\n"
666 << "by nucleons, pions and light ion on any nucleus. The reaction is\n"
667 << "described as an avalanche of binary nucleon-nucleon collisions, which can\n"
668 << "lead to the emission of energetic particles and to the formation of an\n"
669 << "excited thermalised nucleus (remnant). The de-excitation of the remnant is\n"
670 << "outside the scope of INCL++ and is typically described by another model.\n\n"
671 << "INCL++ has been reasonably well tested for nucleon (~50 MeV to ~15 GeV),\n"
672 << "pion (idem) and light-ion projectiles (up to A=18, ~10A MeV to 1A GeV).\n"
673 << "Most tests involved target nuclei close to the stability valley, with\n"
674 << "numbers between 4 and 250.\n\n"
675 << "Reference: D. Mancusi et al., Phys. Rev. C90 (2014) 054602\n\n";
676 }
677
678 G4String const &G4INCLXXInterface::GetDeExcitationModelName() const {
679 return theDeExcitation->GetModelName();
680 }
681