source: trunk/source/processes/hadronic/models/qmd/src/G4QMDGroundStateNucleus.cc@ 1201

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 26.4 KB
RevLine 
[819]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//
[962]26// 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
27//
[819]28#include "G4QMDGroundStateNucleus.hh"
29
30#include "G4NucleiProperties.hh"
31
32#include "G4Proton.hh"
33#include "G4Neutron.hh"
34
35#include "Randomize.hh"
36
37G4QMDGroundStateNucleus::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
[1196]50 if ( z == 1 && a == 1 ) // Hydrogen Case or proton primary
[819]51 {
52 SetParticipant( new G4QMDParticipant( G4Proton::Proton() , G4ThreeVector( 0.0 ) , G4ThreeVector( 0.0 ) ) );
53 return;
54 }
[1196]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 }
[819]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
122void G4QMDGroundStateNucleus::packNucleons()
123{
124
125 //std::cout << "G4QMDGroundStateNucleus::packNucleons" << std::endl;
126
[962]127 ebini = - G4NucleiProperties::GetBindingEnergy( GetMassNumber() , GetAtomicNumber() ) / GetMassNumber();
[819]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
652G4bool 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
739G4bool 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
919void 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}
Note: See TracBrowser for help on using the repository browser.