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