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