| 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(std::abs(projectileDef->GetBaryonNumber()) < 2) // Every non-composite particle. abs()-> in case of anti-nucleus projectile
139 return false; 126 return false;
140 127
141 // Here all projectiles should be nuclei 128 // Here all projectiles should be nuclei
142 const G4int pA = projectileDef->GetAtomicMas 129 const G4int pA = projectileDef->GetAtomicMass();
143 if(pA<=0) { 130 if(pA<=0) {
144 std::stringstream ss; 131 std::stringstream ss;
145 ss << "the model does not know how to hand 132 ss << "the model does not know how to handle a collision between a "
146 << projectileDef->GetParticleName() << " 133 << projectileDef->GetParticleName() << " projectile and a Z="
147 << theNucleus.GetZ_asInt() << ", A=" << 134 << theNucleus.GetZ_asInt() << ", A=" << theNucleus.GetA_asInt();
148 theInterfaceStore->EmitBigWarning(ss.str() 135 theInterfaceStore->EmitBigWarning(ss.str());
149 return true; 136 return true;
150 } 137 }
151 138
152 // If either nucleus is a LCP (A<=4), run th 139 // If either nucleus is a LCP (A<=4), run the collision as light on heavy
153 const G4int tA = theNucleus.GetA_asInt(); 140 const G4int tA = theNucleus.GetA_asInt();
154 if(tA<=4 || pA<=4) { 141 if(tA<=4 || pA<=4) {
155 if(pA<tA) 142 if(pA<tA)
156 return false; 143 return false;
157 else 144 else
158 return true; 145 return true;
159 } 146 }
160 147
161 // If one of the nuclei is heavier than theM 148 // If one of the nuclei is heavier than theMaxProjMassINCL, run the collision
162 // as light on heavy. 149 // as light on heavy.
163 // Note that here we are sure that either th 150 // Note that here we are sure that either the projectile or the target is
164 // smaller than theMaxProjMass; otherwise th 151 // smaller than theMaxProjMass; otherwise theBackupModel would have been
165 // called. 152 // called.
166 const G4int theMaxProjMassINCL = theInterfac 153 const G4int theMaxProjMassINCL = theInterfaceStore->GetMaxProjMassINCL();
167 if(pA > theMaxProjMassINCL) 154 if(pA > theMaxProjMassINCL)
168 return true; 155 return true;
169 else if(tA > theMaxProjMassINCL) 156 else if(tA > theMaxProjMassINCL)
170 return false; 157 return false;
171 else 158 else
172 // In all other cases, use the global sett 159 // In all other cases, use the global setting
173 return theInterfaceStore->GetAccurateProje 160 return theInterfaceStore->GetAccurateProjectile();
174 } 161 }
175 162
176 G4HadFinalState* G4INCLXXInterface::ApplyYours 163 G4HadFinalState* G4INCLXXInterface::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
177 { 164 {
178 G4ParticleDefinition const * const trackDefi 165 G4ParticleDefinition const * const trackDefinition = aTrack.GetDefinition();
179 const G4bool isIonTrack = trackDefinition->G 166 const G4bool isIonTrack = trackDefinition->GetParticleType()==G4GenericIon::GenericIon()->GetParticleType();
180 const G4int trackA = trackDefinition->GetAto 167 const G4int trackA = trackDefinition->GetAtomicMass();
181 const G4int trackZ = (G4int) trackDefinition 168 const G4int trackZ = (G4int) trackDefinition->GetPDGCharge();
182 const G4int trackL = trackDefinition->GetNum <<
183 const G4int nucleusA = theNucleus.GetA_asInt 169 const G4int nucleusA = theNucleus.GetA_asInt();
184 const G4int nucleusZ = theNucleus.GetZ_asInt 170 const G4int nucleusZ = theNucleus.GetZ_asInt();
185 171
186 // For reactions induced by weird projectile 172 // For reactions induced by weird projectiles (e.g. He2), bail out
187 if((isIonTrack && ((trackZ<=0 && trackL==0) << 173 if((isIonTrack && (trackZ<=0 || trackA<=trackZ)) || (nucleusA>1 && (nucleusZ<=0 || nucleusA<=nucleusZ))) {
188 (nucleusA>1 && (nucleusZ<= <<
189 theResult.Clear(); 174 theResult.Clear();
190 theResult.SetStatusChange(isAlive); 175 theResult.SetStatusChange(isAlive);
191 theResult.SetEnergyChange(aTrack.GetKineti 176 theResult.SetEnergyChange(aTrack.GetKineticEnergy());
192 theResult.SetMomentumChange(aTrack.Get4Mom 177 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
193 return &theResult; 178 return &theResult;
194 } 179 }
195 180
196 // For reactions on nucleons, use the backup << 181 // For reactions on nucleons, use the backup model (without complaining)
197 // except for anti_proton projectile (in thi << 182 if(trackA<=1 && nucleusA<=1) {
198 if(trackA<=1 && nucleusA<=1 && (trackZ>=0 || <<
199 return theBackupModelNucleon->ApplyYoursel 183 return theBackupModelNucleon->ApplyYourself(aTrack, theNucleus);
200 } 184 }
201 << 185
202 // For systems heavier than theMaxProjMassIN 186 // For systems heavier than theMaxProjMassINCL, use another model (typically
203 // BIC) 187 // BIC)
204 const G4int theMaxProjMassINCL = theInterfac 188 const G4int theMaxProjMassINCL = theInterfaceStore->GetMaxProjMassINCL();
205 if(trackA > theMaxProjMassINCL && nucleusA > 189 if(trackA > theMaxProjMassINCL && nucleusA > theMaxProjMassINCL) {
206 if(!complainedAboutBackupModel) { 190 if(!complainedAboutBackupModel) {
207 complainedAboutBackupModel = true; 191 complainedAboutBackupModel = true;
208 std::stringstream ss; 192 std::stringstream ss;
209 ss << "INCL++ refuses to handle reaction 193 ss << "INCL++ refuses to handle reactions between nuclei with A>"
210 << theMaxProjMassINCL 194 << theMaxProjMassINCL
211 << ". A backup model (" 195 << ". A backup model ("
212 << theBackupModel->GetModelName() 196 << theBackupModel->GetModelName()
213 << ") will be used instead."; 197 << ") will be used instead.";
214 theInterfaceStore->EmitBigWarning(ss.str 198 theInterfaceStore->EmitBigWarning(ss.str());
215 } 199 }
216 return theBackupModel->ApplyYourself(aTrac 200 return theBackupModel->ApplyYourself(aTrack, theNucleus);
217 } 201 }
218 202
219 // For energies lower than cascadeMinEnergyP 203 // For energies lower than cascadeMinEnergyPerNucleon, use PreCompound
220 const G4double cascadeMinEnergyPerNucleon = 204 const G4double cascadeMinEnergyPerNucleon = theInterfaceStore->GetCascadeMinEnergyPerNucleon();
221 const G4double trackKinE = aTrack.GetKinetic 205 const G4double trackKinE = aTrack.GetKineticEnergy();
222 if((trackDefinition==G4Neutron::NeutronDefin 206 if((trackDefinition==G4Neutron::NeutronDefinition() || trackDefinition==G4Proton::ProtonDefinition())
223 && trackKinE < cascadeMinEnergyPerNucleo 207 && trackKinE < cascadeMinEnergyPerNucleon) {
224 if(!complainedAboutPreCompound) { 208 if(!complainedAboutPreCompound) {
225 complainedAboutPreCompound = true; 209 complainedAboutPreCompound = true;
226 std::stringstream ss; 210 std::stringstream ss;
227 ss << "INCL++ refuses to handle nucleon- 211 ss << "INCL++ refuses to handle nucleon-induced reactions below "
228 << cascadeMinEnergyPerNucleon / MeV 212 << cascadeMinEnergyPerNucleon / MeV
229 << " MeV. A PreCoumpound model (" 213 << " MeV. A PreCoumpound model ("
230 << thePreCompoundModel->GetModelName() 214 << thePreCompoundModel->GetModelName()
231 << ") will be used instead."; 215 << ") will be used instead.";
232 theInterfaceStore->EmitBigWarning(ss.str 216 theInterfaceStore->EmitBigWarning(ss.str());
233 } 217 }
234 return thePreCompoundModel->ApplyYourself( 218 return thePreCompoundModel->ApplyYourself(aTrack, theNucleus);
235 } 219 }
236 220
237 // Calculate the total four-momentum in the 221 // Calculate the total four-momentum in the entrance channel
238 const G4double theNucleusMass = theIonTable- 222 const G4double theNucleusMass = theIonTable->GetIonMass(nucleusZ, nucleusA);
239 const G4double theTrackMass = trackDefinitio 223 const G4double theTrackMass = trackDefinition->GetPDGMass();
240 const G4double theTrackEnergy = trackKinE + 224 const G4double theTrackEnergy = trackKinE + theTrackMass;
241 const G4double theTrackMomentumAbs2 = theTra 225 const G4double theTrackMomentumAbs2 = theTrackEnergy*theTrackEnergy - theTrackMass*theTrackMass;
242 const G4double theTrackMomentumAbs = ((theTr 226 const G4double theTrackMomentumAbs = ((theTrackMomentumAbs2>0.0) ? std::sqrt(theTrackMomentumAbs2) : 0.0);
243 const G4ThreeVector theTrackMomentum = aTrac 227 const G4ThreeVector theTrackMomentum = aTrack.Get4Momentum().getV().unit() * theTrackMomentumAbs;
244 228
245 G4LorentzVector goodTrack4Momentum(theTrackM 229 G4LorentzVector goodTrack4Momentum(theTrackMomentum, theTrackEnergy);
246 G4LorentzVector fourMomentumIn; 230 G4LorentzVector fourMomentumIn;
247 fourMomentumIn.setE(theTrackEnergy + theNucl 231 fourMomentumIn.setE(theTrackEnergy + theNucleusMass);
248 fourMomentumIn.setVect(theTrackMomentum); 232 fourMomentumIn.setVect(theTrackMomentum);
249 233
250 // Check if inverse kinematics should be use 234 // Check if inverse kinematics should be used
251 const G4bool inverseKinematics = AccuratePro 235 const G4bool inverseKinematics = AccurateProjectile(aTrack, theNucleus);
252 236
253 // If we are running in inverse kinematics, 237 // If we are running in inverse kinematics, redefine aTrack and theNucleus
254 G4LorentzRotation *toInverseKinematics = NUL 238 G4LorentzRotation *toInverseKinematics = NULL;
255 G4LorentzRotation *toDirectKinematics = NULL 239 G4LorentzRotation *toDirectKinematics = NULL;
256 G4HadProjectile const *aProjectileTrack = &a 240 G4HadProjectile const *aProjectileTrack = &aTrack;
257 G4Nucleus *theTargetNucleus = &theNucleus; 241 G4Nucleus *theTargetNucleus = &theNucleus;
258 if(inverseKinematics) { 242 if(inverseKinematics) {
259 G4ParticleDefinition *oldTargetDef = theIo 243 G4ParticleDefinition *oldTargetDef = theIonTable->GetIon(nucleusZ, nucleusA, 0);
260 const G4ParticleDefinition *oldProjectileD 244 const G4ParticleDefinition *oldProjectileDef = trackDefinition;
261 245
262 if(oldProjectileDef != 0 && oldTargetDef ! 246 if(oldProjectileDef != 0 && oldTargetDef != 0) {
263 const G4int newTargetA = oldProjectileDe 247 const G4int newTargetA = oldProjectileDef->GetAtomicMass();
264 const G4int newTargetZ = oldProjectileDe 248 const G4int newTargetZ = oldProjectileDef->GetAtomicNumber();
265 const G4int newTargetL = oldProjectileDe << 249
266 if(newTargetA > 0 && newTargetZ > 0) { 250 if(newTargetA > 0 && newTargetZ > 0) {
267 // This should give us the same energy 251 // This should give us the same energy per nucleon
268 theTargetNucleus = new G4Nucleus(newTa << 252 theTargetNucleus = new G4Nucleus(newTargetA, newTargetZ);
269 toInverseKinematics = new G4LorentzRot 253 toInverseKinematics = new G4LorentzRotation(goodTrack4Momentum.boostVector());
270 G4LorentzVector theProjectile4Momentum 254 G4LorentzVector theProjectile4Momentum(0.0, 0.0, 0.0, theNucleusMass);
271 G4DynamicParticle swappedProjectilePar 255 G4DynamicParticle swappedProjectileParticle(oldTargetDef, (*toInverseKinematics) * theProjectile4Momentum);
272 aProjectileTrack = new G4HadProjectile 256 aProjectileTrack = new G4HadProjectile(swappedProjectileParticle);
273 } else { 257 } else {
274 G4String message = "badly defined targ 258 G4String message = "badly defined target after swapping. Falling back to normal (non-swapped) mode.";
275 theInterfaceStore->EmitWarning(message 259 theInterfaceStore->EmitWarning(message);
276 toInverseKinematics = new G4LorentzRot 260 toInverseKinematics = new G4LorentzRotation;
277 } 261 }
278 } else { 262 } else {
279 G4String message = "oldProjectileDef or 263 G4String message = "oldProjectileDef or oldTargetDef was null";
280 theInterfaceStore->EmitWarning(message); 264 theInterfaceStore->EmitWarning(message);
281 toInverseKinematics = new G4LorentzRotat 265 toInverseKinematics = new G4LorentzRotation;
282 } 266 }
283 } 267 }
284 268
285 // INCL assumes the projectile particle is g 269 // INCL assumes the projectile particle is going in the direction of
286 // the Z-axis. Here we construct proper rota 270 // the Z-axis. Here we construct proper rotation to convert the
287 // momentum vectors of the outcoming particl 271 // momentum vectors of the outcoming particles to the original
288 // coordinate system. 272 // coordinate system.
289 G4LorentzVector projectileMomentum = aProjec 273 G4LorentzVector projectileMomentum = aProjectileTrack->Get4Momentum();
290 274
291 // INCL++ assumes that the projectile is goi 275 // INCL++ assumes that the projectile is going in the direction of
292 // the z-axis. In principle, if the coordina 276 // the z-axis. In principle, if the coordinate system used by G4
293 // hadronic framework is defined differently 277 // hadronic framework is defined differently we need a rotation to
294 // transform the INCL++ reaction products to 278 // transform the INCL++ reaction products to the appropriate
295 // frame. Please note that it isn't necessar 279 // frame. Please note that it isn't necessary to apply this
296 // transform to the projectile because when 280 // transform to the projectile because when creating the INCL++
297 // projectile particle; \see{toINCLKineticEn 281 // projectile particle; \see{toINCLKineticEnergy} needs to use only the
298 // projectile energy (direction is simply as 282 // projectile energy (direction is simply assumed to be along z-axis).
299 G4RotationMatrix toZ; 283 G4RotationMatrix toZ;
300 toZ.rotateZ(-projectileMomentum.phi()); 284 toZ.rotateZ(-projectileMomentum.phi());
301 toZ.rotateY(-projectileMomentum.theta()); 285 toZ.rotateY(-projectileMomentum.theta());
302 G4RotationMatrix toLabFrame3 = toZ.inverse() 286 G4RotationMatrix toLabFrame3 = toZ.inverse();
303 G4LorentzRotation toLabFrame(toLabFrame3); 287 G4LorentzRotation toLabFrame(toLabFrame3);
304 // However, it turns out that the projectile 288 // However, it turns out that the projectile given to us by G4
305 // hadronic framework is already going in th 289 // hadronic framework is already going in the direction of the
306 // z-axis so this rotation is actually unnec 290 // z-axis so this rotation is actually unnecessary. Both toZ and
307 // toLabFrame turn out to be unit matrices a 291 // toLabFrame turn out to be unit matrices as can be seen by
308 // uncommenting the folowing two lines: 292 // uncommenting the folowing two lines:
309 // G4cout <<"toZ = " << toZ << G4endl; 293 // G4cout <<"toZ = " << toZ << G4endl;
310 // G4cout <<"toLabFrame = " << toLabFrame << 294 // G4cout <<"toLabFrame = " << toLabFrame << G4endl;
311 295
312 theResult.Clear(); // Make sure the output d 296 theResult.Clear(); // Make sure the output data structure is clean.
313 theResult.SetStatusChange(stopAndKill); 297 theResult.SetStatusChange(stopAndKill);
314 298
315 std::list<G4Fragment> remnants; 299 std::list<G4Fragment> remnants;
316 300
317 const G4int maxTries = 200; 301 const G4int maxTries = 200;
318 G4int nTries = 0; 302 G4int nTries = 0;
319 // INCL can generate transparent events. How 303 // INCL can generate transparent events. However, this is meaningful
320 // only in the standalone code. In Geant4 we 304 // only in the standalone code. In Geant4 we must "force" INCL to
321 // produce a valid cascade. 305 // produce a valid cascade.
322 G4bool eventIsOK = false; 306 G4bool eventIsOK = false;
323 do { 307 do {
324 const G4INCL::ParticleSpecies theSpecies = 308 const G4INCL::ParticleSpecies theSpecies = toINCLParticleSpecies(*aProjectileTrack);
325 const G4double kineticEnergy = toINCLKinet 309 const G4double kineticEnergy = toINCLKineticEnergy(*aProjectileTrack);
326 310
327 // The INCL model will be created at the f 311 // The INCL model will be created at the first use
328 theINCLModel = G4INCLXXInterfaceStore::Get 312 theINCLModel = G4INCLXXInterfaceStore::GetInstance()->GetINCLModel();
329 313
330 const G4INCL::EventInfo eventInfo = theINC << 314 const G4INCL::EventInfo eventInfo = theINCLModel->processEvent(theSpecies, kineticEnergy, theTargetNucleus->GetA_asInt(), theTargetNucleus->GetZ_asInt(),0);
331 theTargetNucleus->GetA_asIn << 315 // eventIsOK = !eventInfo.transparent && nTries < maxTries;
332 theTargetNucleus->GetZ_asIn <<
333 -theTargetNucleus->GetL()); <<
334 // eventIsOK = !eventInfo.transparent & <<
335 eventIsOK = !eventInfo.transparent; 316 eventIsOK = !eventInfo.transparent;
336 if(eventIsOK) { 317 if(eventIsOK) {
337 318
338 // If the collision was run in inverse k 319 // If the collision was run in inverse kinematics, prepare to boost back
339 // all the reaction products 320 // all the reaction products
340 if(inverseKinematics) { 321 if(inverseKinematics) {
341 toDirectKinematics = new G4LorentzRota 322 toDirectKinematics = new G4LorentzRotation(toInverseKinematics->inverse());
342 } 323 }
343 324
344 G4LorentzVector fourMomentumOut; 325 G4LorentzVector fourMomentumOut;
345 326
346 for(G4int i = 0; i < eventInfo.nParticle 327 for(G4int i = 0; i < eventInfo.nParticles; ++i) {
347 G4int A = eventInfo.A[i]; << 328 G4int A = eventInfo.A[i];
348 G4int Z = eventInfo.Z[i]; << 329 G4int Z = eventInfo.Z[i];
349 G4int S = eventInfo.S[i]; // Strangen << 330 G4int PDGCode = eventInfo.PDGCode[i];
350 G4int PDGCode = eventInfo.PDGCode[i]; << 331 // G4cout <<"INCL particle A = " << A << " Z = " << Z << G4endl;
351 G4double kinE = eventInfo.EKin[i]; << 332 G4double kinE = eventInfo.EKin[i];
352 G4double px = eventInfo.px[i]; << 333 G4double px = eventInfo.px[i];
353 G4double py = eventInfo.py[i]; << 334 G4double py = eventInfo.py[i];
354 G4double pz = eventInfo.pz[i]; << 335 G4double pz = eventInfo.pz[i];
355 G4DynamicParticle *p = toG4Particle(A, << 336 G4DynamicParticle *p = toG4Particle(A, Z, PDGCode, kinE, px, py, pz);
356 if(p != 0) { << 337 if(p != 0) {
357 G4LorentzVector momentum = p->Get4Mo << 338 G4LorentzVector momentum = p->Get4Momentum();
358 339
359 // Apply the toLabFrame rotation 340 // Apply the toLabFrame rotation
360 momentum *= toLabFrame; 341 momentum *= toLabFrame;
361 // Apply the inverse-kinematics boos 342 // Apply the inverse-kinematics boost
362 if(inverseKinematics) { 343 if(inverseKinematics) {
363 momentum *= *toDirectKinematics; 344 momentum *= *toDirectKinematics;
364 momentum.setVect(-momentum.vect()) 345 momentum.setVect(-momentum.vect());
365 } 346 }
366 347
367 // Set the four-momentum of the reaction p << 348 // Set the four-momentum of the reaction products
368 p->Set4Momentum(momentum); << 349 p->Set4Momentum(momentum);
369 fourMomentumOut += momentum; 350 fourMomentumOut += momentum;
>> 351 theResult.AddSecondary(p);
370 352
371 // Propagate the particle's parent resonan << 353 } else {
372 G4HadSecondary secondary(p, 1.0, sec << 354 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 355 theInterfaceStore->EmitWarning(message);
385 } << 356 }
386 } 357 }
387 358
388 for(G4int i = 0; i < eventInfo.nRemnants 359 for(G4int i = 0; i < eventInfo.nRemnants; ++i) {
389 const G4int A = eventInfo.ARem[i]; << 360 const G4int A = eventInfo.ARem[i];
390 const G4int Z = eventInfo.ZRem[i]; << 361 const G4int Z = eventInfo.ZRem[i];
391 const G4int S = eventInfo.SRem[i]; << 362 // G4cout <<"INCL particle A = " << A << " Z = " << Z << G4endl;
392 // Check that the remnant is a physica << 363 const G4double kinE = eventInfo.EKinRem[i];
393 if(( Z == 0 && S == 0 && A > 1 ) | << 364 const G4double px = eventInfo.pxRem[i];
394 ( Z == 0 && S != 0 && A < 4 ) | << 365 const G4double py = eventInfo.pyRem[i];
395 ( Z != 0 && S != 0 && A == Z + << 366 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( 367 G4ThreeVector spin(
408 eventInfo.jxRem[i]*hbar_Planck, 368 eventInfo.jxRem[i]*hbar_Planck,
409 eventInfo.jyRem[i]*hbar_Planck, 369 eventInfo.jyRem[i]*hbar_Planck,
410 eventInfo.jzRem[i]*hbar_Planck 370 eventInfo.jzRem[i]*hbar_Planck
411 ); 371 );
412 const G4double excitationE = eventInfo << 372 const G4double excitationE = eventInfo.EStarRem[i];
413 G4double nuclearMass = excitationE; << 373 const G4double nuclearMass = G4NucleiProperties::GetNuclearMass(A, Z) + excitationE;
414 if ( S == 0 ) { << 374 const G4double scaling = remnant4MomentumScaling(nuclearMass,
415 nuclearMass += G4NucleiProperties::G << 375 kinE,
416 } else { << 376 px, py, pz);
417 // Assumed that the opposite of the strang << 377 G4LorentzVector fourMomentum(scaling * px, scaling * py, scaling * pz,
418 nuclearMass += G4HyperNucleiProperti << 378 nuclearMass + kinE);
419 } << 379 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; 380 std::stringstream ss;
425 ss << "momentum scaling = " << scali 381 ss << "momentum scaling = " << scaling
426 << "\n Lorentz vec << 382 << "\n Lorentz vector = " << fourMomentum
427 << ")\n A = " << A << 383 << ")\n A = " << A << ", Z = " << Z
428 << "\n E* = " << e << 384 << "\n E* = " << excitationE << ", nuclearMass = " << nuclearMass
429 << "\n remnant i=" << 385 << "\n remnant i=" << i << ", nRemnants=" << eventInfo.nRemnants
430 << "\n Reaction wa << 386 << "\n Reaction was: " << aTrack.GetKineticEnergy()/MeV
431 << "-MeV " << trackDefinition->Ge << 387 << "-MeV " << trackDefinition->GetParticleName() << " + "
432 << theIonTable->GetIonName(theNuc << 388 << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
433 << ", in " << (inverseKinematics << 389 << ", in " << (inverseKinematics ? "inverse" : "direct") << " kinematics.";
434 theInterfaceStore->EmitWarning(ss.st 390 theInterfaceStore->EmitWarning(ss.str());
435 } << 391 }
436 392
437 // Apply the toLabFrame rotation 393 // Apply the toLabFrame rotation
438 fourMomentum *= toLabFrame; 394 fourMomentum *= toLabFrame;
439 spin *= toLabFrame3; 395 spin *= toLabFrame3;
440 // Apply the inverse-kinematics boost 396 // Apply the inverse-kinematics boost
441 if(inverseKinematics) { 397 if(inverseKinematics) {
442 fourMomentum *= *toDirectKinematics; 398 fourMomentum *= *toDirectKinematics;
443 fourMomentum.setVect(-fourMomentum.v 399 fourMomentum.setVect(-fourMomentum.vect());
444 } 400 }
445 401
446 fourMomentumOut += fourMomentum; 402 fourMomentumOut += fourMomentum;
447 G4Fragment remnant(A, Z, std::abs(S), << 403 G4Fragment remnant(A, Z, fourMomentum);
448 remnant.SetAngularMomentum(spin); 404 remnant.SetAngularMomentum(spin);
449 remnant.SetCreatorModelID(secID); <<
450 if(dumpRemnantInfo) { 405 if(dumpRemnantInfo) {
451 G4cerr << "G4INCLXX_DUMP_REMNANT: " 406 G4cerr << "G4INCLXX_DUMP_REMNANT: " << remnant << " spin: " << spin << G4endl;
452 } 407 }
453 remnants.push_back(remnant); << 408 remnants.push_back(remnant);
454 } 409 }
455 410
456 // Give up is the event is not ok (e.g. << 411 // Check four-momentum conservation
457 if(!eventIsOK) { << 412 const G4LorentzVector violation4Momentum = fourMomentumOut - fourMomentumIn;
458 const G4int nSecondaries = (G4int)theR << 413 const G4double energyViolation = std::abs(violation4Momentum.e());
459 for(G4int j=0; j<nSecondaries; ++j) de << 414 const G4double momentumViolation = violation4Momentum.rho();
>> 415 if(energyViolation > G4INCLXXInterfaceStore::GetInstance()->GetConservationTolerance()) {
>> 416 std::stringstream ss;
>> 417 ss << "energy conservation violated by " << energyViolation/MeV << " MeV in "
>> 418 << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
>> 419 << " + " << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
>> 420 << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics. Will resample.";
>> 421 theInterfaceStore->EmitWarning(ss.str());
>> 422 eventIsOK = false;
>> 423 const G4int nSecondaries = theResult.GetNumberOfSecondaries();
>> 424 for(G4int j=0; j<nSecondaries; ++j)
>> 425 delete theResult.GetSecondary(j)->GetParticle();
>> 426 theResult.Clear();
>> 427 theResult.SetStatusChange(stopAndKill);
>> 428 remnants.clear();
>> 429 } else if(momentumViolation > G4INCLXXInterfaceStore::GetInstance()->GetConservationTolerance()) {
>> 430 std::stringstream ss;
>> 431 ss << "momentum conservation violated by " << momentumViolation/MeV << " MeV in "
>> 432 << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
>> 433 << " + " << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
>> 434 << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics. Will resample.";
>> 435 theInterfaceStore->EmitWarning(ss.str());
>> 436 eventIsOK = false;
>> 437 const G4int nSecondaries = theResult.GetNumberOfSecondaries();
>> 438 for(G4int j=0; j<nSecondaries; ++j)
>> 439 delete theResult.GetSecondary(j)->GetParticle();
460 theResult.Clear(); 440 theResult.Clear();
461 theResult.SetStatusChange(stopAndKill) 441 theResult.SetStatusChange(stopAndKill);
462 remnants.clear(); 442 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 } 443 }
496 } 444 }
497 nTries++; 445 nTries++;
498 } while(!eventIsOK && nTries < maxTries); /* 446 } while(!eventIsOK && nTries < maxTries); /* Loop checking, 10.07.2015, D.Mancusi */
499 447
500 // Clean up the objects that we created for 448 // Clean up the objects that we created for the inverse kinematics
501 if(inverseKinematics) { 449 if(inverseKinematics) {
502 delete aProjectileTrack; 450 delete aProjectileTrack;
503 delete theTargetNucleus; 451 delete theTargetNucleus;
504 delete toInverseKinematics; 452 delete toInverseKinematics;
505 delete toDirectKinematics; 453 delete toDirectKinematics;
506 } 454 }
507 455
508 if(!eventIsOK) { 456 if(!eventIsOK) {
509 std::stringstream ss; 457 std::stringstream ss;
510 ss << "maximum number of tries exceeded fo 458 ss << "maximum number of tries exceeded for the proposed "
511 << aTrack.GetKineticEnergy()/MeV << "-MeV 459 << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
512 << " + " << theIonTable->GetIonName(nucleu 460 << " + " << theIonTable->GetIonName(nucleusZ, nucleusA, 0)
513 << " inelastic reaction, in " << (inverseK 461 << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics.";
514 theInterfaceStore->EmitWarning(ss.str()); 462 theInterfaceStore->EmitWarning(ss.str());
515 theResult.SetStatusChange(isAlive); 463 theResult.SetStatusChange(isAlive);
516 theResult.SetEnergyChange(aTrack.GetKineti 464 theResult.SetEnergyChange(aTrack.GetKineticEnergy());
517 theResult.SetMomentumChange(aTrack.Get4Mom 465 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
518 return &theResult; 466 return &theResult;
519 } 467 }
520 468
521 // De-excitation: 469 // De-excitation:
522 470
523 if(theDeExcitation != 0) { 471 if(theDeExcitation != 0) {
524 for(std::list<G4Fragment>::iterator i = re 472 for(std::list<G4Fragment>::iterator i = remnants.begin();
525 i != remnants.end(); ++i) { << 473 i != remnants.end(); ++i) {
526 G4ReactionProductVector *deExcitationRes 474 G4ReactionProductVector *deExcitationResult = theDeExcitation->DeExcite((*i));
527 475
528 for(G4ReactionProductVector::iterator fr 476 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
529 fragment != deExcitationResult->end(); ++f << 477 fragment != deExcitationResult->end(); ++fragment) {
530 const G4ParticleDefinition *def = (*fragment << 478 const G4ParticleDefinition *def = (*fragment)->GetDefinition();
531 if(def != 0) { << 479 if(def != 0) {
532 G4DynamicParticle *theFragment = new G4Dyn << 480 G4DynamicParticle *theFragment = new G4DynamicParticle(def, (*fragment)->GetMomentum());
533 theResult.AddSecondary(theFragment, (*frag << 481 theResult.AddSecondary(theFragment);
534 } << 482 }
535 } 483 }
536 484
537 for(G4ReactionProductVector::iterator fr 485 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
538 fragment != deExcitationResult->end(); ++f << 486 fragment != deExcitationResult->end(); ++fragment) {
539 delete (*fragment); << 487 delete (*fragment);
540 } 488 }
541 deExcitationResult->clear(); 489 deExcitationResult->clear();
542 delete deExcitationResult; 490 delete deExcitationResult;
543 } 491 }
544 } 492 }
545 493
546 remnants.clear(); 494 remnants.clear();
547 495
548 if((theTally = theInterfaceStore->GetTally() 496 if((theTally = theInterfaceStore->GetTally()))
549 theTally->Tally(aTrack, theNucleus, theRes 497 theTally->Tally(aTrack, theNucleus, theResult);
550 498
551 return &theResult; 499 return &theResult;
552 } 500 }
553 501
554 G4ReactionProductVector* G4INCLXXInterface::Pr 502 G4ReactionProductVector* G4INCLXXInterface::Propagate(G4KineticTrackVector* , G4V3DNucleus* ) {
555 return 0; 503 return 0;
556 } 504 }
557 505
558 G4INCL::ParticleType G4INCLXXInterface::toINCL 506 G4INCL::ParticleType G4INCLXXInterface::toINCLParticleType(G4ParticleDefinition const * const pdef) const {
559 if( pdef == G4Proton::Proton()) << 507 if( pdef == G4Proton::Proton()) return G4INCL::Proton;
560 else if(pdef == G4Neutron::Neutron()) << 508 else if(pdef == G4Neutron::Neutron()) return G4INCL::Neutron;
561 else if(pdef == G4PionPlus::PionPlus()) << 509 else if(pdef == G4PionPlus::PionPlus()) return G4INCL::PiPlus;
562 else if(pdef == G4PionMinus::PionMinus()) << 510 else if(pdef == G4PionMinus::PionMinus()) return G4INCL::PiMinus;
563 else if(pdef == G4PionZero::PionZero()) << 511 else if(pdef == G4PionZero::PionZero()) return G4INCL::PiZero;
564 else if(pdef == G4KaonPlus::KaonPlus()) << 512 else if(pdef == G4KaonPlus::KaonPlus()) return G4INCL::KPlus;
565 else if(pdef == G4KaonZero::KaonZero()) << 513 else if(pdef == G4KaonMinus::KaonMinus()) return G4INCL::KMinus;
566 else if(pdef == G4KaonMinus::KaonMinus()) << 514 else if(pdef == G4Deuteron::Deuteron()) return G4INCL::Composite;
567 else if(pdef == G4AntiKaonZero::AntiKaonZero << 515 else if(pdef == G4Triton::Triton()) return G4INCL::Composite;
568 // For K0L & K0S we do not take into account << 516 else if(pdef == G4He3::He3()) return G4INCL::Composite;
569 else if(pdef == G4KaonZeroLong::KaonZeroLong << 517 else if(pdef == G4Alpha::Alpha()) return G4INCL::Composite;
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 518 else if(pdef->GetParticleType() == G4GenericIon::GenericIon()->GetParticleType()) return G4INCL::Composite;
577 else << 519 else return G4INCL::UnknownParticle;
578 } 520 }
579 521
580 G4INCL::ParticleSpecies G4INCLXXInterface::toI 522 G4INCL::ParticleSpecies G4INCLXXInterface::toINCLParticleSpecies(G4HadProjectile const &aTrack) const {
581 const G4ParticleDefinition *pdef = aTrack.Ge 523 const G4ParticleDefinition *pdef = aTrack.GetDefinition();
582 const G4INCL::ParticleType theType = toINCLP 524 const G4INCL::ParticleType theType = toINCLParticleType(pdef);
583 if(theType!=G4INCL::Composite) 525 if(theType!=G4INCL::Composite)
584 return G4INCL::ParticleSpecies(theType); 526 return G4INCL::ParticleSpecies(theType);
585 else { 527 else {
586 G4INCL::ParticleSpecies theSpecies; 528 G4INCL::ParticleSpecies theSpecies;
587 theSpecies.theType=theType; 529 theSpecies.theType=theType;
588 theSpecies.theA=pdef->GetAtomicMass(); 530 theSpecies.theA=pdef->GetAtomicMass();
589 theSpecies.theZ=pdef->GetAtomicNumber(); 531 theSpecies.theZ=pdef->GetAtomicNumber();
590 return theSpecies; 532 return theSpecies;
591 } 533 }
592 } 534 }
593 535
594 G4double G4INCLXXInterface::toINCLKineticEnerg 536 G4double G4INCLXXInterface::toINCLKineticEnergy(G4HadProjectile const &aTrack) const {
595 return aTrack.GetKineticEnergy(); 537 return aTrack.GetKineticEnergy();
596 } 538 }
597 539
598 G4ParticleDefinition *G4INCLXXInterface::toG4P << 540 G4ParticleDefinition *G4INCLXXInterface::toG4ParticleDefinition(G4int A, G4int Z, G4int PDGCode) const {
599 if (PDGCode == 2212) { return G4Proton << 541 if(PDGCode == 2212) return G4Proton::Proton();
600 } else if(PDGCode == 2112) { return G4Neutro << 542 else if(PDGCode == 2112) return G4Neutron::Neutron();
601 } else if(PDGCode == 211) { return G4PionPl << 543 else if(PDGCode == 211) return G4PionPlus::PionPlus();
602 } else if(PDGCode == 111) { return G4PionZe << 544 else if(PDGCode == 111) return G4PionZero::PionZero();
603 } else if(PDGCode == -211) { return G4PionMi << 545 else if(PDGCode == -211) return G4PionMinus::PionMinus();
604 546
605 } else if(PDGCode == 221) { return G4Eta::E << 547 else if(PDGCode == 221) return G4Eta::Eta();
606 } else if(PDGCode == 22) { return G4Gamma: << 548 else if(PDGCode == 22) return G4Gamma::Gamma();
607 549
608 } else if(PDGCode == 3122) { return G4Lambda << 550 else if(PDGCode == 3122) return G4Lambda::Lambda();
609 } else if(PDGCode == 3222) { return G4SigmaP << 551 else if(PDGCode == 3222) return G4SigmaPlus::SigmaPlus();
610 } else if(PDGCode == 3212) { return G4SigmaZ << 552 else if(PDGCode == 3212) return G4SigmaZero::SigmaZero();
611 } else if(PDGCode == 3112) { return G4SigmaM << 553 else if(PDGCode == 3112) return G4SigmaMinus::SigmaMinus();
612 } else if(PDGCode == 321) { return G4KaonPl << 554 else if(PDGCode == 321) return G4KaonPlus::KaonPlus();
613 } else if(PDGCode == -321) { return G4KaonMi << 555 else if(PDGCode == -321) return G4KaonMinus::KaonMinus();
614 } else if(PDGCode == 130) { return G4KaonZe << 556 else if(PDGCode == 130) return G4KaonZeroLong::KaonZeroLong();
615 } else if(PDGCode == 310) { return G4KaonZe << 557 else if(PDGCode == 310) return G4KaonZeroShort::KaonZeroShort();
616 558
617 } else if(PDGCode == 1002) { return G4Deuter << 559 else if(PDGCode == 1002) return G4Deuteron::Deuteron();
618 } else if(PDGCode == 1003) { return G4Triton << 560 else if(PDGCode == 1003) return G4Triton::Triton();
619 } else if(PDGCode == 2003) { return G4He3::H << 561 else if(PDGCode == 2003) return G4He3::He3();
620 } else if(PDGCode == 2004) { return G4Alpha: << 562 else if(PDGCode == 2004) return G4Alpha::Alpha();
621 << 563 else if(A > 0 && Z > 0 && A > Z) { // Returns ground state ion definition. No hyper-nucleus allows in Geant4
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); 564 return theIonTable->GetIon(Z, A, 0);
>> 565 } else { // Error, unrecognized particle
>> 566 return 0;
632 } 567 }
633 return 0; // Error, unrecognized particle <<
634 } 568 }
635 569
636 G4DynamicParticle *G4INCLXXInterface::toG4Part << 570 G4DynamicParticle *G4INCLXXInterface::toG4Particle(G4int A, G4int Z, G4int PDGCode,
637 G4double kinE, G4double px, << 571 G4double kinE,
638 << 572 G4double px,
639 const G4ParticleDefinition *def = toG4Partic << 573 G4double py, G4double pz) const {
>> 574 const G4ParticleDefinition *def = toG4ParticleDefinition(A, Z, PDGCode);
640 if(def == 0) { // Check if we have a valid p 575 if(def == 0) { // Check if we have a valid particle definition
641 return 0; 576 return 0;
642 } 577 }
643 const G4double energy = kinE * MeV; 578 const G4double energy = kinE * MeV;
644 const G4ThreeVector momentum(px, py, pz); 579 const G4ThreeVector momentum(px, py, pz);
645 const G4ThreeVector momentumDirection = mome 580 const G4ThreeVector momentumDirection = momentum.unit();
646 G4DynamicParticle *p = new G4DynamicParticle 581 G4DynamicParticle *p = new G4DynamicParticle(def, momentumDirection, energy);
647 return p; 582 return p;
648 } 583 }
649 584
650 G4double G4INCLXXInterface::remnant4MomentumSc 585 G4double G4INCLXXInterface::remnant4MomentumScaling(G4double mass,
651 G4double kineticE, << 586 G4double kineticE,
652 G4double px, G4double py, << 587 G4double px, G4double py,
653 G4double pz) const { << 588 G4double pz) const {
654 const G4double p2 = px*px + py*py + pz*pz; 589 const G4double p2 = px*px + py*py + pz*pz;
655 if(p2 > 0.0) { 590 if(p2 > 0.0) {
656 const G4double pnew2 = kineticE*kineticE + 591 const G4double pnew2 = kineticE*kineticE + 2.0*kineticE*mass;
657 return std::sqrt(pnew2)/std::sqrt(p2); 592 return std::sqrt(pnew2)/std::sqrt(p2);
658 } else { 593 } else {
659 return 1.0; 594 return 1.0;
660 } 595 }
661 } 596 }
662 597
663 void G4INCLXXInterface::ModelDescription(std:: 598 void G4INCLXXInterface::ModelDescription(std::ostream& outFile) const {
664 outFile 599 outFile
665 << "The Liège Intranuclear Cascade (INCL << 600 << "The Li�ge Intranuclear Cascade (INCL++) is a model for reactions induced\n"
666 << "by nucleons, pions and light ion on a 601 << "by nucleons, pions and light ion on any nucleus. The reaction is\n"
667 << "described as an avalanche of binary n 602 << "described as an avalanche of binary nucleon-nucleon collisions, which can\n"
668 << "lead to the emission of energetic par 603 << "lead to the emission of energetic particles and to the formation of an\n"
669 << "excited thermalised nucleus (remnant) 604 << "excited thermalised nucleus (remnant). The de-excitation of the remnant is\n"
670 << "outside the scope of INCL++ and is ty 605 << "outside the scope of INCL++ and is typically described by another model.\n\n"
671 << "INCL++ has been reasonably well teste 606 << "INCL++ has been reasonably well tested for nucleon (~50 MeV to ~15 GeV),\n"
672 << "pion (idem) and light-ion projectiles 607 << "pion (idem) and light-ion projectiles (up to A=18, ~10A MeV to 1A GeV).\n"
673 << "Most tests involved target nuclei clo 608 << "Most tests involved target nuclei close to the stability valley, with\n"
674 << "numbers between 4 and 250.\n\n" 609 << "numbers between 4 and 250.\n\n"
675 << "Reference: D. Mancusi et al., Phys. R 610 << "Reference: D. Mancusi et al., Phys. Rev. C90 (2014) 054602\n\n";
676 } 611 }
677 612
678 G4String const &G4INCLXXInterface::GetDeExcita 613 G4String const &G4INCLXXInterface::GetDeExcitationModelName() const {
679 return theDeExcitation->GetModelName(); 614 return theDeExcitation->GetModelName();
680 } 615 }
681 616