| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
27 //
28 // 230309 Skyrme-QMD parameters added by Y-H. Sato and A. Haga
29 // 230309 Total energy evaluated by Lorentz covariant version by Y-H. Sato and A. Haga
30
31 #include <numeric>
32
33 #include "G4LightIonQMDNucleus.hh"
34 #include "G4Pow.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4Proton.hh"
37 #include "G4Neutron.hh"
38 #include "G4NucleiProperties.hh"
39 #include "G4HadronicException.hh"
40
41 #include "G4LightIonQMDParameters.hh" // 20230309
42 #include "G4PhysicalConstants.hh" // 20230309
43 #include <cmath> // 20230309
44 #include <CLHEP/Random/Stat.h> // 20230309
45
46 G4LightIonQMDNucleus::G4LightIonQMDNucleus()
47 {
48 G4LightIonQMDParameters* parameters = G4LightIonQMDParameters::GetInstance();
49 hbc = parameters->Get_hbc();
50
51 jj = 0; // will be calcualted in CalEnergyAndAngularMomentumInCM;
52 potentialEnergy = 0.0; // will be set through set method
53 excitationEnergy = 0.0;
54
55 // Following Parameters are added (20230309)
56 wl = parameters->Get_wl();
57 cl = parameters->Get_cl();
58 rho0 = parameters->Get_rho0();
59 gamm = parameters->Get_gamm();
60 eta = parameters->Get_eta(); // Skyrme-QMD
61 kappas = parameters->Get_kappas(); // Skyrme-QMD
62
63 cpw = parameters->Get_cpw();
64 cph = parameters->Get_cph();
65 cpc = parameters->Get_cpc();
66
67 c0 = parameters->Get_c0();
68 c3 = parameters->Get_c3();
69 cs = parameters->Get_cs();
70 g0 = parameters->Get_g0(); // Skyrme-QMD
71 g0iso = parameters->Get_g0iso(); // Skyrme-QMD
72 gtau0 = parameters->Get_gtau0(); // Skyrme-QMD
73
74 // distance
75 c0w = 1.0/4.0/wl;
76 //c3w = 1.0/4.0/wl; //no need
77 c0sw = std::sqrt( c0w );
78 clw = 2.0 / std::sqrt ( 4.0 * pi * wl );
79
80 // graduate
81 c0g = - c0 / ( 2.0 * wl );
82 c3g = - c3 / ( 4.0 * wl ) * gamm;
83 csg = - cs / ( 2.0 * wl );
84 pag = gamm - 1;
85 pag_tau = eta - 1; // Skyrme-QMD
86 cg0 = - g0 / ( 2.0 * wl ); // Skyrme-QMD
87 cgtau0 = - gtau0 / ( 4.0 * wl ) * eta; // Skyrme-QMD
88 }
89
90
91
92 //G4LightIonQMDNucleus::~G4LightIonQMDNucleus()
93 //{
94 // ;
95 //}
96
97
98 G4LorentzVector G4LightIonQMDNucleus::Get4Momentum()
99 {
100 G4LorentzVector p( 0 );
101 std::vector< G4QMDParticipant* >::iterator it;
102 for ( it = participants.begin() ; it != participants.end() ; it++ )
103 p += (*it)->Get4Momentum();
104
105 return p;
106 }
107
108
109
110 G4int G4LightIonQMDNucleus::GetMassNumber()
111 {
112
113 G4int A = 0;
114 std::vector< G4QMDParticipant* >::iterator it;
115 for ( it = participants.begin() ; it != participants.end() ; it++ )
116 {
117 if ( (*it)->GetDefinition() == G4Proton::Proton()
118 || (*it)->GetDefinition() == G4Neutron::Neutron() )
119 A++;
120 }
121
122 if ( A == 0 ) {
123 throw G4HadronicException(__FILE__, __LINE__, "G4LightIonQMDNucleus has the mass number of 0!");
124 }
125
126 return A;
127 }
128
129
130
131 G4int G4LightIonQMDNucleus::GetAtomicNumber()
132 {
133 G4int Z = 0;
134 std::vector< G4QMDParticipant* >::iterator it;
135 for ( it = participants.begin() ; it != participants.end() ; it++ )
136 {
137 if ( (*it)->GetDefinition() == G4Proton::Proton() )
138 Z++;
139 }
140 return Z;
141 }
142
143
144
145 G4double G4LightIonQMDNucleus::GetNuclearMass()
146 {
147
148 G4double mass = G4NucleiProperties::GetNuclearMass( GetMassNumber() , GetAtomicNumber() );
149
150 if ( mass == 0.0 )
151 {
152
153 G4int Z = GetAtomicNumber();
154 G4int A = GetMassNumber();
155 G4int N = A - Z;
156
157 // Weizsacker-Bethe
158
159 G4double Av = 16*MeV;
160 G4double As = 17*MeV;
161 G4double Ac = 0.7*MeV;
162 G4double Asym = 23*MeV;
163
164 G4double BE = Av * A
165 - As * G4Pow::GetInstance()->A23 ( G4double ( A ) )
166 - Ac * Z*Z/G4Pow::GetInstance()->A13 ( G4double ( A ) )
167 - Asym * ( N - Z )* ( N - Z ) / A;
168
169 mass = Z * G4Proton::Proton()->GetPDGMass()
170 + N * G4Neutron::Neutron()->GetPDGMass()
171 - BE;
172
173 }
174
175 return mass;
176 }
177
178
179
180 void G4LightIonQMDNucleus::CalEnergyAndAngularMomentumInCM()
181 {
182
183 //G4cout << "CalEnergyAndAngularMomentumInCM " << this->GetAtomicNumber() << " " << GetMassNumber() << G4endl;
184
185 G4double gamma = Get4Momentum().gamma();
186 G4ThreeVector beta = Get4Momentum().v()/ Get4Momentum().e();
187
188 G4ThreeVector pcm0( 0.0 ) ;
189
190 G4int n = GetTotalNumberOfParticipant();
191 pcm.resize( n );
192
193 for ( G4int i= 0; i < n ; i++ )
194 {
195 G4ThreeVector p_i = GetParticipant( i )->GetMomentum();
196
197 G4double trans = gamma / ( gamma + 1.0 ) * p_i * beta;
198 pcm[i] = p_i - trans*beta;
199
200 pcm0 += pcm[i];
201 }
202
203 pcm0 = pcm0 / double ( n );
204
205 //G4cout << "pcm0 " << pcm0 << G4endl;
206
207 for ( G4int i= 0; i < n ; i++ )
208 {
209 pcm[i] += -pcm0;
210 //G4cout << "pcm " << i << " " << pcm[i] << G4endl;
211 }
212
213
214 G4double tmass = 0;
215 G4ThreeVector rcm0( 0.0 ) ;
216 rcm.resize( n );
217 es.resize( n );
218
219 // binding energy should be evaluated with a relativistic version: 20230308 by Y-H. Sato and A. Haga
220 for ( G4int i= 0; i < n ; i++ )
221 {
222 G4ThreeVector ri = GetParticipant( i )->GetPosition();
223 G4double trans = gamma / ( gamma + 1.0 ) * ri * beta;
224 G4double nucpote = GetNuclPotential( i );
225
226 es[i] = std::sqrt ( G4Pow::GetInstance()->powN ( GetParticipant( i )->GetMass() , 2 ) + pcm[i]*pcm[i] + 2.0*GetParticipant( i )->GetMass()*nucpote) - GetParticipant( i )->GetMass(); //R-JQMD
227
228 rcm[i] = ri + trans*beta;
229
230 rcm0 += rcm[i]*es[i];
231
232 tmass += es[i];
233 }
234
235 rcm0 = rcm0/tmass;
236
237 for ( G4int i= 0; i < n ; i++ )
238 {
239 rcm[i] += -rcm0;
240 //G4cout << "rcm " << i << " " << rcm[i] << G4endl;
241 }
242
243 // Angular momentum
244
245 G4ThreeVector rl ( 0.0 );
246 for ( G4int i= 0; i < n ; i++ )
247 {
248 rl += rcm[i].cross ( pcm[i] );
249 }
250
251 // DHW: move hbc outside of sqrt to get correct units
252 // jj = int ( std::sqrt ( rl*rl / hbc ) + 0.5 );
253
254 jj = int (std::sqrt(rl*rl)/hbc + 0.5);
255
256 // kinetic energy per nucleon in CM
257
258 /*
259 G4double totalMass = 0.0;
260 for ( G4int i= 0; i < n ; i++ )
261 {
262 // following two lines are equivalent
263 //totalMass += GetParticipant( i )->GetDefinition()->GetPDGMass()/GeV;
264 totalMass += GetParticipant( i )->GetMass();
265 }
266 */
267
268 //G4double kineticEnergyPerNucleon = ( std::accumulate ( es.begin() , es.end() , 0.0 ) - totalMass )/n;
269
270 // Total (not per nucleion ) Binding Energy
271 // relativistic version Y-H. Sato and A. Haga 20230309
272 G4double bindingEnergy = ( std::accumulate ( es.begin() , es.end() , 0.0 ) );
273
274 //G4cout << "n " << n << "totalpote " << totalpote << " " << potentialEnergy << " " << bindingEnergy << G4endl;
275 //G4cout << "KineticEnergyPerNucleon in GeV " << kineticEnergyPerNucleon << G4endl;
276 //G4cout << "KineticEnergySum in GeV " << std::accumulate ( es.begin() , es.end() , 0.0 ) - totalMass << G4endl;
277 //G4cout << "PotentialEnergy in GeV " << potentialEnergy << G4endl;
278 //G4cout << "BindingEnergy in GeV " << bindingEnergy << G4endl;
279 //G4cout << "G4BindingEnergy in GeV " << G4NucleiProperties::GetBindingEnergy( GetAtomicNumber() , GetMassNumber() )/GeV << G4endl;
280
281 excitationEnergy = bindingEnergy + G4NucleiProperties::GetBindingEnergy( GetMassNumber() , GetAtomicNumber() )/GeV;
282 if ( excitationEnergy < 0 ) excitationEnergy = 0.0;
283
284 }
285
286 // Get potential with a relativistic version added by Y-H. Sato and A. Haga 20230309
287 G4double G4LightIonQMDNucleus::GetNuclPotential( G4int i )
288 {
289 epsx = -20.0;
290 epscl = 0.0001; // coulomb term
291 irelcr = 1;
292 G4int n = GetTotalNumberOfParticipant();
293
294 G4double rhoa = 0.0;
295 G4double rho3 = 0.0;
296 G4double fsij_rhoa = 0.0; // Skyrme-QMD
297 // G4double fsij_rhos = 0.0; // Skyrme-QMD
298 G4double rho3_tau = 0.0; // Skyrme-QMD
299 G4double rhos = 0.0;
300 G4double rhoc = 0.0;
301
302
303 G4int icharge = GetParticipant(i)->GetChargeInUnitOfEplus();
304 G4int inuc = GetParticipant(i)->GetNuc();
305 G4int ibry = GetParticipant(i)->GetBaryonNumber();
306
307 G4ThreeVector ri = GetParticipant( i )->GetPosition();
308 G4LorentzVector p4i = GetParticipant( i )->Get4Momentum();
309
310 for ( G4int j = 0 ; j < n ; j ++ )
311 {
312 G4double cef = 1.0;
313 if (i == j)
314 {
315 cef = 0.0;
316 }
317
318
319 G4int jcharge = GetParticipant(j)->GetChargeInUnitOfEplus();
320 G4int jnuc = GetParticipant(j)->GetNuc();
321 G4int jbry = GetParticipant(j)->GetBaryonNumber();
322
323 G4ThreeVector rj = GetParticipant( j )->GetPosition();
324 G4LorentzVector p4j = GetParticipant( j )->Get4Momentum();
325
326 G4ThreeVector rij = ri - rj;
327 G4ThreeVector pij = (p4i - p4j).v();
328 G4LorentzVector p4ij = p4i - p4j;
329 G4ThreeVector bij = ( p4i + p4j ).boostVector();
330 G4double gammaij = ( p4i + p4j ).gamma();
331
332 //G4double eij = ( p4i + p4j ).e();
333
334 G4double rbrb = rij*bij;
335 // G4double bij2 = bij*bij;
336 G4double rij2 = rij*rij;
337 //G4double pij2 = pij*pij;
338
339
340 rbrb = irelcr * rbrb;
341 G4double gamma2_ij = gammaij*gammaij;
342
343 G4double rr2 = rij2 + gamma2_ij * rbrb*rbrb;
344
345 G4double expa1 = - (rij2 + gamma2_ij * rbrb*rbrb) * c0w;
346 G4double rh1;
347 if ( expa1 > epsx )
348 {
349 rh1 = G4Exp( expa1 );
350 }
351 else
352 {
353 rh1 = 0.0;
354 }
355
356 G4double rrs2 = (rij2 + gamma2_ij * rbrb*rbrb) + epscl;
357 G4double rrs = std::sqrt ( rrs2 );
358
359 G4double xerf = 0.0;
360 // T. K. add this protection. 5.8 is good enough for double
361 if ( rrs*c0sw < 5.8 ) {
362 //erf = G4RandStat::erf ( rrs*c0sw );
363 //Restore to CLHEP for avoiding compilation error in MT
364 //erf = CLHEP::HepStat::erf ( rrs*c0sw );
365 //Use cmath
366 #if defined WIN32-VC
367 xerf = CLHEP::HepStat::erf ( rrs*c0sw );
368 #else
369 xerf = std::erf ( rrs*c0sw );
370 #endif
371 } else {
372 xerf = 1.0;
373 }
374
375 G4double erfij = xerf/rrs;
376
377 G4double fsij = 3.0/(2*wl) - rr2/(2*wl)/(2*wl); // Add for Skyrme-QMD
378
379 rhoa += ibry*jbry*rh1*cef;
380 fsij_rhoa += fsij * ibry*jbry*rh1*cef; // Skyrme-QMD
381 rhoc += icharge*jcharge * erfij * cef;
382 rhos += ibry*jbry*rh1 * jnuc * inuc * cef
383 * ( 1 - 2 * std::abs ( jcharge - icharge ) )
384 * (1. - kappas * fsij);
385
386 //G4cout << i << " " << j << " " << ( - erfij ) << " " << clw << G4endl;
387
388
389 }
390
391 rho3 = G4Pow::GetInstance()->powA ( rhoa , gamm );
392 rho3_tau = G4Pow::GetInstance()->powA ( rhoa , eta );
393
394 G4double potential = c0 * rhoa
395 + c3 * rho3
396 + g0 * fsij_rhoa // Skyrme-QMD
397 // + g0iso * fsij_rhos // Skyrme-QMD
398 + gtau0 * rho3_tau // Skyrme-QMD
399 + cs * rhos
400 + cl * rhoc;
401
402 //G4cout << "n " << n << " " << rho3 << G4endl;
403 return potential;
404 }
405