| 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 // * 20 // * *
21 // * Parts of this code which have been devel 21 // * Parts of this code which have been developed by QinetiQ Ltd *
22 // * under contract to the European Space Agen 22 // * under contract to the European Space Agency (ESA) are the *
23 // * intellectual property of ESA. Rights to u 23 // * intellectual property of ESA. Rights to use, copy, modify and *
24 // * redistribute this software for general pu 24 // * redistribute this software for general public use are granted *
25 // * in compliance with any licensing, distrib 25 // * in compliance with any licensing, distribution and development *
26 // * policy adopted by the Geant4 Collaboratio 26 // * policy adopted by the Geant4 Collaboration. This code has been *
27 // * written by QinetiQ Ltd for the European S 27 // * written by QinetiQ Ltd for the European Space Agency, under ESA *
28 // * contract 17191/03/NL/LvH (Aurora Programm 28 // * contract 17191/03/NL/LvH (Aurora Programme). *
29 // * 29 // * *
30 // * By using, copying, modifying or distri 30 // * By using, copying, modifying or distributing the software (or *
31 // * any work based on the software) you ag 31 // * any work based on the software) you agree to acknowledge its *
32 // * use in resulting scientific publicati 32 // * use in resulting scientific publications, and indicate your *
33 // * acceptance of all terms of the Geant4 Sof 33 // * acceptance of all terms of the Geant4 Software license. *
34 // ******************************************* 34 // ********************************************************************
35 // 35 //
36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37 // 37 //
38 // MODULE: G4WilsonAbrasionModel. 38 // MODULE: G4WilsonAbrasionModel.cc
39 // 39 //
40 // Version: 1.0 40 // Version: 1.0
41 // Date: 08/12/2009 41 // Date: 08/12/2009
42 // Author: P R Truscott 42 // Author: P R Truscott
43 // Organisation: QinetiQ Ltd, UK 43 // Organisation: QinetiQ Ltd, UK
44 // Customer: ESA/ESTEC, NOORDWIJK 44 // Customer: ESA/ESTEC, NOORDWIJK
45 // Contract: 17191/03/NL/LvH 45 // Contract: 17191/03/NL/LvH
46 // 46 //
47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48 // 48 //
49 // CHANGE HISTORY 49 // CHANGE HISTORY
50 // -------------- 50 // --------------
51 // 51 //
52 // 6 October 2003, P R Truscott, QinetiQ Ltd, 52 // 6 October 2003, P R Truscott, QinetiQ Ltd, UK
53 // Created. 53 // Created.
54 // 54 //
55 // 15 March 2004, P R Truscott, QinetiQ Ltd, U 55 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK
56 // Beta release 56 // Beta release
57 // 57 //
58 // 18 January 2005, M H Mendenhall, Vanderbilt 58 // 18 January 2005, M H Mendenhall, Vanderbilt University, US
59 // Pointers to theAbrasionGeometry and product 59 // Pointers to theAbrasionGeometry and products generated by the deexcitation
60 // handler deleted to prevent memory leaks. A 60 // handler deleted to prevent memory leaks. Also particle change of projectile
61 // fragment previously not properly defined. 61 // fragment previously not properly defined.
62 // 62 //
63 // 08 December 2009, P R Truscott, QinetiQ Ltd 63 // 08 December 2009, P R Truscott, QinetiQ Ltd, Ltd
64 // ver 1.0 64 // ver 1.0
65 // There was originally a possibility of the m 65 // There was originally a possibility of the minimum impact parameter AFTER
66 // considering Couloumb repulsion, rm, being t 66 // considering Couloumb repulsion, rm, being too large. Now if:
67 // rm >= fradius * (rP + rT) 67 // rm >= fradius * (rP + rT)
68 // where fradius is currently 0.99, then it is 68 // where fradius is currently 0.99, then it is assumed the primary track is
69 // unchanged. Additional conditions to escape 69 // unchanged. Additional conditions to escape from while-loop over impact
70 // parameter: if the loop counter evtcnt reach 70 // parameter: if the loop counter evtcnt reaches 1000, the primary track is
71 // continued. 71 // continued.
72 // Additional clauses have been included in 72 // Additional clauses have been included in
73 // G4WilsonAbrasionModel::GetNucleonInduced 73 // G4WilsonAbrasionModel::GetNucleonInducedExcitation
74 // Previously it was possible to get sqrt of n 74 // Previously it was possible to get sqrt of negative number as Wilson
75 // algorithm not properly defined if either: 75 // algorithm not properly defined if either:
76 // rT > rP && rsq < rTsq - rPsq) or (rP > r 76 // rT > rP && rsq < rTsq - rPsq) or (rP > rT && rsq < rPsq - rTsq)
77 // 77 //
78 // 12 June 2012, A. Ribon, CERN, Switzerland 78 // 12 June 2012, A. Ribon, CERN, Switzerland
79 // Fixing trivial warning errors of shadowed v 79 // Fixing trivial warning errors of shadowed variables.
80 // 80 //
81 // 4 August 2015, A. Ribon, CERN, Switzerland 81 // 4 August 2015, A. Ribon, CERN, Switzerland
82 // Replacing std::exp and std::pow with the fa 82 // Replacing std::exp and std::pow with the faster versions G4Exp and G4Pow.
83 // 83 //
84 // 7 August 2015, A. Ribon, CERN, Switzerland 84 // 7 August 2015, A. Ribon, CERN, Switzerland
85 // Checking of 'while' loops. 85 // Checking of 'while' loops.
86 // 86 //
87 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 87 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
88 ////////////////////////////////////////////// 88 ///////////////////////////////////////////////////////////////////////////////
89 89
90 #include "G4WilsonAbrasionModel.hh" 90 #include "G4WilsonAbrasionModel.hh"
91 #include "G4WilsonRadius.hh" 91 #include "G4WilsonRadius.hh"
92 #include "G4NuclearAbrasionGeometry.hh" 92 #include "G4NuclearAbrasionGeometry.hh"
93 #include "G4WilsonAblationModel.hh" 93 #include "G4WilsonAblationModel.hh"
94 94
95 #include "G4PhysicalConstants.hh" 95 #include "G4PhysicalConstants.hh"
96 #include "G4SystemOfUnits.hh" 96 #include "G4SystemOfUnits.hh"
97 #include "G4ExcitationHandler.hh" 97 #include "G4ExcitationHandler.hh"
98 #include "G4Evaporation.hh" 98 #include "G4Evaporation.hh"
>> 99 #include "G4FermiBreakUp.hh"
>> 100 #include "G4StatMF.hh"
99 #include "G4ParticleDefinition.hh" 101 #include "G4ParticleDefinition.hh"
100 #include "G4DynamicParticle.hh" 102 #include "G4DynamicParticle.hh"
101 #include "Randomize.hh" 103 #include "Randomize.hh"
102 #include "G4Fragment.hh" 104 #include "G4Fragment.hh"
103 #include "G4ReactionProductVector.hh" 105 #include "G4ReactionProductVector.hh"
104 #include "G4LorentzVector.hh" 106 #include "G4LorentzVector.hh"
105 #include "G4ParticleMomentum.hh" 107 #include "G4ParticleMomentum.hh"
106 #include "G4Poisson.hh" 108 #include "G4Poisson.hh"
107 #include "G4ParticleTable.hh" 109 #include "G4ParticleTable.hh"
108 #include "G4IonTable.hh" 110 #include "G4IonTable.hh"
109 #include "globals.hh" 111 #include "globals.hh"
110 112
111 #include "G4Exp.hh" 113 #include "G4Exp.hh"
112 #include "G4Pow.hh" 114 #include "G4Pow.hh"
113 115
114 #include "G4PhysicsModelCatalog.hh" <<
115 <<
116 116
117 G4WilsonAbrasionModel::G4WilsonAbrasionModel(G 117 G4WilsonAbrasionModel::G4WilsonAbrasionModel(G4bool useAblation1)
118 : G4HadronicInteraction("G4WilsonAbrasion"), << 118 :G4HadronicInteraction("G4WilsonAbrasion")
119 { 119 {
120 // Send message to stdout to advise that the 120 // Send message to stdout to advise that the G4Abrasion model is being used.
121 PrintWelcomeMessage(); 121 PrintWelcomeMessage();
122 122
123 // Set the default verbose level to 0 - no o 123 // Set the default verbose level to 0 - no output.
124 verboseLevel = 0; 124 verboseLevel = 0;
125 useAblation = useAblation1; 125 useAblation = useAblation1;
126 theAblation = nullptr; <<
127 126
128 // No de-excitation handler has been supplie 127 // No de-excitation handler has been supplied - define the default handler.
129 128
130 theExcitationHandler = new G4ExcitationHandl << 129 theExcitationHandler = new G4ExcitationHandler;
>> 130 theExcitationHandlerx = new G4ExcitationHandler;
131 if (useAblation) 131 if (useAblation)
132 { 132 {
133 theAblation = new G4WilsonAblationModel; 133 theAblation = new G4WilsonAblationModel;
134 theAblation->SetVerboseLevel(verboseLevel) 134 theAblation->SetVerboseLevel(verboseLevel);
135 theExcitationHandler->SetEvaporation(theAb << 135 theExcitationHandler->SetEvaporation(theAblation);
>> 136 theExcitationHandlerx->SetEvaporation(theAblation);
>> 137 }
>> 138 else
>> 139 {
>> 140 theAblation = NULL;
>> 141 G4Evaporation * theEvaporation = new G4Evaporation;
>> 142 G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp;
>> 143 G4StatMF * theMF = new G4StatMF;
>> 144 theExcitationHandler->SetEvaporation(theEvaporation);
>> 145 theExcitationHandler->SetFermiModel(theFermiBreakUp);
>> 146 theExcitationHandler->SetMultiFragmentation(theMF);
>> 147 theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6);
>> 148 theExcitationHandler->SetMinEForMultiFrag(5.0*MeV);
>> 149
>> 150 theEvaporation = new G4Evaporation;
>> 151 theFermiBreakUp = new G4FermiBreakUp;
>> 152 theExcitationHandlerx->SetEvaporation(theEvaporation);
>> 153 theExcitationHandlerx->SetFermiModel(theFermiBreakUp);
>> 154 theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6);
136 } 155 }
137 156
138 // Set the minimum and maximum range for the 157 // Set the minimum and maximum range for the model (despite nomanclature,
139 // this is in energy per nucleon number). 158 // this is in energy per nucleon number).
140 159
141 SetMinEnergy(70.0*MeV); 160 SetMinEnergy(70.0*MeV);
142 SetMaxEnergy(10.1*GeV); 161 SetMaxEnergy(10.1*GeV);
143 isBlocked = false; 162 isBlocked = false;
144 163
145 // npK, when mutiplied by the nuclear Fermi 164 // npK, when mutiplied by the nuclear Fermi momentum, determines the range of
146 // momentum over which the secondary nucleon 165 // momentum over which the secondary nucleon momentum is sampled.
147 166
148 r0sq = 0.0; 167 r0sq = 0.0;
149 npK = 5.0; 168 npK = 5.0;
150 B = 10.0 * MeV; 169 B = 10.0 * MeV;
151 third = 1.0 / 3.0; 170 third = 1.0 / 3.0;
152 fradius = 0.99; 171 fradius = 0.99;
153 conserveEnergy = false; 172 conserveEnergy = false;
154 conserveMomentum = true; 173 conserveMomentum = true;
155 <<
156 // Creator model ID for the secondaries crea <<
157 secID = G4PhysicsModelCatalog::GetModelID( " <<
158 } 174 }
159 175
160 void G4WilsonAbrasionModel::ModelDescription(s 176 void G4WilsonAbrasionModel::ModelDescription(std::ostream& outFile) const
161 { 177 {
162 outFile << "G4WilsonAbrasionModel is a macro 178 outFile << "G4WilsonAbrasionModel is a macroscopic treatment of\n"
163 << "nucleus-nucleus collisions using 179 << "nucleus-nucleus collisions using simple geometric arguments.\n"
164 << "The smaller projectile nucleus g 180 << "The smaller projectile nucleus gouges out a part of the larger\n"
165 << "target nucleus, leaving a residu 181 << "target nucleus, leaving a residual nucleus and a fireball\n"
166 << "region where the projectile and 182 << "region where the projectile and target intersect. The fireball"
167 << "is then treated as a highly exci 183 << "is then treated as a highly excited nuclear fragment. This\n"
168 << "model is based on the NUCFRG2 mo 184 << "model is based on the NUCFRG2 model and is valid for all\n"
169 << "projectile energies between 70 M 185 << "projectile energies between 70 MeV/n and 10.1 GeV/n. \n";
170 } 186 }
171 187
172 G4WilsonAbrasionModel::G4WilsonAbrasionModel(G << 188 G4WilsonAbrasionModel::G4WilsonAbrasionModel(G4ExcitationHandler* aExcitationHandler)
173 G4HadronicInteraction("G4WilsonAbrasion"), se <<
174 { 189 {
175 // Send message to stdout to advise that the G 190 // Send message to stdout to advise that the G4Abrasion model is being used.
176 191
177 PrintWelcomeMessage(); 192 PrintWelcomeMessage();
178 193
179 // Set the default verbose level to 0 - no out 194 // Set the default verbose level to 0 - no output.
180 195
181 verboseLevel = 0; 196 verboseLevel = 0;
182 197
183 theAblation = nullptr; //A.R. 26-Jul-2012 << 198 theAblation = NULL; //A.R. 26-Jul-2012 Coverity fix.
184 useAblation = false; //A.R. 14-Aug-2012 << 199 useAblation = false; //A.R. 14-Aug-2012 Coverity fix.
185 200
186 // 201 //
187 // The user is able to provide the excitation 202 // The user is able to provide the excitation handler as well as an argument
188 // which is provided in this instantiation is 203 // which is provided in this instantiation is used to determine
189 // whether the spectators of the interaction a 204 // whether the spectators of the interaction are free following the abrasion.
190 // 205 //
191 theExcitationHandler = aExcitationHandler; << 206 theExcitationHandler = aExcitationHandler;
>> 207 theExcitationHandlerx = new G4ExcitationHandler;
>> 208 G4Evaporation * theEvaporation = new G4Evaporation;
>> 209 G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp;
>> 210 theExcitationHandlerx->SetEvaporation(theEvaporation);
>> 211 theExcitationHandlerx->SetFermiModel(theFermiBreakUp);
>> 212 theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6);
192 // 213 //
193 // 214 //
194 // Set the minimum and maximum range for the m 215 // Set the minimum and maximum range for the model (despite nomanclature, this
195 // is in energy per nucleon number). 216 // is in energy per nucleon number).
196 // 217 //
197 SetMinEnergy(70.0*MeV); 218 SetMinEnergy(70.0*MeV);
198 SetMaxEnergy(10.1*GeV); 219 SetMaxEnergy(10.1*GeV);
199 isBlocked = false; 220 isBlocked = false;
200 // 221 //
201 // 222 //
202 // npK, when mutiplied by the nuclear Fermi mo 223 // npK, when mutiplied by the nuclear Fermi momentum, determines the range of
203 // momentum over which the secondary nucleon m 224 // momentum over which the secondary nucleon momentum is sampled.
204 // 225 //
205 r0sq = 0.0; << 226 r0sq = 0.0; //A.R. 14-Aug-2012 Coverity fix.
206 npK = 5.0; 227 npK = 5.0;
207 B = 10.0 * MeV; 228 B = 10.0 * MeV;
208 third = 1.0 / 3.0; 229 third = 1.0 / 3.0;
209 fradius = 0.99; 230 fradius = 0.99;
210 conserveEnergy = false; 231 conserveEnergy = false;
211 conserveMomentum = true; 232 conserveMomentum = true;
212 <<
213 // Creator model ID for the secondaries crea <<
214 secID = G4PhysicsModelCatalog::GetModelID( " <<
215 } 233 }
216 ////////////////////////////////////////////// 234 ////////////////////////////////////////////////////////////////////////////////
217 // 235 //
218 G4WilsonAbrasionModel::~G4WilsonAbrasionModel( << 236 G4WilsonAbrasionModel::~G4WilsonAbrasionModel ()
219 { 237 {
220 delete theExcitationHandler; << 238 //
>> 239 //
>> 240 // The destructor doesn't have to do a great deal!
>> 241 //
>> 242 if (theExcitationHandler) delete theExcitationHandler;
>> 243 if (theExcitationHandlerx) delete theExcitationHandlerx;
>> 244 if (useAblation) delete theAblation;
>> 245 // delete theExcitationHandler;
>> 246 // delete theExcitationHandlerx;
221 } 247 }
222 ////////////////////////////////////////////// 248 ////////////////////////////////////////////////////////////////////////////////
223 // 249 //
224 G4HadFinalState *G4WilsonAbrasionModel::ApplyY 250 G4HadFinalState *G4WilsonAbrasionModel::ApplyYourself (
225 const G4HadProjectile &theTrack, G4Nucleus & 251 const G4HadProjectile &theTrack, G4Nucleus &theTarget)
226 { 252 {
227 // 253 //
228 // 254 //
229 // The secondaries will be returned in G4HadFi 255 // The secondaries will be returned in G4HadFinalState &theParticleChange -
230 // initialise this. The original track will a 256 // initialise this. The original track will always be discontinued and
231 // secondaries followed. 257 // secondaries followed.
232 // 258 //
233 theParticleChange.Clear(); 259 theParticleChange.Clear();
234 theParticleChange.SetStatusChange(stopAndKil 260 theParticleChange.SetStatusChange(stopAndKill);
235 // 261 //
236 // 262 //
237 // Get relevant information about the projecti 263 // Get relevant information about the projectile and target (A, Z, energy/nuc,
238 // momentum, etc). 264 // momentum, etc).
239 // 265 //
240 const G4ParticleDefinition *definitionP = th 266 const G4ParticleDefinition *definitionP = theTrack.GetDefinition();
241 const G4double AP = definitionP->GetBaryonN 267 const G4double AP = definitionP->GetBaryonNumber();
242 const G4double ZP = definitionP->GetPDGChar 268 const G4double ZP = definitionP->GetPDGCharge();
243 G4LorentzVector pP = theTrack.Get4Momentum() 269 G4LorentzVector pP = theTrack.Get4Momentum();
244 G4double E = theTrack.GetKineticEner 270 G4double E = theTrack.GetKineticEnergy()/AP;
245 G4double AT = theTarget.GetA_asInt(); 271 G4double AT = theTarget.GetA_asInt();
246 G4double ZT = theTarget.GetZ_asInt(); 272 G4double ZT = theTarget.GetZ_asInt();
247 G4double TotalEPre = theTrack.GetTotalEnergy 273 G4double TotalEPre = theTrack.GetTotalEnergy() +
248 theTarget.AtomicMass(AT, ZT) + theTarget.G 274 theTarget.AtomicMass(AT, ZT) + theTarget.GetEnergyDeposit();
249 G4double TotalEPost = 0.0; 275 G4double TotalEPost = 0.0;
250 // 276 //
251 // 277 //
252 // Determine the radii of the projectile and t 278 // Determine the radii of the projectile and target nuclei.
253 // 279 //
254 G4WilsonRadius aR; 280 G4WilsonRadius aR;
255 G4double rP = aR.GetWilsonRadius(AP); 281 G4double rP = aR.GetWilsonRadius(AP);
256 G4double rT = aR.GetWilsonRadius(AT); 282 G4double rT = aR.GetWilsonRadius(AT);
257 G4double rPsq = rP * rP; 283 G4double rPsq = rP * rP;
258 G4double rTsq = rT * rT; 284 G4double rTsq = rT * rT;
259 if (verboseLevel >= 2) 285 if (verboseLevel >= 2)
260 { 286 {
261 G4cout <<"################################ 287 G4cout <<"########################################"
262 <<"################################ 288 <<"########################################"
263 <<G4endl; 289 <<G4endl;
264 G4cout.precision(6); 290 G4cout.precision(6);
265 G4cout <<"IN G4WilsonAbrasionModel" <<G4en 291 G4cout <<"IN G4WilsonAbrasionModel" <<G4endl;
266 G4cout <<"Initial projectile A=" <<AP 292 G4cout <<"Initial projectile A=" <<AP
267 <<", Z=" <<ZP 293 <<", Z=" <<ZP
268 <<", radius = " <<rP/fermi <<" fm" 294 <<", radius = " <<rP/fermi <<" fm"
269 <<G4endl; 295 <<G4endl;
270 G4cout <<"Initial target A=" <<AT 296 G4cout <<"Initial target A=" <<AT
271 <<", Z=" <<ZT 297 <<", Z=" <<ZT
272 <<", radius = " <<rT/fermi <<" fm" 298 <<", radius = " <<rT/fermi <<" fm"
273 <<G4endl; 299 <<G4endl;
274 G4cout <<"Projectile momentum and Energy/n 300 G4cout <<"Projectile momentum and Energy/nuc = " <<pP <<" ," <<E <<G4endl;
275 } 301 }
276 // 302 //
277 // 303 //
278 // The following variables are used to determi 304 // The following variables are used to determine the impact parameter in the
279 // near-field (i.e. taking into consideration 305 // near-field (i.e. taking into consideration the electrostatic repulsion).
280 // 306 //
281 G4double rm = ZP * ZT * elm_coupling / (E 307 G4double rm = ZP * ZT * elm_coupling / (E * AP);
282 G4double r = 0.0; 308 G4double r = 0.0;
283 G4double rsq = 0.0; 309 G4double rsq = 0.0;
284 // 310 //
285 // 311 //
286 // Initialise some of the variables which wll 312 // Initialise some of the variables which wll be used to calculate the chord-
287 // length for nucleons in the projectile and t 313 // length for nucleons in the projectile and target, and hence calculate the
288 // number of abraded nucleons and the excitati 314 // number of abraded nucleons and the excitation energy.
289 // 315 //
290 G4NuclearAbrasionGeometry *theAbrasionGeomet << 316 G4NuclearAbrasionGeometry *theAbrasionGeometry = NULL;
291 G4double CT = 0.0; 317 G4double CT = 0.0;
292 G4double F = 0.0; 318 G4double F = 0.0;
293 G4int Dabr = 0; 319 G4int Dabr = 0;
294 // 320 //
295 // 321 //
296 // The following loop is performed until the n 322 // The following loop is performed until the number of nucleons which are
297 // abraded by the process is >1, i.e. an inter 323 // abraded by the process is >1, i.e. an interaction MUST occur.
298 // 324 //
299 G4bool skipInteraction = false; // It will 325 G4bool skipInteraction = false; // It will be set true if the two nuclei fail to collide
300 const G4int maxNumberOfLoops = 1000; 326 const G4int maxNumberOfLoops = 1000;
301 G4int loopCounter = -1; 327 G4int loopCounter = -1;
302 while (Dabr == 0 && ++loopCounter < maxNumbe 328 while (Dabr == 0 && ++loopCounter < maxNumberOfLoops) /* Loop checking, 07.08.2015, A.Ribon */
303 { 329 {
>> 330 // Added by MHM 20050119 to fix leaking memory on second pass through this loop
>> 331 if (theAbrasionGeometry)
>> 332 {
>> 333 delete theAbrasionGeometry;
>> 334 theAbrasionGeometry = NULL;
>> 335 }
304 // 336 //
305 // 337 //
306 // Sample the impact parameter. For the momen 338 // Sample the impact parameter. For the moment, this class takes account of
307 // electrostatic effects on the impact paramet 339 // electrostatic effects on the impact parameter, but (like HZETRN AND NUCFRG2)
308 // does not make any correction for the effect 340 // does not make any correction for the effects of nuclear-nuclear repulsion.
309 // 341 //
310 G4double rPT = rP + rT; 342 G4double rPT = rP + rT;
311 G4double rPTsq = rPT * rPT; 343 G4double rPTsq = rPT * rPT;
312 // 344 //
313 // 345 //
314 // This is a "catch" to make sure we don't go 346 // This is a "catch" to make sure we don't go into an infinite loop because the
315 // energy is too low to overcome nuclear repul 347 // energy is too low to overcome nuclear repulsion. PRT 20091023. If the
316 // value of rm < fradius * rPT then we're unli 348 // value of rm < fradius * rPT then we're unlikely to sample a small enough
317 // impact parameter (energy of incident partic 349 // impact parameter (energy of incident particle is too low).
318 // 350 //
319 if (rm >= fradius * rPT) { 351 if (rm >= fradius * rPT) {
320 skipInteraction = true; 352 skipInteraction = true;
321 } 353 }
322 // 354 //
323 // 355 //
324 // Now sample impact parameter until the crite 356 // Now sample impact parameter until the criterion is met that projectile
325 // and target overlap, but repulsion is taken 357 // and target overlap, but repulsion is taken into consideration.
326 // 358 //
327 G4int evtcnt = 0; 359 G4int evtcnt = 0;
328 r = 1.1 * rPT; 360 r = 1.1 * rPT;
329 while (r > rPT && ++evtcnt < 1000) /* Loo 361 while (r > rPT && ++evtcnt < 1000) /* Loop checking, 07.08.2015, A.Ribon */
330 { 362 {
331 G4double bsq = rPTsq * G4UniformRand(); 363 G4double bsq = rPTsq * G4UniformRand();
332 r = (rm + std::sqrt(rm*rm + 4 364 r = (rm + std::sqrt(rm*rm + 4.0*bsq)) / 2.0;
333 } 365 }
334 // 366 //
335 // 367 //
336 // We've tried to sample this 1000 times, but 368 // We've tried to sample this 1000 times, but failed.
337 // 369 //
338 if (evtcnt >= 1000) { 370 if (evtcnt >= 1000) {
339 skipInteraction = true; 371 skipInteraction = true;
340 } 372 }
341 373
342 rsq = r * r; 374 rsq = r * r;
343 // 375 //
344 // 376 //
345 // Now determine the chord-length through the 377 // Now determine the chord-length through the target nucleus.
346 // 378 //
347 if (rT > rP) 379 if (rT > rP)
348 { 380 {
349 G4double x = (rPsq + rsq - rTsq) / 2.0 / 381 G4double x = (rPsq + rsq - rTsq) / 2.0 / r;
350 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - 382 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
351 else CT = 2.0 * std::sqrt(rTsq - 383 else CT = 2.0 * std::sqrt(rTsq - rsq);
352 } 384 }
353 else 385 else
354 { 386 {
355 G4double x = (rTsq + rsq - rPsq) / 2.0 / 387 G4double x = (rTsq + rsq - rPsq) / 2.0 / r;
356 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - 388 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
357 else CT = 2.0 * rT; 389 else CT = 2.0 * rT;
358 } 390 }
359 // 391 //
360 // 392 //
361 // Determine the number of abraded nucleons. 393 // Determine the number of abraded nucleons. Note that the mean number of
362 // abraded nucleons is used to sample the Pois 394 // abraded nucleons is used to sample the Poisson distribution. The Poisson
363 // distribution is sampled only ten times with 395 // distribution is sampled only ten times with the current impact parameter,
364 // and if it fails after this to find a case f 396 // and if it fails after this to find a case for which the number of abraded
365 // nucleons >1, the impact parameter is re-sam 397 // nucleons >1, the impact parameter is re-sampled.
366 // 398 //
367 delete theAbrasionGeometry; <<
368 theAbrasionGeometry = new G4NuclearAbrasio 399 theAbrasionGeometry = new G4NuclearAbrasionGeometry(AP,AT,r);
369 F = theAbrasionGeometry- 400 F = theAbrasionGeometry->F();
370 G4double lambda = 16.6*fermi / G4Pow:: 401 G4double lambda = 16.6*fermi / G4Pow::GetInstance()->powA(E/MeV,0.26);
371 G4double Mabr = F * AP * (1.0 - G4Ex 402 G4double Mabr = F * AP * (1.0 - G4Exp(-CT/lambda));
372 G4long n = 0; 403 G4long n = 0;
373 for (G4int i = 0; i<10; ++i) << 404 for (G4int i = 0; i<10; i++)
374 { 405 {
375 n = G4Poisson(Mabr); 406 n = G4Poisson(Mabr);
376 if (n > 0) 407 if (n > 0)
377 { 408 {
378 if (n>AP) Dabr = (G4int) AP; 409 if (n>AP) Dabr = (G4int) AP;
379 else Dabr = (G4int) n; 410 else Dabr = (G4int) n;
380 break; 411 break;
381 } 412 }
382 } 413 }
383 } // End of while loop 414 } // End of while loop
384 415
385 if ( loopCounter >= maxNumberOfLoops || skip 416 if ( loopCounter >= maxNumberOfLoops || skipInteraction ) {
386 // Assume nuclei do not collide and return 417 // Assume nuclei do not collide and return unchanged primary.
387 theParticleChange.SetStatusChange(isAlive) 418 theParticleChange.SetStatusChange(isAlive);
388 theParticleChange.SetEnergyChange(theTrack 419 theParticleChange.SetEnergyChange(theTrack.GetKineticEnergy());
389 theParticleChange.SetMomentumChange(theTra 420 theParticleChange.SetMomentumChange(theTrack.Get4Momentum().vect().unit());
390 if (verboseLevel >= 2) { 421 if (verboseLevel >= 2) {
391 G4cout <<"Particle energy too low to ove 422 G4cout <<"Particle energy too low to overcome repulsion." <<G4endl;
392 G4cout <<"Event rejected and original tr 423 G4cout <<"Event rejected and original track maintained" <<G4endl;
393 G4cout <<"############################## 424 G4cout <<"########################################"
394 <<"############################## 425 <<"########################################"
395 <<G4endl; 426 <<G4endl;
396 } 427 }
397 delete theAbrasionGeometry; <<
398 return &theParticleChange; 428 return &theParticleChange;
399 } 429 }
400 430
401 if (verboseLevel >= 2) 431 if (verboseLevel >= 2)
402 { 432 {
403 G4cout <<G4endl; 433 G4cout <<G4endl;
404 G4cout <<"Impact parameter = " <<r/ferm 434 G4cout <<"Impact parameter = " <<r/fermi <<" fm" <<G4endl;
405 G4cout <<"# Abraded nucleons = " <<Dabr < 435 G4cout <<"# Abraded nucleons = " <<Dabr <<G4endl;
406 } 436 }
407 // 437 //
408 // 438 //
409 // The number of abraded nucleons must be no g 439 // The number of abraded nucleons must be no greater than the number of
410 // nucleons in either the projectile or the ta 440 // nucleons in either the projectile or the target. If AP - Dabr < 2 or
411 // AT - Dabr < 2 then either we have only a nu 441 // AT - Dabr < 2 then either we have only a nucleon left behind in the
412 // projectile/target or we've tried to abrade 442 // projectile/target or we've tried to abrade too many nucleons - and Dabr
413 // should be limited. 443 // should be limited.
414 // 444 //
415 if (AP - (G4double) Dabr < 2.0) Dabr = (G4in 445 if (AP - (G4double) Dabr < 2.0) Dabr = (G4int) AP;
416 if (AT - (G4double) Dabr < 2.0) Dabr = (G4in 446 if (AT - (G4double) Dabr < 2.0) Dabr = (G4int) AT;
417 // 447 //
418 // 448 //
419 // Determine the abraded secondary nucleons fr 449 // Determine the abraded secondary nucleons from the projectile. *fragmentP
420 // is a pointer to the prefragment from the pr 450 // is a pointer to the prefragment from the projectile and nSecP is the number
421 // of nucleons in theParticleChange which have 451 // of nucleons in theParticleChange which have been abraded. The total energy
422 // from these is determined. 452 // from these is determined.
423 // 453 //
424 G4ThreeVector boost = pP.findBoostToCM(); 454 G4ThreeVector boost = pP.findBoostToCM();
425 G4Fragment *fragmentP = GetAbradedNucleons ( 455 G4Fragment *fragmentP = GetAbradedNucleons (Dabr, AP, ZP, rP);
426 G4int nSecP = (G4int)theParticleCh << 456 G4int nSecP = theParticleChange.GetNumberOfSecondaries();
427 G4int i = 0; 457 G4int i = 0;
428 for (i=0; i<nSecP; ++i) << 458 for (i=0; i<nSecP; i++)
429 { 459 {
430 TotalEPost += theParticleChange.GetSeconda 460 TotalEPost += theParticleChange.GetSecondary(i)->
431 GetParticle()->GetTotalEnergy(); 461 GetParticle()->GetTotalEnergy();
432 } 462 }
433 // 463 //
434 // 464 //
435 // Determine the number of spectators in the i 465 // Determine the number of spectators in the interaction region for the
436 // projectile. 466 // projectile.
437 // 467 //
438 G4int DspcP = (G4int) (AP*F) - Dabr; 468 G4int DspcP = (G4int) (AP*F) - Dabr;
439 if (DspcP <= 0) DspcP = 0; 469 if (DspcP <= 0) DspcP = 0;
440 else if (DspcP > AP-Dabr) DspcP = ((G4int) A 470 else if (DspcP > AP-Dabr) DspcP = ((G4int) AP) - Dabr;
441 // 471 //
442 // 472 //
443 // Determine excitation energy associated with 473 // Determine excitation energy associated with excess surface area of the
444 // projectile (EsP) and the excitation due to 474 // projectile (EsP) and the excitation due to scattering of nucleons which are
445 // retained within the projectile (ExP). Add 475 // retained within the projectile (ExP). Add the total energy from the excited
446 // nucleus to the total energy of the secondar 476 // nucleus to the total energy of the secondaries.
447 // 477 //
448 G4bool excitationAbsorbedByProjectile = fals 478 G4bool excitationAbsorbedByProjectile = false;
449 if (fragmentP != nullptr) << 479 if (fragmentP != NULL)
450 { 480 {
451 G4double EsP = theAbrasionGeometry->GetExc 481 G4double EsP = theAbrasionGeometry->GetExcitationEnergyOfProjectile();
452 G4double ExP = 0.0; 482 G4double ExP = 0.0;
453 if (Dabr < AT) 483 if (Dabr < AT)
454 excitationAbsorbedByProjectile = G4Unifo 484 excitationAbsorbedByProjectile = G4UniformRand() < 0.5;
455 if (excitationAbsorbedByProjectile) 485 if (excitationAbsorbedByProjectile)
456 ExP = GetNucleonInducedExcitation(rP, rT 486 ExP = GetNucleonInducedExcitation(rP, rT, r);
457 G4double xP = EsP + ExP; 487 G4double xP = EsP + ExP;
458 if (xP > B*(AP-Dabr)) xP = B*(AP-Dabr); 488 if (xP > B*(AP-Dabr)) xP = B*(AP-Dabr);
459 G4LorentzVector lorentzVector = fragmentP- 489 G4LorentzVector lorentzVector = fragmentP->GetMomentum();
460 lorentzVector.setE(lorentzVector.e()+xP); 490 lorentzVector.setE(lorentzVector.e()+xP);
461 fragmentP->SetMomentum(lorentzVector); 491 fragmentP->SetMomentum(lorentzVector);
462 TotalEPost += lorentzVector.e(); 492 TotalEPost += lorentzVector.e();
463 } 493 }
464 G4double EMassP = TotalEPost; 494 G4double EMassP = TotalEPost;
465 // 495 //
466 // 496 //
467 // Determine the abraded secondary nucleons fr 497 // Determine the abraded secondary nucleons from the target. Note that it's
468 // assumed that the same number of nucleons ar 498 // assumed that the same number of nucleons are abraded from the target as for
469 // the projectile, and obviously no boost is a 499 // the projectile, and obviously no boost is applied to the products. *fragmentT
470 // is a pointer to the prefragment from the ta 500 // is a pointer to the prefragment from the target and nSec is the total number
471 // of nucleons in theParticleChange which have 501 // of nucleons in theParticleChange which have been abraded. The total energy
472 // from these is determined. 502 // from these is determined.
473 // 503 //
474 G4Fragment *fragmentT = GetAbradedNucleons ( 504 G4Fragment *fragmentT = GetAbradedNucleons (Dabr, AT, ZT, rT);
475 G4int nSec = (G4int)theParticleChange.GetNum << 505 G4int nSec = theParticleChange.GetNumberOfSecondaries();
476 for (i=nSecP; i<nSec; ++i) << 506 for (i=nSecP; i<nSec; i++)
477 { 507 {
478 TotalEPost += theParticleChange.GetSeconda 508 TotalEPost += theParticleChange.GetSecondary(i)->
479 GetParticle()->GetTotalEnergy(); 509 GetParticle()->GetTotalEnergy();
480 } 510 }
481 // 511 //
482 // 512 //
483 // Determine the number of spectators in the i 513 // Determine the number of spectators in the interaction region for the
484 // target. 514 // target.
485 // 515 //
486 G4int DspcT = (G4int) (AT*F) - Dabr; 516 G4int DspcT = (G4int) (AT*F) - Dabr;
487 if (DspcT <= 0) DspcT = 0; 517 if (DspcT <= 0) DspcT = 0;
488 else if (DspcT > AP-Dabr) DspcT = ((G4int) A 518 else if (DspcT > AP-Dabr) DspcT = ((G4int) AT) - Dabr;
489 // 519 //
490 // 520 //
491 // Determine excitation energy associated with 521 // Determine excitation energy associated with excess surface area of the
492 // target (EsT) and the excitation due to scat 522 // target (EsT) and the excitation due to scattering of nucleons which are
493 // retained within the target (ExT). Add the 523 // retained within the target (ExT). Add the total energy from the excited
494 // nucleus to the total energy of the secondar 524 // nucleus to the total energy of the secondaries.
495 // 525 //
496 if (fragmentT != nullptr) << 526 if (fragmentT != NULL)
497 { 527 {
498 G4double EsT = theAbrasionGeometry->GetExc 528 G4double EsT = theAbrasionGeometry->GetExcitationEnergyOfTarget();
499 G4double ExT = 0.0; 529 G4double ExT = 0.0;
500 if (!excitationAbsorbedByProjectile) 530 if (!excitationAbsorbedByProjectile)
501 ExT = GetNucleonInducedExcitation(rT, rP 531 ExT = GetNucleonInducedExcitation(rT, rP, r);
502 G4double xT = EsT + ExT; 532 G4double xT = EsT + ExT;
503 if (xT > B*(AT-Dabr)) xT = B*(AT-Dabr); 533 if (xT > B*(AT-Dabr)) xT = B*(AT-Dabr);
504 G4LorentzVector lorentzVector = fragmentT- 534 G4LorentzVector lorentzVector = fragmentT->GetMomentum();
505 lorentzVector.setE(lorentzVector.e()+xT); 535 lorentzVector.setE(lorentzVector.e()+xT);
506 fragmentT->SetMomentum(lorentzVector); 536 fragmentT->SetMomentum(lorentzVector);
507 TotalEPost += lorentzVector.e(); 537 TotalEPost += lorentzVector.e();
508 } 538 }
509 // 539 //
510 // 540 //
511 // Now determine the difference between the pr 541 // Now determine the difference between the pre and post interaction
512 // energy - this will be used to determine the 542 // energy - this will be used to determine the Lorentz boost if conservation
513 // of energy is to be imposed/attempted. 543 // of energy is to be imposed/attempted.
514 // 544 //
515 G4double deltaE = TotalEPre - TotalEPost; 545 G4double deltaE = TotalEPre - TotalEPost;
516 if (deltaE > 0.0 && conserveEnergy) 546 if (deltaE > 0.0 && conserveEnergy)
517 { 547 {
518 G4double beta = std::sqrt(1.0 - EMassP*EMa 548 G4double beta = std::sqrt(1.0 - EMassP*EMassP/G4Pow::GetInstance()->powN(deltaE+EMassP,2));
519 boost = boost / boost.mag() * beta; 549 boost = boost / boost.mag() * beta;
520 } 550 }
521 // 551 //
522 // 552 //
523 // Now boost the secondaries from the projecti 553 // Now boost the secondaries from the projectile.
524 // 554 //
525 G4ThreeVector pBalance = pP.vect(); 555 G4ThreeVector pBalance = pP.vect();
526 for (i=0; i<nSecP; ++i) << 556 for (i=0; i<nSecP; i++)
527 { 557 {
528 G4DynamicParticle *dynamicP = theParticleC 558 G4DynamicParticle *dynamicP = theParticleChange.GetSecondary(i)->
529 GetParticle(); 559 GetParticle();
530 G4LorentzVector lorentzVector = dynamicP-> 560 G4LorentzVector lorentzVector = dynamicP->Get4Momentum();
531 lorentzVector.boost(-boost); 561 lorentzVector.boost(-boost);
532 dynamicP->Set4Momentum(lorentzVector); 562 dynamicP->Set4Momentum(lorentzVector);
533 pBalance -= lorentzVector.vect(); 563 pBalance -= lorentzVector.vect();
534 } 564 }
535 // 565 //
536 // 566 //
537 // Set the boost for the projectile prefragmen 567 // Set the boost for the projectile prefragment. This is now based on the
538 // conservation of momentum. However, if the 568 // conservation of momentum. However, if the user selected momentum of the
539 // prefragment is not to be conserved this sim 569 // prefragment is not to be conserved this simply boosted to the velocity of the
540 // original projectile times the ratio of the 570 // original projectile times the ratio of the unexcited to the excited mass
541 // of the prefragment (the excitation increase 571 // of the prefragment (the excitation increases the effective mass of the
542 // prefragment, and therefore modifying the bo 572 // prefragment, and therefore modifying the boost is an attempt to prevent
543 // the momentum of the prefragment being exces 573 // the momentum of the prefragment being excessive).
544 // 574 //
545 if (fragmentP != nullptr) << 575 if (fragmentP != NULL)
546 { 576 {
547 G4LorentzVector lorentzVector = fragmentP- 577 G4LorentzVector lorentzVector = fragmentP->GetMomentum();
548 G4double fragmentM = lorentzVec 578 G4double fragmentM = lorentzVector.m();
549 if (conserveMomentum) 579 if (conserveMomentum)
550 fragmentP->SetMomentum 580 fragmentP->SetMomentum
551 (G4LorentzVector(pBalance,std::sqrt(pB 581 (G4LorentzVector(pBalance,std::sqrt(pBalance.mag2()+fragmentM*fragmentM+1.0*eV*eV)));
552 else 582 else
553 { 583 {
554 G4double fragmentGroundStateM = fragment 584 G4double fragmentGroundStateM = fragmentP->GetGroundStateMass();
555 fragmentP->SetMomentum(lorentzVector.boo 585 fragmentP->SetMomentum(lorentzVector.boost(-boost * fragmentGroundStateM/fragmentM));
556 } 586 }
557 } 587 }
558 // 588 //
559 // 589 //
560 // Output information to user if verbose infor 590 // Output information to user if verbose information requested.
561 // 591 //
562 if (verboseLevel >= 2) 592 if (verboseLevel >= 2)
563 { 593 {
564 G4cout <<G4endl; 594 G4cout <<G4endl;
565 G4cout <<"-------------------------------- 595 G4cout <<"-----------------------------------" <<G4endl;
566 G4cout <<"Secondary nucleons from projecti 596 G4cout <<"Secondary nucleons from projectile:" <<G4endl;
567 G4cout <<"-------------------------------- 597 G4cout <<"-----------------------------------" <<G4endl;
568 G4cout.precision(7); 598 G4cout.precision(7);
569 for (i=0; i<nSecP; ++i) << 599 for (i=0; i<nSecP; i++)
570 { 600 {
571 G4cout <<"Particle # " <<i <<G4endl; 601 G4cout <<"Particle # " <<i <<G4endl;
572 theParticleChange.GetSecondary(i)->GetPa 602 theParticleChange.GetSecondary(i)->GetParticle()->DumpInfo();
573 G4DynamicParticle *dyn = theParticleChan 603 G4DynamicParticle *dyn = theParticleChange.GetSecondary(i)->GetParticle();
574 G4cout <<"New nucleon (P) " <<dyn->GetDe 604 G4cout <<"New nucleon (P) " <<dyn->GetDefinition()->GetParticleName()
575 <<" : " <<dyn->Get4M 605 <<" : " <<dyn->Get4Momentum()
576 <<G4endl; 606 <<G4endl;
577 } 607 }
578 G4cout <<"---------------------------" <<G 608 G4cout <<"---------------------------" <<G4endl;
579 G4cout <<"The projectile prefragment:" <<G 609 G4cout <<"The projectile prefragment:" <<G4endl;
580 G4cout <<"---------------------------" <<G 610 G4cout <<"---------------------------" <<G4endl;
581 if (fragmentP != nullptr) << 611 if (fragmentP != NULL)
582 G4cout <<*fragmentP <<G4endl; 612 G4cout <<*fragmentP <<G4endl;
583 else 613 else
584 G4cout <<"(No residual prefragment)" <<G 614 G4cout <<"(No residual prefragment)" <<G4endl;
585 G4cout <<G4endl; 615 G4cout <<G4endl;
586 G4cout <<"-------------------------------" 616 G4cout <<"-------------------------------" <<G4endl;
587 G4cout <<"Secondary nucleons from target:" 617 G4cout <<"Secondary nucleons from target:" <<G4endl;
588 G4cout <<"-------------------------------" 618 G4cout <<"-------------------------------" <<G4endl;
589 G4cout.precision(7); 619 G4cout.precision(7);
590 for (i=nSecP; i<nSec; ++i) << 620 for (i=nSecP; i<nSec; i++)
591 { 621 {
592 G4cout <<"Particle # " <<i <<G4endl; 622 G4cout <<"Particle # " <<i <<G4endl;
593 theParticleChange.GetSecondary(i)->GetPa 623 theParticleChange.GetSecondary(i)->GetParticle()->DumpInfo();
594 G4DynamicParticle *dyn = theParticleChan 624 G4DynamicParticle *dyn = theParticleChange.GetSecondary(i)->GetParticle();
595 G4cout <<"New nucleon (T) " <<dyn->GetDe 625 G4cout <<"New nucleon (T) " <<dyn->GetDefinition()->GetParticleName()
596 <<" : " <<dyn->Get4M 626 <<" : " <<dyn->Get4Momentum()
597 <<G4endl; 627 <<G4endl;
598 } 628 }
599 G4cout <<"-----------------------" <<G4end 629 G4cout <<"-----------------------" <<G4endl;
600 G4cout <<"The target prefragment:" <<G4end 630 G4cout <<"The target prefragment:" <<G4endl;
601 G4cout <<"-----------------------" <<G4end 631 G4cout <<"-----------------------" <<G4endl;
602 if (fragmentT != nullptr) << 632 if (fragmentT != NULL)
603 G4cout <<*fragmentT <<G4endl; 633 G4cout <<*fragmentT <<G4endl;
604 else 634 else
605 G4cout <<"(No residual prefragment)" <<G 635 G4cout <<"(No residual prefragment)" <<G4endl;
606 } 636 }
607 // 637 //
608 // 638 //
609 // Now we can decay the nuclear fragments if p 639 // Now we can decay the nuclear fragments if present. The secondaries are
610 // collected and boosted as well. This is per 640 // collected and boosted as well. This is performed first for the projectile...
611 // 641 //
612 if (fragmentP !=nullptr) << 642 if (fragmentP !=NULL)
613 { 643 {
614 G4ReactionProductVector *products = nullpt << 644 G4ReactionProductVector *products = NULL;
615 // if (fragmentP->GetZ_asInt() != fragm << 645 if (fragmentP->GetZ_asInt() != fragmentP->GetA_asInt())
616 products = theExcitationHandler->BreakItUp << 646 products = theExcitationHandler->BreakItUp(*fragmentP);
617 // else << 647 else
618 // products = theExcitationHandlerx->Bre << 648 products = theExcitationHandlerx->BreakItUp(*fragmentP);
619 delete fragmentP; 649 delete fragmentP;
620 fragmentP = nullptr; << 650 fragmentP = NULL;
621 651
622 G4ReactionProductVector::iterator iter; 652 G4ReactionProductVector::iterator iter;
623 for (iter = products->begin(); iter != pro 653 for (iter = products->begin(); iter != products->end(); ++iter)
624 { 654 {
625 G4DynamicParticle *secondary = 655 G4DynamicParticle *secondary =
626 new G4DynamicParticle((*iter)->GetDefi 656 new G4DynamicParticle((*iter)->GetDefinition(),
627 (*iter)->GetTotalEnergy(), (*iter)->Ge 657 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
628 theParticleChange.AddSecondary (secondar << 658 theParticleChange.AddSecondary (secondary); // Added MHM 20050118
629 G4String particleName = (*iter)->GetDefi 659 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
630 delete (*iter); // get rid of leftover p << 660 delete (*iter); // get rid of leftover particle def! // Added MHM 20050118
631 if (verboseLevel >= 2 && particleName.fi 661 if (verboseLevel >= 2 && particleName.find("[",0) < particleName.size())
632 { 662 {
633 G4cout <<"------------------------" << 663 G4cout <<"------------------------" <<G4endl;
634 G4cout <<"The projectile fragment:" << 664 G4cout <<"The projectile fragment:" <<G4endl;
635 G4cout <<"------------------------" << 665 G4cout <<"------------------------" <<G4endl;
636 G4cout <<" fragmentP = " <<particleNam 666 G4cout <<" fragmentP = " <<particleName
637 <<" Energy = " <<secondary-> 667 <<" Energy = " <<secondary->GetKineticEnergy()
638 <<G4endl; 668 <<G4endl;
639 } 669 }
640 } 670 }
641 delete products; << 671 delete products; // Added MHM 20050118
642 } 672 }
643 // 673 //
644 // 674 //
645 // Now decay the target nucleus - no boost is 675 // Now decay the target nucleus - no boost is applied since in this
646 // approximation it is assumed that there is n 676 // approximation it is assumed that there is negligible momentum transfer from
647 // the projectile. 677 // the projectile.
648 // 678 //
649 if (fragmentT != nullptr) << 679 if (fragmentT != NULL)
650 { 680 {
651 G4ReactionProductVector *products = nullpt << 681 G4ReactionProductVector *products = NULL;
652 // if (fragmentT->GetZ_asInt() != fragm << 682 if (fragmentT->GetZ_asInt() != fragmentT->GetA_asInt())
653 products = theExcitationHandler->BreakIt 683 products = theExcitationHandler->BreakItUp(*fragmentT);
654 // else << 684 else
655 // products = theExcitationHandlerx->Bre << 685 products = theExcitationHandlerx->BreakItUp(*fragmentT);
656 // delete fragmentT; << 686 delete fragmentT;
657 fragmentT = nullptr; << 687 fragmentT = NULL;
658 688
659 G4ReactionProductVector::iterator iter; 689 G4ReactionProductVector::iterator iter;
660 for (iter = products->begin(); iter != pro 690 for (iter = products->begin(); iter != products->end(); ++iter)
661 { 691 {
662 G4DynamicParticle *secondary = 692 G4DynamicParticle *secondary =
663 new G4DynamicParticle((*iter)->GetDefi 693 new G4DynamicParticle((*iter)->GetDefinition(),
664 (*iter)->GetTotalEnergy(), (*iter)->Ge 694 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
665 theParticleChange.AddSecondary (secondar << 695 theParticleChange.AddSecondary (secondary);
666 G4String particleName = (*iter)->GetDefi 696 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
667 delete (*iter); // get rid of leftover p << 697 delete (*iter); // get rid of leftover particle def! // Added MHM 20050118
668 if (verboseLevel >= 2 && particleName.fi 698 if (verboseLevel >= 2 && particleName.find("[",0) < particleName.size())
669 { 699 {
670 G4cout <<"--------------------" <<G4en 700 G4cout <<"--------------------" <<G4endl;
671 G4cout <<"The target fragment:" <<G4en 701 G4cout <<"The target fragment:" <<G4endl;
672 G4cout <<"--------------------" <<G4en 702 G4cout <<"--------------------" <<G4endl;
673 G4cout <<" fragmentT = " <<particleNam 703 G4cout <<" fragmentT = " <<particleName
674 <<" Energy = " <<secondary-> 704 <<" Energy = " <<secondary->GetKineticEnergy()
675 <<G4endl; 705 <<G4endl;
676 } 706 }
677 } 707 }
678 delete products; << 708 delete products; // Added MHM 20050118
679 } 709 }
680 710
681 if (verboseLevel >= 2) 711 if (verboseLevel >= 2)
682 G4cout <<"############################### 712 G4cout <<"########################################"
683 <<"############################### 713 <<"########################################"
684 <<G4endl; 714 <<G4endl;
685 715
686 delete theAbrasionGeometry; 716 delete theAbrasionGeometry;
>> 717
687 return &theParticleChange; 718 return &theParticleChange;
688 } 719 }
689 ////////////////////////////////////////////// 720 ////////////////////////////////////////////////////////////////////////////////
690 // 721 //
691 G4Fragment *G4WilsonAbrasionModel::GetAbradedN 722 G4Fragment *G4WilsonAbrasionModel::GetAbradedNucleons (G4int Dabr, G4double A,
692 G4double Z, G4double r) 723 G4double Z, G4double r)
693 { 724 {
694 // 725 //
695 // 726 //
696 // Initialise variables. tau is the Fermi rad 727 // Initialise variables. tau is the Fermi radius of the nucleus. The variables
697 // p..., C... and gamma are used to help sampl 728 // p..., C... and gamma are used to help sample the secondary nucleon
698 // spectrum. 729 // spectrum.
699 // 730 //
700 731
701 G4double pK = hbarc * G4Pow::GetInstance() 732 G4double pK = hbarc * G4Pow::GetInstance()->A13(9.0 * pi / 4.0 * A) / (1.29 * r);
702 if (A <= 24.0) pK *= -0.229*G4Pow::GetInstan 733 if (A <= 24.0) pK *= -0.229*G4Pow::GetInstance()->A13(A) + 1.62;
703 G4double pKsq = pK * pK; 734 G4double pKsq = pK * pK;
704 G4double p1sq = 2.0/5.0 * pKsq; 735 G4double p1sq = 2.0/5.0 * pKsq;
705 G4double p2sq = 6.0/5.0 * pKsq; 736 G4double p2sq = 6.0/5.0 * pKsq;
706 G4double p3sq = 500.0 * 500.0; 737 G4double p3sq = 500.0 * 500.0;
707 G4double C1 = 1.0; 738 G4double C1 = 1.0;
708 G4double C2 = 0.03; 739 G4double C2 = 0.03;
709 G4double C3 = 0.0002; 740 G4double C3 = 0.0002;
710 G4double gamma = 90.0 * MeV; 741 G4double gamma = 90.0 * MeV;
711 G4double maxn = C1 + C2 + C3; 742 G4double maxn = C1 + C2 + C3;
712 // 743 //
713 // 744 //
714 // initialise the number of secondary nucleons 745 // initialise the number of secondary nucleons abraded to zero, and initially set
715 // the type of nucleon abraded to proton ... j 746 // the type of nucleon abraded to proton ... just for now.
716 // 747 //
717 G4double Aabr = 0.0; 748 G4double Aabr = 0.0;
718 G4double Zabr = 0.0; 749 G4double Zabr = 0.0;
719 G4ParticleDefinition *typeNucleon = G4Proton 750 G4ParticleDefinition *typeNucleon = G4Proton::ProtonDefinition();
720 G4DynamicParticle *dynamicNucleon = nullptr; << 751 G4DynamicParticle *dynamicNucleon = NULL;
721 G4ParticleMomentum pabr(0.0, 0.0, 0.0); 752 G4ParticleMomentum pabr(0.0, 0.0, 0.0);
722 // 753 //
723 // 754 //
724 // Now go through each abraded nucleon and sam 755 // Now go through each abraded nucleon and sample type, spectrum and angle.
725 // 756 //
726 G4bool isForLoopExitAnticipated = false; 757 G4bool isForLoopExitAnticipated = false;
727 for (G4int i=0; i<Dabr; ++i) << 758 for (G4int i=0; i<Dabr; i++)
728 { 759 {
729 // 760 //
730 // 761 //
731 // Sample the nucleon momentum distribution by 762 // Sample the nucleon momentum distribution by simple rejection techniques. We
732 // reject values of p == 0.0 since this causes 763 // reject values of p == 0.0 since this causes bad behaviour in the sinh term.
733 // 764 //
734 G4double p = 0.0; 765 G4double p = 0.0;
735 G4bool found = false; 766 G4bool found = false;
736 const G4int maxNumberOfLoops = 100000; 767 const G4int maxNumberOfLoops = 100000;
737 G4int loopCounter = -1; 768 G4int loopCounter = -1;
738 while (!found && ++loopCounter < maxNumber 769 while (!found && ++loopCounter < maxNumberOfLoops) /* Loop checking, 07.08.2015, A.Ribon */
739 { 770 {
740 while (p <= 0.0) p = npK * pK * G4Unifor 771 while (p <= 0.0) p = npK * pK * G4UniformRand(); /* Loop checking, 07.08.2015, A.Ribon */
741 G4double psq = p * p; 772 G4double psq = p * p;
742 found = maxn * G4UniformRand() < C1*G4Ex 773 found = maxn * G4UniformRand() < C1*G4Exp(-psq/p1sq/2.0) +
743 C2*G4Exp(-psq/p2sq/2.0) + C3*G4Exp(-ps 774 C2*G4Exp(-psq/p2sq/2.0) + C3*G4Exp(-psq/p3sq/2.0) + p/gamma/(0.5*(G4Exp(p/gamma)-G4Exp(-p/gamma)));
744 } 775 }
745 if ( loopCounter >= maxNumberOfLoops ) 776 if ( loopCounter >= maxNumberOfLoops )
746 { 777 {
747 isForLoopExitAnticipated = true; 778 isForLoopExitAnticipated = true;
748 break; 779 break;
749 } 780 }
750 // 781 //
751 // 782 //
752 // Determine the type of particle abraded. Ca 783 // Determine the type of particle abraded. Can only be proton or neutron,
753 // and the probability is determine to be prop 784 // and the probability is determine to be proportional to the ratio as found
754 // in the nucleus at each stage. 785 // in the nucleus at each stage.
755 // 786 //
756 G4double prob = (Z-Zabr)/(A-Aabr); 787 G4double prob = (Z-Zabr)/(A-Aabr);
757 if (G4UniformRand()<prob) 788 if (G4UniformRand()<prob)
758 { 789 {
759 Zabr++; 790 Zabr++;
760 typeNucleon = G4Proton::ProtonDefinition 791 typeNucleon = G4Proton::ProtonDefinition();
761 } 792 }
762 else 793 else
763 typeNucleon = G4Neutron::NeutronDefiniti 794 typeNucleon = G4Neutron::NeutronDefinition();
764 Aabr++; 795 Aabr++;
765 // 796 //
766 // 797 //
767 // The angular distribution of the secondary n 798 // The angular distribution of the secondary nucleons is approximated to an
768 // isotropic distribution in the rest frame of 799 // isotropic distribution in the rest frame of the nucleus (this will be Lorentz
769 // boosted later. 800 // boosted later.
770 // 801 //
771 G4double costheta = 2.*G4UniformRand()-1.0 802 G4double costheta = 2.*G4UniformRand()-1.0;
772 G4double sintheta = std::sqrt((1.0 - costh 803 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
773 G4double phi = 2.0*pi*G4UniformRand() 804 G4double phi = 2.0*pi*G4UniformRand()*rad;
774 G4ThreeVector direction(sintheta*std::cos( 805 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
775 G4double nucleonMass = typeNucleon->GetPDG 806 G4double nucleonMass = typeNucleon->GetPDGMass();
776 G4double E = std::sqrt(p*p + nuc 807 G4double E = std::sqrt(p*p + nucleonMass*nucleonMass)-nucleonMass;
777 dynamicNucleon = new G4DynamicParticle(typ 808 dynamicNucleon = new G4DynamicParticle(typeNucleon,direction,E);
778 theParticleChange.AddSecondary (dynamicNuc << 809 theParticleChange.AddSecondary (dynamicNucleon);
779 pabr += p*direction; 810 pabr += p*direction;
780 } 811 }
781 // 812 //
782 // 813 //
783 // Next determine the details of the nuclear p 814 // Next determine the details of the nuclear prefragment .. that is if there
784 // is one or more protons in the residue. (No 815 // is one or more protons in the residue. (Note that the 1 eV in the total
785 // energy is a safety factor to avoid any poss 816 // energy is a safety factor to avoid any possibility of negative rest mass
786 // energy.) 817 // energy.)
787 // 818 //
788 G4Fragment *fragment = nullptr; << 819 G4Fragment *fragment = NULL;
789 if ( ! isForLoopExitAnticipated && Z-Zabr>=1 820 if ( ! isForLoopExitAnticipated && Z-Zabr>=1.0 )
790 { 821 {
791 G4double ionMass = G4ParticleTable::GetPar 822 G4double ionMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
792 GetIonMass(G4lrint(Z-Za 823 GetIonMass(G4lrint(Z-Zabr),G4lrint(A-Aabr));
793 G4double E = std::sqrt(pabr.mag2() + 824 G4double E = std::sqrt(pabr.mag2() + ionMass*ionMass);
794 G4LorentzVector lorentzVector = G4LorentzV 825 G4LorentzVector lorentzVector = G4LorentzVector(-pabr, E + 1.0*eV);
795 fragment = 826 fragment =
796 new G4Fragment((G4int) (A-Aabr), (G4int) 827 new G4Fragment((G4int) (A-Aabr), (G4int) (Z-Zabr), lorentzVector);
797 } 828 }
798 829
799 return fragment; 830 return fragment;
800 } 831 }
801 ////////////////////////////////////////////// 832 ////////////////////////////////////////////////////////////////////////////////
802 // 833 //
803 G4double G4WilsonAbrasionModel::GetNucleonIndu 834 G4double G4WilsonAbrasionModel::GetNucleonInducedExcitation
804 (G4double rP, G4double rT, G4double r) 835 (G4double rP, G4double rT, G4double r)
805 { 836 {
806 // 837 //
807 // 838 //
808 // Initialise variables. 839 // Initialise variables.
809 // 840 //
810 G4double Cl = 0.0; 841 G4double Cl = 0.0;
811 G4double rPsq = rP * rP; 842 G4double rPsq = rP * rP;
812 G4double rTsq = rT * rT; 843 G4double rTsq = rT * rT;
813 G4double rsq = r * r; 844 G4double rsq = r * r;
814 // 845 //
815 // 846 //
816 // Depending upon the impact parameter, a diff 847 // Depending upon the impact parameter, a different form of the chord length is
817 // is used. 848 // is used.
818 // 849 //
819 if (r > rT) Cl = 2.0*std::sqrt(rPsq + 2.0*r* 850 if (r > rT) Cl = 2.0*std::sqrt(rPsq + 2.0*r*rT - rsq - rTsq);
820 else Cl = 2.0*rP; 851 else Cl = 2.0*rP;
821 // 852 //
822 // 853 //
823 // The next lines have been changed to include 854 // The next lines have been changed to include a "catch" to make sure if the
824 // projectile and target are too close, Ct is 855 // projectile and target are too close, Ct is set to twice rP or twice rT.
825 // Otherwise the standard Wilson algorithm sho 856 // Otherwise the standard Wilson algorithm should work fine.
826 // PRT 20091023. 857 // PRT 20091023.
827 // 858 //
828 G4double Ct = 0.0; 859 G4double Ct = 0.0;
829 if (rT > rP && rsq < rTsq - rPsq) Ct = 860 if (rT > rP && rsq < rTsq - rPsq) Ct = 2.0 * rP;
830 else if (rP > rT && rsq < rPsq - rTsq) Ct = 861 else if (rP > rT && rsq < rPsq - rTsq) Ct = 2.0 * rT;
831 else { 862 else {
832 G4double bP = (rPsq+rsq-rTsq)/2.0/r; 863 G4double bP = (rPsq+rsq-rTsq)/2.0/r;
833 G4double x = rPsq - bP*bP; 864 G4double x = rPsq - bP*bP;
834 if (x < 0.0) { 865 if (x < 0.0) {
835 G4cerr <<"############################## 866 G4cerr <<"########################################"
836 <<"############################## 867 <<"########################################"
837 <<G4endl; 868 <<G4endl;
838 G4cerr <<"ERROR IN G4WilsonAbrasionModel 869 G4cerr <<"ERROR IN G4WilsonAbrasionModel::GetNucleonInducedExcitation"
839 <<G4endl; 870 <<G4endl;
840 G4cerr <<"rPsq - bP*bP < 0.0 and cannot 871 G4cerr <<"rPsq - bP*bP < 0.0 and cannot be square-rooted" <<G4endl;
841 G4cerr <<"Set to zero instead" <<G4endl; 872 G4cerr <<"Set to zero instead" <<G4endl;
842 G4cerr <<"############################## 873 G4cerr <<"########################################"
843 <<"############################## 874 <<"########################################"
844 <<G4endl; 875 <<G4endl;
845 } 876 }
846 Ct = 2.0*std::sqrt(x); 877 Ct = 2.0*std::sqrt(x);
847 } 878 }
848 879
849 G4double Ex = 13.0 * Cl / fermi; 880 G4double Ex = 13.0 * Cl / fermi;
850 if (Ct > 1.5*fermi) 881 if (Ct > 1.5*fermi)
851 Ex += 13.0 * Cl / fermi /3.0 * (Ct/fermi - 882 Ex += 13.0 * Cl / fermi /3.0 * (Ct/fermi - 1.5);
852 883
853 return Ex; 884 return Ex;
854 } 885 }
855 ////////////////////////////////////////////// 886 ////////////////////////////////////////////////////////////////////////////////
856 // 887 //
857 void G4WilsonAbrasionModel::SetUseAblation (G4 888 void G4WilsonAbrasionModel::SetUseAblation (G4bool useAblation1)
858 { 889 {
859 if (useAblation != useAblation1) 890 if (useAblation != useAblation1)
860 { 891 {
861 useAblation = useAblation1; 892 useAblation = useAblation1;
>> 893 delete theExcitationHandler;
>> 894 delete theExcitationHandlerx;
>> 895 theExcitationHandler = new G4ExcitationHandler;
>> 896 theExcitationHandlerx = new G4ExcitationHandler;
862 if (useAblation) 897 if (useAblation)
863 { 898 {
864 theAblation = new G4WilsonAblationModel; 899 theAblation = new G4WilsonAblationModel;
865 theAblation->SetVerboseLevel(verboseLeve 900 theAblation->SetVerboseLevel(verboseLevel);
866 theExcitationHandler->SetEvaporation(the 901 theExcitationHandler->SetEvaporation(theAblation);
>> 902 theExcitationHandlerx->SetEvaporation(theAblation);
867 } 903 }
868 else 904 else
869 { 905 {
870 delete theExcitationHandler; << 906 theAblation = NULL;
871 theAblation = nullp << 907 G4Evaporation * theEvaporation = new G4Evaporation;
872 theExcitationHandler = new G4Excitation << 908 G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp;
>> 909 G4StatMF * theMF = new G4StatMF;
>> 910 theExcitationHandler->SetEvaporation(theEvaporation);
>> 911 theExcitationHandler->SetFermiModel(theFermiBreakUp);
>> 912 theExcitationHandler->SetMultiFragmentation(theMF);
>> 913 theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6);
>> 914 theExcitationHandler->SetMinEForMultiFrag(5.0*MeV);
>> 915
>> 916 theEvaporation = new G4Evaporation;
>> 917 theFermiBreakUp = new G4FermiBreakUp;
>> 918 theExcitationHandlerx->SetEvaporation(theEvaporation);
>> 919 theExcitationHandlerx->SetFermiModel(theFermiBreakUp);
>> 920 theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6);
873 } 921 }
874 } 922 }
875 return; 923 return;
876 } 924 }
877 ////////////////////////////////////////////// 925 ////////////////////////////////////////////////////////////////////////////////
878 // 926 //
879 void G4WilsonAbrasionModel::PrintWelcomeMessag 927 void G4WilsonAbrasionModel::PrintWelcomeMessage ()
880 { 928 {
881 G4cout <<G4endl; 929 G4cout <<G4endl;
882 G4cout <<" ********************************* 930 G4cout <<" *****************************************************************"
883 <<G4endl; 931 <<G4endl;
884 G4cout <<" Nuclear abrasion model for nuclea 932 G4cout <<" Nuclear abrasion model for nuclear-nuclear interactions activated"
885 <<G4endl; 933 <<G4endl;
886 G4cout <<" (Written by QinetiQ Ltd for the E 934 G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)"
887 <<G4endl; 935 <<G4endl;
888 G4cout <<" ********************************* 936 G4cout <<" *****************************************************************"
889 <<G4endl; 937 <<G4endl;
890 G4cout << G4endl; 938 G4cout << G4endl;
891 939
892 return; 940 return;
893 } 941 }
894 ////////////////////////////////////////////// 942 ////////////////////////////////////////////////////////////////////////////////
895 // 943 //
896 944