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 | #include "G4QMDGroundStateNucleus.hh" |
---|
29 | |
---|
30 | #include "G4NucleiProperties.hh" |
---|
31 | |
---|
32 | #include "G4Proton.hh" |
---|
33 | #include "G4Neutron.hh" |
---|
34 | |
---|
35 | #include "Randomize.hh" |
---|
36 | |
---|
37 | G4QMDGroundStateNucleus::G4QMDGroundStateNucleus( G4int z , G4int a ) |
---|
38 | : r00 ( 1.124 ) // radius parameter for Woods-Saxon [fm] |
---|
39 | , r01 ( 0.5 ) // radius parameter for Woods-Saxon |
---|
40 | , saa ( 0.2 ) // diffuse parameter for initial Woods-Saxon shape |
---|
41 | , rada ( 0.9 ) // cutoff parameter |
---|
42 | , radb ( 0.3 ) // cutoff parameter |
---|
43 | , dsam ( 1.5 ) // minimum distance for same particle [fm] |
---|
44 | , ddif ( 1.0 ) // minimum distance for different particle |
---|
45 | , epse ( 0.000001 ) // torelance for energy in [GeV] |
---|
46 | { |
---|
47 | |
---|
48 | //std::cout << " G4QMDGroundStateNucleus( G4int z , G4int a ) Begin " << z << " " << a << std::endl; |
---|
49 | |
---|
50 | if ( z == 1 && a == 1 ) // Hydrogen Case or proton primary |
---|
51 | { |
---|
52 | SetParticipant( new G4QMDParticipant( G4Proton::Proton() , G4ThreeVector( 0.0 ) , G4ThreeVector( 0.0 ) ) ); |
---|
53 | return; |
---|
54 | } |
---|
55 | else if ( z == 0 && a == 1 ) // Neutron primary |
---|
56 | { |
---|
57 | SetParticipant( new G4QMDParticipant( G4Neutron::Neutron() , G4ThreeVector( 0.0 ) , G4ThreeVector( 0.0 ) ) ); |
---|
58 | return; |
---|
59 | } |
---|
60 | |
---|
61 | dsam2 = dsam*dsam; |
---|
62 | ddif2 = ddif*ddif; |
---|
63 | |
---|
64 | G4QMDParameters* parameters = G4QMDParameters::GetInstance(); |
---|
65 | |
---|
66 | hbc = parameters->Get_hbc(); |
---|
67 | gamm = parameters->Get_gamm(); |
---|
68 | cpw = parameters->Get_cpw(); |
---|
69 | cph = parameters->Get_cph(); |
---|
70 | epsx = parameters->Get_epsx(); |
---|
71 | cpc = parameters->Get_cpc(); |
---|
72 | |
---|
73 | cdp = parameters->Get_cdp(); |
---|
74 | c0p = parameters->Get_c0p(); |
---|
75 | c3p = parameters->Get_c3p(); |
---|
76 | csp = parameters->Get_csp(); |
---|
77 | clp = parameters->Get_clp(); |
---|
78 | |
---|
79 | edepth = 0.0; |
---|
80 | |
---|
81 | for ( int i = 0 ; i < a ; i++ ) |
---|
82 | { |
---|
83 | |
---|
84 | G4ParticleDefinition* pd; |
---|
85 | |
---|
86 | if ( i < z ) |
---|
87 | { |
---|
88 | pd = G4Proton::Proton(); |
---|
89 | } |
---|
90 | else |
---|
91 | { |
---|
92 | pd = G4Neutron::Neutron(); |
---|
93 | } |
---|
94 | |
---|
95 | G4ThreeVector p( 0.0 ); |
---|
96 | G4ThreeVector r( 0.0 ); |
---|
97 | G4QMDParticipant* aParticipant = new G4QMDParticipant( pd , p , r ); |
---|
98 | SetParticipant( aParticipant ); |
---|
99 | |
---|
100 | } |
---|
101 | |
---|
102 | G4double radious = r00 * std::pow ( double ( GetMassNumber() ) , 1.0/3.0 ); |
---|
103 | |
---|
104 | rt00 = radious - r01; |
---|
105 | radm = radious - rada * ( gamm - 1.0 ) + radb; |
---|
106 | rmax = 1.0 / ( 1.0 + std::exp ( -rt00/saa ) ); |
---|
107 | |
---|
108 | maxTrial = 1000; |
---|
109 | meanfield = new G4QMDMeanField(); |
---|
110 | meanfield->SetSystem( this ); |
---|
111 | |
---|
112 | //std::cout << "G4QMDGroundStateNucleus( G4int z , G4int a ) packNucleons Begin ( z , a ) ( " << z << ", " << a << " )" << std::endl; |
---|
113 | packNucleons(); |
---|
114 | //std::cout << "G4QMDGroundStateNucleus( G4int z , G4int a ) packNucleons End" << std::endl; |
---|
115 | |
---|
116 | delete meanfield; |
---|
117 | |
---|
118 | } |
---|
119 | |
---|
120 | |
---|
121 | |
---|
122 | void G4QMDGroundStateNucleus::packNucleons() |
---|
123 | { |
---|
124 | |
---|
125 | //std::cout << "G4QMDGroundStateNucleus::packNucleons" << std::endl; |
---|
126 | |
---|
127 | ebini = - G4NucleiProperties::GetBindingEnergy( GetMassNumber() , GetAtomicNumber() ) / GetMassNumber(); |
---|
128 | |
---|
129 | G4double ebin00 = ebini * 0.001; |
---|
130 | |
---|
131 | G4double ebin0 = 0.0; |
---|
132 | G4double ebin1 = 0.0; |
---|
133 | |
---|
134 | if ( GetMassNumber() != 4 ) |
---|
135 | { |
---|
136 | ebin0 = ( ebini - 0.5 ) * 0.001; |
---|
137 | ebin1 = ( ebini + 0.5 ) * 0.001; |
---|
138 | } |
---|
139 | else |
---|
140 | { |
---|
141 | ebin0 = ( ebini - 1.5 ) * 0.001; |
---|
142 | ebin1 = ( ebini + 1.5 ) * 0.001; |
---|
143 | } |
---|
144 | |
---|
145 | { |
---|
146 | G4int n0Try = 0; |
---|
147 | G4bool isThisOK = false; |
---|
148 | while ( n0Try < maxTrial ) |
---|
149 | { |
---|
150 | n0Try++; |
---|
151 | //std::cout << "TKDB packNucleons n0Try " << n0Try << std::endl; |
---|
152 | |
---|
153 | // Sampling Position |
---|
154 | |
---|
155 | //std::cout << "TKDB Sampling Position " << std::endl; |
---|
156 | |
---|
157 | G4bool areThesePsOK = false; |
---|
158 | G4int npTry = 0; |
---|
159 | while ( npTry < maxTrial ) |
---|
160 | { |
---|
161 | //std::cout << "TKDB Sampling Position npTry " << npTry << std::endl; |
---|
162 | npTry++; |
---|
163 | G4int i = 0; |
---|
164 | if ( samplingPosition( i ) ) |
---|
165 | { |
---|
166 | //std::cout << "packNucleons samplingPosition 0 succeed " << std::endl; |
---|
167 | for ( i = 1 ; i < GetMassNumber() ; i++ ) |
---|
168 | { |
---|
169 | //std::cout << "packNucleons samplingPosition " << i << " trying " << std::endl; |
---|
170 | if ( !( samplingPosition( i ) ) ) |
---|
171 | { |
---|
172 | //std::cout << "packNucleons samplingPosition " << i << " failed" << std::endl; |
---|
173 | break; |
---|
174 | } |
---|
175 | } |
---|
176 | if ( i == GetMassNumber() ) |
---|
177 | { |
---|
178 | //std::cout << "packNucleons samplingPosition all scucceed " << std::endl; |
---|
179 | areThesePsOK = true; |
---|
180 | break; |
---|
181 | } |
---|
182 | } |
---|
183 | } |
---|
184 | if ( areThesePsOK == false ) continue; |
---|
185 | |
---|
186 | //std::cout << "TKDB Sampling Position End" << std::endl; |
---|
187 | |
---|
188 | // Calculate Two-body quantities |
---|
189 | |
---|
190 | meanfield->Cal2BodyQuantities(); |
---|
191 | std::vector< G4double > rho_a ( GetMassNumber() , 0.0 ); |
---|
192 | std::vector< G4double > rho_s ( GetMassNumber() , 0.0 ); |
---|
193 | std::vector< G4double > rho_c ( GetMassNumber() , 0.0 ); |
---|
194 | |
---|
195 | rho_l.resize ( GetMassNumber() , 0.0 ); |
---|
196 | d_pot.resize ( GetMassNumber() , 0.0 ); |
---|
197 | |
---|
198 | for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
199 | { |
---|
200 | for ( G4int j = 0 ; j < GetMassNumber() ; j++ ) |
---|
201 | { |
---|
202 | |
---|
203 | rho_a[ i ] += meanfield->GetRHA( j , i ); |
---|
204 | G4int k = 0; |
---|
205 | if ( participants[j]->GetDefinition() != participants[i]->GetDefinition() ) |
---|
206 | { |
---|
207 | k = 1; |
---|
208 | } |
---|
209 | |
---|
210 | rho_s[ i ] += meanfield->GetRHA( j , i )*( 1.0 - 2.0 * k ); // OK? |
---|
211 | |
---|
212 | rho_c[ i ] += meanfield->GetRHE( j , i ); |
---|
213 | } |
---|
214 | |
---|
215 | } |
---|
216 | |
---|
217 | for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
218 | { |
---|
219 | rho_l[i] = cdp * ( rho_a[i] + 1.0 ); |
---|
220 | d_pot[i] = c0p * rho_a[i] |
---|
221 | + c3p * std::pow ( rho_a[i] , gamm ) |
---|
222 | + csp * rho_s[i] |
---|
223 | + clp * rho_c[i]; |
---|
224 | |
---|
225 | //std::cout << "d_pot[i] " << i << " " << d_pot[i] << std::endl; |
---|
226 | } |
---|
227 | |
---|
228 | |
---|
229 | // Sampling Momentum |
---|
230 | |
---|
231 | //std::cout << "TKDB Sampling Momentum " << std::endl; |
---|
232 | |
---|
233 | phase_g.clear(); |
---|
234 | phase_g.resize( GetMassNumber() , 0.0 ); |
---|
235 | |
---|
236 | //std::cout << "TKDB Sampling Momentum 1st " << std::endl; |
---|
237 | G4bool isThis1stMOK = false; |
---|
238 | G4int nmTry = 0; |
---|
239 | while ( nmTry < maxTrial ) |
---|
240 | { |
---|
241 | nmTry++; |
---|
242 | G4int i = 0; |
---|
243 | if ( samplingMomentum( i ) ) |
---|
244 | { |
---|
245 | isThis1stMOK = true; |
---|
246 | break; |
---|
247 | } |
---|
248 | } |
---|
249 | if ( isThis1stMOK == false ) continue; |
---|
250 | |
---|
251 | //std::cout << "TKDB Sampling Momentum 2nd so on" << std::endl; |
---|
252 | |
---|
253 | G4bool areTheseMsOK = true; |
---|
254 | nmTry = 0; |
---|
255 | while ( nmTry < maxTrial ) |
---|
256 | { |
---|
257 | nmTry++; |
---|
258 | G4int i = 0; |
---|
259 | for ( i = 1 ; i < GetMassNumber() ; i++ ) |
---|
260 | { |
---|
261 | //std::cout << "TKDB packNucleons samplingMomentum try " << i << std::endl; |
---|
262 | if ( !( samplingMomentum( i ) ) ) |
---|
263 | { |
---|
264 | //std::cout << "TKDB packNucleons samplingMomentum " << i << " failed" << std::endl; |
---|
265 | areTheseMsOK = false; |
---|
266 | break; |
---|
267 | } |
---|
268 | } |
---|
269 | if ( i == GetMassNumber() ) |
---|
270 | { |
---|
271 | areTheseMsOK = true; |
---|
272 | } |
---|
273 | |
---|
274 | break; |
---|
275 | } |
---|
276 | if ( areTheseMsOK == false ) continue; |
---|
277 | |
---|
278 | // Kill Angluar Momentum |
---|
279 | |
---|
280 | //std::cout << "TKDB Sampling Kill Angluar Momentum " << std::endl; |
---|
281 | |
---|
282 | killCMMotionAndAngularM(); |
---|
283 | |
---|
284 | |
---|
285 | // Check Binding Energy |
---|
286 | |
---|
287 | //std::cout << "packNucleons Check Binding Energy Begin " << std::endl; |
---|
288 | |
---|
289 | G4double ekinal = 0.0; |
---|
290 | for ( int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
291 | { |
---|
292 | ekinal += participants[i]->GetKineticEnergy(); |
---|
293 | } |
---|
294 | |
---|
295 | meanfield->Cal2BodyQuantities(); |
---|
296 | |
---|
297 | G4double totalPotentialE = meanfield->GetTotalPotential(); |
---|
298 | G4double ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() ); |
---|
299 | |
---|
300 | //std::cout << "packNucleons totalPotentialE " << totalPotentialE << std::endl; |
---|
301 | //std::cout << "packNucleons ebinal " << ebinal << std::endl; |
---|
302 | //std::cout << "packNucleons ekinal " << ekinal << std::endl; |
---|
303 | |
---|
304 | if ( ebinal < ebin0 || ebinal > ebin1 ) |
---|
305 | { |
---|
306 | //std::cout << "packNucleons ebin0 " << ebin0 << std::endl; |
---|
307 | //std::cout << "packNucleons ebin1 " << ebin1 << std::endl; |
---|
308 | //std::cout << "packNucleons ebinal " << ebinal << std::endl; |
---|
309 | //std::cout << "packNucleons Check Binding Energy Failed " << std::endl; |
---|
310 | continue; |
---|
311 | } |
---|
312 | |
---|
313 | //std::cout << "packNucleons Check Binding Energy End = OK " << std::endl; |
---|
314 | |
---|
315 | |
---|
316 | // Energy Adujstment |
---|
317 | |
---|
318 | G4double dtc = 1.0; |
---|
319 | G4double frg = -0.1; |
---|
320 | G4double rdf0 = 0.5; |
---|
321 | |
---|
322 | G4double edif0 = ebinal - ebin00; |
---|
323 | |
---|
324 | G4double cfrc = 0.0; |
---|
325 | if ( 0 < edif0 ) |
---|
326 | { |
---|
327 | cfrc = frg; |
---|
328 | } |
---|
329 | else |
---|
330 | { |
---|
331 | cfrc = -frg; |
---|
332 | } |
---|
333 | |
---|
334 | G4int ifrc = 1; |
---|
335 | |
---|
336 | G4int neaTry = 0; |
---|
337 | |
---|
338 | G4bool isThisEAOK = false; |
---|
339 | while ( neaTry < maxTrial ) |
---|
340 | { |
---|
341 | neaTry++; |
---|
342 | |
---|
343 | G4double edif = ebinal - ebin00; |
---|
344 | |
---|
345 | //std::cout << "TKDB edif " << edif << std::endl; |
---|
346 | if ( std::abs ( edif ) < epse ) |
---|
347 | { |
---|
348 | |
---|
349 | isThisEAOK = true; |
---|
350 | //std::cout << "isThisEAOK " << isThisEAOK << std::endl; |
---|
351 | break; |
---|
352 | } |
---|
353 | |
---|
354 | G4int jfrc = 0; |
---|
355 | if ( edif < 0.0 ) |
---|
356 | { |
---|
357 | jfrc = 1; |
---|
358 | } |
---|
359 | else |
---|
360 | { |
---|
361 | jfrc = -1; |
---|
362 | } |
---|
363 | |
---|
364 | if ( jfrc != ifrc ) |
---|
365 | { |
---|
366 | cfrc = -rdf0 * cfrc; |
---|
367 | dtc = rdf0 * dtc; |
---|
368 | } |
---|
369 | |
---|
370 | if ( jfrc == ifrc && std::abs( edif0 ) < std::abs( edif ) ) |
---|
371 | { |
---|
372 | cfrc = -rdf0 * cfrc; |
---|
373 | dtc = rdf0 * dtc; |
---|
374 | } |
---|
375 | |
---|
376 | ifrc = jfrc; |
---|
377 | edif0 = edif; |
---|
378 | |
---|
379 | meanfield->CalGraduate(); |
---|
380 | |
---|
381 | for ( int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
382 | { |
---|
383 | G4ThreeVector ri = participants[i]->GetPosition(); |
---|
384 | G4ThreeVector p_i = participants[i]->GetMomentum(); |
---|
385 | |
---|
386 | ri += dtc * ( meanfield->GetFFr(i) - cfrc * ( meanfield->GetFFp(i) ) ); |
---|
387 | p_i += dtc * ( meanfield->GetFFp(i) + cfrc * ( meanfield->GetFFr(i) ) ); |
---|
388 | |
---|
389 | participants[i]->SetPosition( ri ); |
---|
390 | participants[i]->SetMomentum( p_i ); |
---|
391 | } |
---|
392 | |
---|
393 | ekinal = 0.0; |
---|
394 | |
---|
395 | for ( int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
396 | { |
---|
397 | ekinal += participants[i]->GetKineticEnergy(); |
---|
398 | } |
---|
399 | |
---|
400 | meanfield->Cal2BodyQuantities(); |
---|
401 | totalPotentialE = meanfield->GetTotalPotential(); |
---|
402 | |
---|
403 | ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() ); |
---|
404 | |
---|
405 | } |
---|
406 | //std::cout << "isThisEAOK " << isThisEAOK << std::endl; |
---|
407 | if ( isThisEAOK == false ) continue; |
---|
408 | |
---|
409 | isThisOK = true; |
---|
410 | //std::cout << "isThisOK " << isThisOK << std::endl; |
---|
411 | break; |
---|
412 | |
---|
413 | } |
---|
414 | |
---|
415 | if ( isThisOK == false ) |
---|
416 | { |
---|
417 | std::cout << "GroundStateNucleus state cannot be created. Try again with another parameters." << std::endl; |
---|
418 | } |
---|
419 | |
---|
420 | //std::cout << "packNucleons End" << std::endl; |
---|
421 | return; |
---|
422 | |
---|
423 | } |
---|
424 | |
---|
425 | |
---|
426 | |
---|
427 | |
---|
428 | |
---|
429 | // Start packing |
---|
430 | // position |
---|
431 | |
---|
432 | G4int n0Try = 0; |
---|
433 | |
---|
434 | while ( n0Try < maxTrial ) |
---|
435 | { |
---|
436 | if ( samplingPosition( 0 ) ) |
---|
437 | { |
---|
438 | G4int i = 0; |
---|
439 | for ( i = 1 ; i < GetMassNumber() ; i++ ) |
---|
440 | { |
---|
441 | if ( !( samplingPosition( i ) ) ) |
---|
442 | { |
---|
443 | break; |
---|
444 | } |
---|
445 | } |
---|
446 | if ( i == GetMassNumber() ) break; |
---|
447 | } |
---|
448 | n0Try++; |
---|
449 | } |
---|
450 | |
---|
451 | if ( n0Try > maxTrial ) |
---|
452 | { |
---|
453 | std::cout << "GroundStateNucleus state cannot be created. Try again with another parameters." << std::endl; |
---|
454 | return; |
---|
455 | } |
---|
456 | |
---|
457 | meanfield->Cal2BodyQuantities(); |
---|
458 | std::vector< G4double > rho_a ( GetMassNumber() , 0.0 ); |
---|
459 | std::vector< G4double > rho_s ( GetMassNumber() , 0.0 ); |
---|
460 | std::vector< G4double > rho_c ( GetMassNumber() , 0.0 ); |
---|
461 | |
---|
462 | rho_l.resize ( GetMassNumber() , 0.0 ); |
---|
463 | d_pot.resize ( GetMassNumber() , 0.0 ); |
---|
464 | |
---|
465 | for ( int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
466 | { |
---|
467 | for ( int j = 0 ; j < GetMassNumber() ; j++ ) |
---|
468 | { |
---|
469 | |
---|
470 | rho_a[ i ] += meanfield->GetRHA( j , i ); |
---|
471 | G4int k = 0; |
---|
472 | if ( participants[i]->GetDefinition() != participants[i]->GetDefinition() ) |
---|
473 | { |
---|
474 | k = 1; |
---|
475 | } |
---|
476 | |
---|
477 | rho_s[ i ] += meanfield->GetRHA( j , i )*( 1.0 - 2.0 * k ); // OK? |
---|
478 | |
---|
479 | rho_c[ i ] += meanfield->GetRHE( j , i ); |
---|
480 | } |
---|
481 | } |
---|
482 | |
---|
483 | for ( int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
484 | { |
---|
485 | rho_l[i] = cdp * ( rho_a[i] + 1.0 ); |
---|
486 | d_pot[i] = c0p * rho_a[i] |
---|
487 | + c3p * std::pow ( rho_a[i] , gamm ) |
---|
488 | + csp * rho_s[i] |
---|
489 | + clp * rho_c[i]; |
---|
490 | } |
---|
491 | |
---|
492 | |
---|
493 | |
---|
494 | |
---|
495 | |
---|
496 | |
---|
497 | // momentum |
---|
498 | |
---|
499 | phase_g.resize( GetMassNumber() , 0.0 ); |
---|
500 | |
---|
501 | //G4int i = 0; |
---|
502 | samplingMomentum( 0 ); |
---|
503 | |
---|
504 | G4int n1Try = 0; |
---|
505 | //G4int maxTry = 1000; |
---|
506 | |
---|
507 | while ( n1Try < maxTrial ) |
---|
508 | { |
---|
509 | if ( samplingPosition( 0 ) ) |
---|
510 | { |
---|
511 | G4int i = 0; |
---|
512 | G4bool isThisOK = true; |
---|
513 | for ( i = 1 ; i < GetMassNumber() ; i++ ) |
---|
514 | { |
---|
515 | if ( !( samplingMomentum( i ) ) ) |
---|
516 | { |
---|
517 | isThisOK = false; |
---|
518 | break; |
---|
519 | } |
---|
520 | } |
---|
521 | if ( isThisOK == true ) break; |
---|
522 | //if ( i == GetMassNumber() ) break; |
---|
523 | } |
---|
524 | n1Try++; |
---|
525 | } |
---|
526 | |
---|
527 | if ( n1Try > maxTrial ) |
---|
528 | { |
---|
529 | std::cout << "GroundStateNucleus state cannot be created. Try again with another parameters." << std::endl; |
---|
530 | return; |
---|
531 | } |
---|
532 | |
---|
533 | // |
---|
534 | |
---|
535 | // Shift nucleus to thier CM frame and kill angular momentum |
---|
536 | killCMMotionAndAngularM(); |
---|
537 | |
---|
538 | // Check binding energy |
---|
539 | |
---|
540 | G4double ekinal = 0.0; |
---|
541 | for ( int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
542 | { |
---|
543 | ekinal += participants[i]->GetKineticEnergy(); |
---|
544 | } |
---|
545 | |
---|
546 | meanfield->Cal2BodyQuantities(); |
---|
547 | G4double totalPotentialE = meanfield->GetTotalPotential(); |
---|
548 | G4double ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() ); |
---|
549 | |
---|
550 | if ( ebinal < ebin0 || ebinal > ebin1 ) |
---|
551 | { |
---|
552 | // Retry From Position |
---|
553 | } |
---|
554 | |
---|
555 | |
---|
556 | // Adjust by frictional cooling or heating |
---|
557 | |
---|
558 | G4double dtc = 1.0; |
---|
559 | G4double frg = -0.1; |
---|
560 | G4double rdf0 = 0.5; |
---|
561 | |
---|
562 | G4double edif0 = ebinal - ebin00; |
---|
563 | |
---|
564 | G4double cfrc = 0.0; |
---|
565 | if ( 0 < edif0 ) |
---|
566 | { |
---|
567 | cfrc = frg; |
---|
568 | } |
---|
569 | else |
---|
570 | { |
---|
571 | cfrc = -frg; |
---|
572 | } |
---|
573 | |
---|
574 | G4int ifrc = 1; |
---|
575 | |
---|
576 | G4int ntryACH = 0; |
---|
577 | |
---|
578 | G4bool isThisOK = false; |
---|
579 | while ( ntryACH < maxTrial ) |
---|
580 | { |
---|
581 | |
---|
582 | G4double edif = ebinal - ebin00; |
---|
583 | if ( std::abs ( edif ) < epse ) |
---|
584 | { |
---|
585 | isThisOK = true; |
---|
586 | break; |
---|
587 | } |
---|
588 | |
---|
589 | G4int jfrc = 0; |
---|
590 | if ( edif < 0.0 ) |
---|
591 | { |
---|
592 | jfrc = 1; |
---|
593 | } |
---|
594 | else |
---|
595 | { |
---|
596 | jfrc = -1; |
---|
597 | } |
---|
598 | |
---|
599 | if ( jfrc != ifrc ) |
---|
600 | { |
---|
601 | cfrc = -rdf0 * cfrc; |
---|
602 | dtc = rdf0 * dtc; |
---|
603 | } |
---|
604 | |
---|
605 | if ( jfrc == ifrc && std::abs( edif0 ) < std::abs( edif ) ) |
---|
606 | { |
---|
607 | cfrc = -rdf0 * cfrc; |
---|
608 | dtc = rdf0 * dtc; |
---|
609 | } |
---|
610 | |
---|
611 | ifrc = jfrc; |
---|
612 | edif0 = edif; |
---|
613 | |
---|
614 | meanfield->CalGraduate(); |
---|
615 | |
---|
616 | for ( int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
617 | { |
---|
618 | G4ThreeVector ri = participants[i]->GetPosition(); |
---|
619 | G4ThreeVector p_i = participants[i]->GetMomentum(); |
---|
620 | |
---|
621 | ri += dtc * ( meanfield->GetFFr(i) - cfrc * ( meanfield->GetFFp(i) ) ); |
---|
622 | p_i += dtc * ( meanfield->GetFFp(i) + cfrc * ( meanfield->GetFFr(i) ) ); |
---|
623 | |
---|
624 | participants[i]->SetPosition( ri ); |
---|
625 | participants[i]->SetMomentum( p_i ); |
---|
626 | } |
---|
627 | |
---|
628 | ekinal = 0.0; |
---|
629 | |
---|
630 | for ( int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
631 | { |
---|
632 | ekinal += participants[i]->GetKineticEnergy(); |
---|
633 | } |
---|
634 | |
---|
635 | meanfield->Cal2BodyQuantities(); |
---|
636 | totalPotentialE = meanfield->GetTotalPotential(); |
---|
637 | |
---|
638 | ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() ); |
---|
639 | |
---|
640 | |
---|
641 | ntryACH++; |
---|
642 | } |
---|
643 | |
---|
644 | if ( isThisOK ) |
---|
645 | { |
---|
646 | return; |
---|
647 | } |
---|
648 | |
---|
649 | } |
---|
650 | |
---|
651 | |
---|
652 | G4bool G4QMDGroundStateNucleus::samplingPosition( G4int i ) |
---|
653 | { |
---|
654 | |
---|
655 | G4bool result = false; |
---|
656 | |
---|
657 | G4int nTry = 0; |
---|
658 | |
---|
659 | while ( nTry < maxTrial ) |
---|
660 | { |
---|
661 | |
---|
662 | //std::cout << "samplingPosition i th particle, nTtry " << i << " " << nTry << std::endl; |
---|
663 | |
---|
664 | G4double rwod = -1.0; |
---|
665 | G4double rrr = 0.0; |
---|
666 | |
---|
667 | G4double rx = 0.0; |
---|
668 | G4double ry = 0.0; |
---|
669 | G4double rz = 0.0; |
---|
670 | |
---|
671 | while ( G4UniformRand() * rmax > rwod ) |
---|
672 | { |
---|
673 | |
---|
674 | G4double rsqr = 10.0; |
---|
675 | while ( rsqr > 1.0 ) |
---|
676 | { |
---|
677 | rx = 1.0 - 2.0 * G4UniformRand(); |
---|
678 | ry = 1.0 - 2.0 * G4UniformRand(); |
---|
679 | rz = 1.0 - 2.0 * G4UniformRand(); |
---|
680 | rsqr = rx*rx + ry*ry + rz*rz; |
---|
681 | } |
---|
682 | rrr = radm * std::sqrt ( rsqr ); |
---|
683 | rwod = 1.0 / ( 1.0 + std::exp ( ( rrr - rt00 ) / saa ) ); |
---|
684 | |
---|
685 | } |
---|
686 | |
---|
687 | participants[i]->SetPosition( G4ThreeVector( rx , ry , rz )*radm ); |
---|
688 | |
---|
689 | if ( i == 0 ) |
---|
690 | { |
---|
691 | result = true; |
---|
692 | return result; |
---|
693 | } |
---|
694 | |
---|
695 | // i > 1 ( Second Particle or later ) |
---|
696 | // Check Distance to others |
---|
697 | |
---|
698 | G4bool isThisOK = true; |
---|
699 | for ( G4int j = 0 ; j < i ; j++ ) |
---|
700 | { |
---|
701 | |
---|
702 | G4double r2 = participants[j]->GetPosition().diff2( participants[i]->GetPosition() ); |
---|
703 | G4double dmin2 = 0.0; |
---|
704 | |
---|
705 | if ( participants[j]->GetDefinition() == participants[i]->GetDefinition() ) |
---|
706 | { |
---|
707 | dmin2 = dsam2; |
---|
708 | } |
---|
709 | else |
---|
710 | { |
---|
711 | dmin2 = ddif2; |
---|
712 | } |
---|
713 | |
---|
714 | //std::cout << "distance between j and i " << j << " " << i << " " << r2 << " " << dmin2 << std::endl; |
---|
715 | if ( r2 < dmin2 ) |
---|
716 | { |
---|
717 | isThisOK = false; |
---|
718 | break; |
---|
719 | } |
---|
720 | |
---|
721 | } |
---|
722 | |
---|
723 | if ( isThisOK == true ) |
---|
724 | { |
---|
725 | result = true; |
---|
726 | return result; |
---|
727 | } |
---|
728 | |
---|
729 | nTry++; |
---|
730 | |
---|
731 | } |
---|
732 | |
---|
733 | // Here return "false" |
---|
734 | return result; |
---|
735 | } |
---|
736 | |
---|
737 | |
---|
738 | |
---|
739 | G4bool G4QMDGroundStateNucleus::samplingMomentum( G4int i ) |
---|
740 | { |
---|
741 | |
---|
742 | //std::cout << "TKDB samplingMomentum for " << i << std::endl; |
---|
743 | |
---|
744 | G4bool result = false; |
---|
745 | |
---|
746 | G4double pfm = hbc * std::pow ( ( 3.0 / 2.0 * pi*pi * rho_l[i] ) , 1./3. ); |
---|
747 | |
---|
748 | if ( 10 < GetMassNumber() && -5.5 < ebini ) |
---|
749 | { |
---|
750 | pfm = pfm * ( 1.0 + 0.2 * std::sqrt( std::abs( 8.0 + ebini ) / 8.0 ) ); |
---|
751 | } |
---|
752 | |
---|
753 | //std::cout << "TKDB samplingMomentum pfm " << pfm << std::endl; |
---|
754 | |
---|
755 | std::vector< G4double > phase; |
---|
756 | phase.resize( i+1 ); // i start from 0 |
---|
757 | |
---|
758 | G4int ntry = 0; |
---|
759 | // 710 |
---|
760 | while ( ntry < maxTrial ) |
---|
761 | { |
---|
762 | |
---|
763 | //std::cout << " TKDB ntry " << ntry << std::endl; |
---|
764 | ntry++; |
---|
765 | |
---|
766 | G4double ke = DBL_MAX; |
---|
767 | |
---|
768 | G4int tkdb_i =0; |
---|
769 | // 700 |
---|
770 | while ( ke + d_pot [i] > edepth ) |
---|
771 | { |
---|
772 | |
---|
773 | G4double psqr = 10.0; |
---|
774 | G4double px = 0.0; |
---|
775 | G4double py = 0.0; |
---|
776 | G4double pz = 0.0; |
---|
777 | |
---|
778 | while ( psqr > 1.0 ) |
---|
779 | { |
---|
780 | px = 1.0 - 2.0*G4UniformRand(); |
---|
781 | py = 1.0 - 2.0*G4UniformRand(); |
---|
782 | pz = 1.0 - 2.0*G4UniformRand(); |
---|
783 | |
---|
784 | psqr = px*px + py*py + pz*pz; |
---|
785 | } |
---|
786 | |
---|
787 | G4ThreeVector p ( px , py , pz ); |
---|
788 | p = pfm * p; |
---|
789 | participants[i]->SetMomentum( p ); |
---|
790 | G4LorentzVector p4 = participants[i]->Get4Momentum(); |
---|
791 | //ke = p4.e() - p4.restMass(); |
---|
792 | ke = participants[i]->GetKineticEnergy(); |
---|
793 | |
---|
794 | |
---|
795 | tkdb_i++; |
---|
796 | if ( tkdb_i > maxTrial ) return result; // return false |
---|
797 | |
---|
798 | } |
---|
799 | |
---|
800 | //std::cout << "TKDB ke d_pot[i] " << ke << " " << d_pot[i] << std::endl; |
---|
801 | |
---|
802 | |
---|
803 | if ( i == 0 ) |
---|
804 | { |
---|
805 | result = true; |
---|
806 | return result; |
---|
807 | } |
---|
808 | |
---|
809 | G4bool isThisOK = true; |
---|
810 | |
---|
811 | // Check Pauli principle |
---|
812 | |
---|
813 | phase[ i ] = 0.0; |
---|
814 | |
---|
815 | //std::cout << "TKDB Check Pauli Principle " << i << std::endl; |
---|
816 | |
---|
817 | for ( G4int j = 0 ; j < i ; j++ ) |
---|
818 | { |
---|
819 | phase[ j ] = 0.0; |
---|
820 | //std::cout << "TKDB Check Pauli Principle i , j " << i << " , " << j << std::endl; |
---|
821 | G4double expa = 0.0; |
---|
822 | if ( participants[j]->GetDefinition() == participants[i]->GetDefinition() ) |
---|
823 | { |
---|
824 | |
---|
825 | expa = - meanfield->GetRR2(i,j) * cpw; |
---|
826 | |
---|
827 | if ( expa > epsx ) |
---|
828 | { |
---|
829 | G4ThreeVector p_i = participants[i]->GetMomentum(); |
---|
830 | G4ThreeVector pj = participants[j]->GetMomentum(); |
---|
831 | G4double dist2_p = p_i.diff2( pj ); |
---|
832 | |
---|
833 | dist2_p = dist2_p*cph; |
---|
834 | expa = expa - dist2_p; |
---|
835 | |
---|
836 | if ( expa > epsx ) |
---|
837 | { |
---|
838 | |
---|
839 | phase[j] = std::exp ( expa ); |
---|
840 | |
---|
841 | if ( phase[j] * cpc > 0.2 ) |
---|
842 | { |
---|
843 | /* |
---|
844 | std::cout << "TKDB Check Pauli Principle A i , j " << i << " , " << j << std::endl; |
---|
845 | std::cout << "TKDB Check Pauli Principle phase[j] " << phase[j] << std::endl; |
---|
846 | std::cout << "TKDB Check Pauli Principle phase[j]*cpc > 0.2 " << phase[j]*cpc << std::endl; |
---|
847 | */ |
---|
848 | isThisOK = false; |
---|
849 | break; |
---|
850 | } |
---|
851 | if ( ( phase_g[j] + phase[j] ) * cpc > 0.5 ) |
---|
852 | { |
---|
853 | /* |
---|
854 | std::cout << "TKDB Check Pauli Principle B i , j " << i << " , " << j << std::endl; |
---|
855 | std::cout << "TKDB Check Pauli Principle B phase_g[j] " << phase_g[j] << std::endl; |
---|
856 | std::cout << "TKDB Check Pauli Principle B phase[j] " << phase[j] << std::endl; |
---|
857 | std::cout << "TKDB Check Pauli Principle B phase_g[j] + phase[j] ) * cpc > 0.5 " << ( phase_g[j] + phase[j] ) * cpc << std::endl; |
---|
858 | */ |
---|
859 | isThisOK = false; |
---|
860 | break; |
---|
861 | } |
---|
862 | |
---|
863 | phase[i] += phase[j]; |
---|
864 | if ( phase[i] * cpc > 0.3 ) |
---|
865 | { |
---|
866 | /* |
---|
867 | std::cout << "TKDB Check Pauli Principle C i , j " << i << " , " << j << std::endl; |
---|
868 | std::cout << "TKDB Check Pauli Principle C phase[i] " << phase[i] << std::endl; |
---|
869 | std::cout << "TKDB Check Pauli Principle C phase[i] * cpc > 0.3 " << phase[i] * cpc << std::endl; |
---|
870 | */ |
---|
871 | isThisOK = false; |
---|
872 | break; |
---|
873 | } |
---|
874 | |
---|
875 | //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl; |
---|
876 | |
---|
877 | } |
---|
878 | else |
---|
879 | { |
---|
880 | //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl; |
---|
881 | } |
---|
882 | |
---|
883 | } |
---|
884 | else |
---|
885 | { |
---|
886 | //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl; |
---|
887 | } |
---|
888 | |
---|
889 | } |
---|
890 | else |
---|
891 | { |
---|
892 | //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl; |
---|
893 | } |
---|
894 | |
---|
895 | } |
---|
896 | |
---|
897 | if ( isThisOK == true ) |
---|
898 | { |
---|
899 | |
---|
900 | phase_g[i] = phase[i]; |
---|
901 | |
---|
902 | for ( int j = 0 ; j < i ; j++ ) |
---|
903 | { |
---|
904 | phase_g[j] += phase[j]; |
---|
905 | } |
---|
906 | |
---|
907 | result = true; |
---|
908 | return result; |
---|
909 | } |
---|
910 | |
---|
911 | } |
---|
912 | |
---|
913 | return result; |
---|
914 | |
---|
915 | } |
---|
916 | |
---|
917 | |
---|
918 | |
---|
919 | void G4QMDGroundStateNucleus::killCMMotionAndAngularM() |
---|
920 | { |
---|
921 | |
---|
922 | // CalEnergyAndAngularMomentumInCM(); |
---|
923 | |
---|
924 | //std::vector< G4ThreeVector > p ( GetMassNumber() , 0.0 ); |
---|
925 | //std::vector< G4ThreeVector > r ( GetMassNumber() , 0.0 ); |
---|
926 | |
---|
927 | // Move to cm system |
---|
928 | |
---|
929 | G4ThreeVector pcm ( 0.0 ); |
---|
930 | G4ThreeVector rcm ( 0.0 ); |
---|
931 | G4double sumMass = 0.0; |
---|
932 | |
---|
933 | for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
934 | { |
---|
935 | pcm += participants[i]->GetMomentum(); |
---|
936 | rcm += participants[i]->GetPosition() * participants[i]->GetMass(); |
---|
937 | sumMass += participants[i]->GetMass(); |
---|
938 | } |
---|
939 | |
---|
940 | pcm = pcm/GetMassNumber(); |
---|
941 | rcm = rcm/sumMass; |
---|
942 | |
---|
943 | for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
944 | { |
---|
945 | participants[i]->SetMomentum( participants[i]->GetMomentum() - pcm ); |
---|
946 | participants[i]->SetPosition( participants[i]->GetPosition() - rcm ); |
---|
947 | } |
---|
948 | |
---|
949 | // kill the angular momentum |
---|
950 | |
---|
951 | G4ThreeVector ll ( 0.0 ); |
---|
952 | for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
953 | { |
---|
954 | ll += participants[i]->GetPosition().cross ( participants[i]->GetMomentum() ); |
---|
955 | } |
---|
956 | |
---|
957 | G4double rr[3][3]; |
---|
958 | G4double ss[3][3]; |
---|
959 | for ( G4int i = 0 ; i < 3 ; i++ ) |
---|
960 | { |
---|
961 | for ( G4int j = 0 ; j < 3 ; j++ ) |
---|
962 | { |
---|
963 | rr [i][j] = 0.0; |
---|
964 | |
---|
965 | if ( i == j ) |
---|
966 | { |
---|
967 | ss [i][j] = 1.0; |
---|
968 | } |
---|
969 | else |
---|
970 | { |
---|
971 | ss [i][j] = 0.0; |
---|
972 | } |
---|
973 | } |
---|
974 | } |
---|
975 | |
---|
976 | for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
977 | { |
---|
978 | G4ThreeVector r = participants[i]->GetPosition(); |
---|
979 | rr[0][0] += r.y() * r.y() + r.z() * r.z(); |
---|
980 | rr[1][0] += - r.y() * r.x(); |
---|
981 | rr[2][0] += - r.z() * r.x(); |
---|
982 | rr[0][1] += - r.x() * r.y(); |
---|
983 | rr[1][1] += r.z() * r.z() + r.x() * r.x(); |
---|
984 | rr[2][1] += - r.z() * r.y(); |
---|
985 | rr[2][0] += - r.x() * r.z(); |
---|
986 | rr[2][1] += - r.y() * r.z(); |
---|
987 | rr[2][2] += r.x() * r.x() + r.y() * r.y(); |
---|
988 | } |
---|
989 | |
---|
990 | for ( G4int i = 0 ; i < 3 ; i++ ) |
---|
991 | { |
---|
992 | G4double x = rr [i][i]; |
---|
993 | for ( G4int j = 0 ; j < 3 ; j++ ) |
---|
994 | { |
---|
995 | rr[i][j] = rr[i][j] / x; |
---|
996 | ss[i][j] = ss[i][j] / x; |
---|
997 | } |
---|
998 | for ( G4int j = 0 ; j < 3 ; j++ ) |
---|
999 | { |
---|
1000 | if ( j != i ) |
---|
1001 | { |
---|
1002 | G4double y = rr [j][i]; |
---|
1003 | for ( G4int k = 0 ; k < 3 ; k++ ) |
---|
1004 | { |
---|
1005 | rr[j][k] += -y * rr[i][k]; |
---|
1006 | ss[j][k] += -y * ss[i][k]; |
---|
1007 | } |
---|
1008 | } |
---|
1009 | } |
---|
1010 | } |
---|
1011 | |
---|
1012 | G4double opl[3]; |
---|
1013 | G4double rll[3]; |
---|
1014 | |
---|
1015 | rll[0] = ll.x(); |
---|
1016 | rll[1] = ll.y(); |
---|
1017 | rll[2] = ll.z(); |
---|
1018 | |
---|
1019 | for ( G4int i = 0 ; i < 3 ; i++ ) |
---|
1020 | { |
---|
1021 | opl[i] = 0.0; |
---|
1022 | |
---|
1023 | for ( G4int j = 0 ; j < 3 ; j++ ) |
---|
1024 | { |
---|
1025 | opl[i] += ss[i][j]*rll[j]; |
---|
1026 | } |
---|
1027 | } |
---|
1028 | |
---|
1029 | for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) |
---|
1030 | { |
---|
1031 | G4ThreeVector p_i = participants[i]->GetMomentum() ; |
---|
1032 | G4ThreeVector ri = participants[i]->GetPosition() ; |
---|
1033 | G4ThreeVector opl_v ( opl[0] , opl[1] , opl[2] ); |
---|
1034 | |
---|
1035 | p_i += ri.cross(opl_v); |
---|
1036 | |
---|
1037 | participants[i]->SetMomentum( p_i ); |
---|
1038 | } |
---|
1039 | |
---|
1040 | } |
---|