| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // 081024 G4NucleiPropertiesTable:: to G4Nucle
27 //
28 #include "G4QMDGroundStateNucleus.hh"
29
30 #include "G4NucleiProperties.hh"
31 #include "G4Proton.hh"
32 #include "G4Neutron.hh"
33 #include "G4Exp.hh"
34 #include "G4Pow.hh"
35 #include "G4PhysicalConstants.hh"
36 #include "Randomize.hh"
37
38 G4QMDGroundStateNucleus::G4QMDGroundStateNucle
39 : maxTrial ( 1000 )
40 , r00 ( 1.124 ) // radius parameter for Woods
41 , r01 ( 0.5 ) // radius parameter for Woods
42 , saa ( 0.2 ) // diffuse parameter for init
43 , rada ( 0.9 ) // cutoff parameter
44 , radb ( 0.3 ) // cutoff parameter
45 , dsam ( 1.5 ) // minimum distance for same
46 , ddif ( 1.0 ) // minimum distance for diffe
47 , edepth ( 0.0 )
48 , epse ( 0.000001 ) // torelance for energy i
49 , meanfield ( NULL )
50 {
51
52 //std::cout << " G4QMDGroundStateNucleus( G
53
54 dsam2 = dsam*dsam;
55 ddif2 = ddif*ddif;
56
57 G4QMDParameters* parameters = G4QMDParamete
58
59 hbc = parameters->Get_hbc();
60 gamm = parameters->Get_gamm();
61 cpw = parameters->Get_cpw();
62 cph = parameters->Get_cph();
63 epsx = parameters->Get_epsx();
64 cpc = parameters->Get_cpc();
65
66 cdp = parameters->Get_cdp();
67 c0p = parameters->Get_c0p();
68 c3p = parameters->Get_c3p();
69 csp = parameters->Get_csp();
70 clp = parameters->Get_clp();
71
72 ebini = 0.0;
73 // Following 10 lines should be here, right
74 // Otherwise, mass number cannot be conserv
75 // the target are nucleons.
76 //Nucleon primary or target case;
77 if ( z == 1 && a == 1 ) { // Hydrogen Cas
78 SetParticipant( new G4QMDParticipant( G4
79 // ebini = 0.0;
80 return;
81 }
82 else if ( z == 0 && a == 1 ) { // Neutron p
83 SetParticipant( new G4QMDParticipant( G4
84 // ebini = 0.0;
85 return;
86 }
87
88
89 //edepth = 0.0;
90
91 for ( G4int i = 0 ; i < a ; ++i )
92 {
93
94 G4ParticleDefinition* pd;
95
96 if ( i < z )
97 {
98 pd = G4Proton::Proton();
99 }
100 else
101 {
102 pd = G4Neutron::Neutron();
103 }
104
105 G4ThreeVector p( 0.0 );
106 G4ThreeVector r( 0.0 );
107 G4QMDParticipant* aParticipant = new G4Q
108 SetParticipant( aParticipant );
109
110 }
111
112 G4double radious = r00 * G4Pow::GetInstance
113
114 rt00 = radious - r01;
115 radm = radious - rada * ( gamm - 1.0 ) + ra
116 rmax = 1.0 / ( 1.0 + G4Exp ( -rt00/saa ) );
117
118 //maxTrial = 1000;
119
120
121 meanfield = new G4QMDMeanField();
122 meanfield->SetSystem( this );
123
124 //std::cout << "G4QMDGroundStateNucleus( G4
125 packNucleons();
126 //std::cout << "G4QMDGroundStateNucleus( G4
127
128 delete meanfield;
129
130 }
131
132
133
134 void G4QMDGroundStateNucleus::packNucleons()
135 {
136
137 //std::cout << "G4QMDGroundStateNucleus::pa
138
139 ebini = - G4NucleiProperties::GetBindingEne
140
141 G4double ebin00 = ebini * 0.001;
142
143 G4double ebin0 = 0.0;
144 G4double ebin1 = 0.0;
145
146 if ( GetMassNumber() != 4 )
147 {
148 ebin0 = ( ebini - 0.5 ) * 0.001;
149 ebin1 = ( ebini + 0.5 ) * 0.001;
150 }
151 else
152 {
153 ebin0 = ( ebini - 1.5 ) * 0.001;
154 ebin1 = ( ebini + 1.5 ) * 0.001;
155 }
156
157 G4int n0Try = 0;
158 G4bool isThisOK = false;
159 while ( n0Try < maxTrial ) // Loop checking
160 {
161 n0Try++;
162 //std::cout << "TKDB packNucleons n0Try
163
164 // Sampling Position
165
166 //std::cout << "TKDB Sampling Position "
167
168 G4bool areThesePsOK = false;
169 G4int npTry = 0;
170 while ( npTry < maxTrial ) // Loop check
171 {
172 //std::cout << "TKDB Sampling Position n
173 npTry++;
174 G4int i = 0;
175 if ( samplingPosition( i ) )
176 {
177 //std::cout << "packNucleons sampl
178 for ( i = 1 ; i < GetMassNumber()
179 {
180 //std::cout << "packNucleons sa
181 if ( !( samplingPosition( i ) )
182 {
183 //std::cout << "packNucleons
184 break;
185 }
186 }
187 if ( i == GetMassNumber() )
188 {
189 //std::cout << "packNucleons sa
190 areThesePsOK = true;
191 break;
192 }
193 }
194 }
195 if ( areThesePsOK == false ) continue;
196
197 //std::cout << "TKDB Sampling Position E
198
199 // Calculate Two-body quantities
200
201 meanfield->Cal2BodyQuantities();
202 std::vector< G4double > rho_a ( GetMassN
203 std::vector< G4double > rho_s ( GetMassN
204 std::vector< G4double > rho_c ( GetMassN
205
206 rho_l.resize ( GetMassNumber() , 0.0 );
207 d_pot.resize ( GetMassNumber() , 0.0 );
208
209 for ( G4int i = 0 ; i < GetMassNumber()
210 {
211 for ( G4int j = 0 ; j < GetMassNumber
212 {
213
214 rho_a[ i ] += meanfield->GetRHA( j
215 G4int k = 0;
216 if ( participants[j]->GetDefinitio
217 {
218 k = 1;
219 }
220
221 rho_s[ i ] += meanfield->GetRHA( j
222
223 rho_c[ i ] += meanfield->GetRHE( j
224 }
225
226 }
227
228 for ( G4int i = 0 ; i < GetMassNumber()
229 {
230 rho_l[i] = cdp * ( rho_a[i] + 1.0 );
231 d_pot[i] = c0p * rho_a[i]
232 + c3p * G4Pow::GetInstance()
233 + csp * rho_s[i]
234 + clp * rho_c[i];
235
236 //std::cout << "d_pot[i] " << i << "
237 }
238
239
240 // Sampling Momentum
241
242 //std::cout << "TKDB Sampling Momentum "
243
244 phase_g.clear();
245 phase_g.resize( GetMassNumber() , 0.0 );
246
247 //std::cout << "TKDB Sampling Momentum 1
248 G4bool isThis1stMOK = false;
249 G4int nmTry = 0;
250 while ( nmTry < maxTrial ) // Loop check
251 {
252 nmTry++;
253 G4int i = 0;
254 if ( samplingMomentum( i ) )
255 {
256 isThis1stMOK = true;
257 break;
258 }
259 }
260 if ( isThis1stMOK == false ) continue;
261
262 //std::cout << "TKDB Sampling Momentum 2
263
264 G4bool areTheseMsOK = true;
265 nmTry = 0;
266 while ( nmTry < maxTrial ) // Loop check
267 {
268 nmTry++;
269 G4int i = 0;
270 for ( i = 1 ; i < GetMassNumber()
271 {
272 //std::cout << "TKDB packNucleo
273 if ( !( samplingMomentum( i ) )
274 {
275 //std::cout << "TKDB packNuc
276 areTheseMsOK = false;
277 break;
278 }
279 }
280 if ( i == GetMassNumber() )
281 {
282 areTheseMsOK = true;
283 }
284
285 break;
286 }
287 if ( areTheseMsOK == false ) continue;
288
289 // Kill Angluar Momentum
290
291 //std::cout << "TKDB Sampling Kill Anglu
292
293 killCMMotionAndAngularM();
294
295
296 // Check Binding Energy
297
298 //std::cout << "packNucleons Check Bindi
299
300 G4double ekinal = 0.0;
301 for ( int i = 0 ; i < GetMassNumber() ;
302 {
303 ekinal += participants[i]->GetKinetic
304 }
305
306 meanfield->Cal2BodyQuantities();
307
308 G4double totalPotentialE = meanfield->Ge
309 G4double ebinal = ( totalPotentialE + ek
310
311 //std::cout << "packNucleons totalPotent
312 //std::cout << "packNucleons ebinal " <<
313 //std::cout << "packNucleons ekinal " <<
314
315 if ( ebinal < ebin0 || ebinal > ebin1 )
316 {
317 //std::cout << "packNucleons ebin0 "
318 //std::cout << "packNucleons ebin1 "
319 //std::cout << "packNucleons ebinal "
320 //std::cout << "packNucleons Check Bindi
321 continue;
322 }
323
324 //std::cout << "packNucleons Check Bindi
325
326
327 // Energy Adujstment
328
329 G4double dtc = 1.0;
330 G4double frg = -0.1;
331 G4double rdf0 = 0.5;
332
333 G4double edif0 = ebinal - ebin00;
334
335 G4double cfrc = 0.0;
336 if ( 0 < edif0 )
337 {
338 cfrc = frg;
339 }
340 else
341 {
342 cfrc = -frg;
343 }
344
345 G4int ifrc = 1;
346
347 G4int neaTry = 0;
348
349 G4bool isThisEAOK = false;
350 while ( neaTry < maxTrial ) // Loop che
351 {
352 neaTry++;
353
354 G4double edif = ebinal - ebin00;
355
356 //std::cout << "TKDB edif " << edif <
357 if ( std::abs ( edif ) < epse )
358 {
359
360 isThisEAOK = true;
361 //std::cout << "isThisEAOK " << is
362 break;
363 }
364
365 G4int jfrc = 0;
366 if ( edif < 0.0 )
367 {
368 jfrc = 1;
369 }
370 else
371 {
372 jfrc = -1;
373 }
374
375 if ( jfrc != ifrc )
376 {
377 cfrc = -rdf0 * cfrc;
378 dtc = rdf0 * dtc;
379 }
380
381 if ( jfrc == ifrc && std::abs( edif0
382 {
383 cfrc = -rdf0 * cfrc;
384 dtc = rdf0 * dtc;
385 }
386
387 ifrc = jfrc;
388 edif0 = edif;
389
390 meanfield->CalGraduate();
391
392 for ( int i = 0 ; i < GetMassNumber()
393 {
394 G4ThreeVector ri = participants[i]
395 G4ThreeVector p_i = participants[i
396
397 ri += dtc * ( meanfield->GetFFr(i)
398 p_i += dtc * ( meanfield->GetFFp(i
399
400 participants[i]->SetPosition( ri )
401 participants[i]->SetMomentum( p_i
402 }
403
404 ekinal = 0.0;
405
406 for ( int i = 0 ; i < GetMassNumber()
407 {
408 ekinal += participants[i]->GetKine
409 }
410
411 meanfield->Cal2BodyQuantities();
412 totalPotentialE = meanfield->GetTotal
413
414 ebinal = ( totalPotentialE + ekinal )
415
416 }
417 //std::cout << "isThisEAOK " << isThisEA
418 if ( isThisEAOK == false ) continue;
419
420 isThisOK = true;
421 //std::cout << "isThisOK " << isThisOK <
422 break;
423
424 }
425
426 if ( isThisOK == false )
427 {
428 G4cout << "GroundStateNucleus state cann
429 }
430
431 //std::cout << "packNucleons End" << std::e
432 return;
433 }
434
435
436 G4bool G4QMDGroundStateNucleus::samplingPositi
437 {
438
439 G4bool result = false;
440
441 G4int nTry = 0;
442
443 while ( nTry < maxTrial ) // Loop checking
444 {
445
446 //std::cout << "samplingPosition i th pa
447
448 G4double rwod = -1.0;
449 G4double rrr = 0.0;
450
451 G4double rx = 0.0;
452 G4double ry = 0.0;
453 G4double rz = 0.0;
454
455 G4int icounter = 0;
456 G4int icounter_max = 1024;
457 while ( G4UniformRand() * rmax > rwod )
458 {
459 icounter++;
460 if ( icounter > icounter_max ) {
461 G4cout << "Loop-counter exceeded the thr
462 break;
463 }
464
465 G4double rsqr = 10.0;
466 G4int jcounter = 0;
467 G4int jcounter_max = 1024;
468 while ( rsqr > 1.0 ) // Loop checking
469 {
470 jcounter++;
471 if ( jcounter > jcounter_max ) {
472 G4cout << "Loop-counter exceeded the
473 break;
474 }
475 rx = 1.0 - 2.0 * G4UniformRand();
476 ry = 1.0 - 2.0 * G4UniformRand();
477 rz = 1.0 - 2.0 * G4UniformRand();
478 rsqr = rx*rx + ry*ry + rz*rz;
479 }
480 rrr = radm * std::sqrt ( rsqr );
481 rwod = 1.0 / ( 1.0 + G4Exp ( ( rrr -
482
483 }
484
485 participants[i]->SetPosition( G4ThreeVec
486
487 if ( i == 0 )
488 {
489 result = true;
490 return result;
491 }
492
493 // i > 1 ( Second Particle or later )
494 // Check Distance to others
495
496 G4bool isThisOK = true;
497 for ( G4int j = 0 ; j < i ; j++ )
498 {
499
500 G4double r2 = participants[j]->GetPo
501 G4double dmin2 = 0.0;
502
503 if ( participants[j]->GetDefinition()
504 {
505 dmin2 = dsam2;
506 }
507 else
508 {
509 dmin2 = ddif2;
510 }
511
512 //std::cout << "distance between j an
513 if ( r2 < dmin2 )
514 {
515 isThisOK = false;
516 break;
517 }
518
519 }
520
521 if ( isThisOK == true )
522 {
523 result = true;
524 return result;
525 }
526
527 nTry++;
528
529 }
530
531 // Here return "false"
532 return result;
533 }
534
535
536
537 G4bool G4QMDGroundStateNucleus::samplingMoment
538 {
539
540 //std::cout << "TKDB samplingMomentum for "
541
542 G4bool result = false;
543
544 G4double pfm = hbc * G4Pow::GetInstance()->
545
546 if ( 10 < GetMassNumber() && -5.5 < ebini
547 {
548 pfm = pfm * ( 1.0 + 0.2 * std::sqrt( std
549 }
550
551 //std::cout << "TKDB samplingMomentum pfm "
552
553 std::vector< G4double > phase;
554 phase.resize( i+1 ); // i start from 0
555
556 G4int ntry = 0;
557 // 710
558 while ( ntry < maxTrial ) // Loop checking
559 {
560
561 //std::cout << " TKDB ntry " << ntry <<
562 ntry++;
563
564 G4double ke = DBL_MAX;
565
566 G4int tkdb_i =0;
567 // 700
568 G4int icounter = 0;
569 G4int icounter_max = 1024;
570 while ( ke + d_pot [i] > edepth ) // Loo
571 {
572 icounter++;
573 if ( icounter > icounter_max ) {
574 G4cout << "Loop-counter exceeded the thr
575 break;
576 }
577
578 G4double psqr = 10.0;
579 G4double px = 0.0;
580 G4double py = 0.0;
581 G4double pz = 0.0;
582
583 G4int jcounter = 0;
584 G4int jcounter_max = 1024;
585 while ( psqr > 1.0 ) // Loop checking
586 {
587 jcounter++;
588 if ( jcounter > jcounter_max ) {
589 G4cout << "Loop-counter exceeded the
590 break;
591 }
592 px = 1.0 - 2.0*G4UniformRand();
593 py = 1.0 - 2.0*G4UniformRand();
594 pz = 1.0 - 2.0*G4UniformRand();
595
596 psqr = px*px + py*py + pz*pz;
597 }
598
599 G4ThreeVector p ( px , py , pz );
600 p = pfm * p;
601 participants[i]->SetMomentum( p );
602 G4LorentzVector p4 = participants[i]-
603 //ke = p4.e() - p4.restMass();
604 ke = participants[i]->GetKineticEnerg
605
606
607 tkdb_i++;
608 if ( tkdb_i > maxTrial ) return resul
609
610 }
611
612 //std::cout << "TKDB ke d_pot[i] " << ke
613
614
615 if ( i == 0 )
616 {
617 result = true;
618 return result;
619 }
620
621 G4bool isThisOK = true;
622
623 // Check Pauli principle
624
625 phase[ i ] = 0.0;
626
627 //std::cout << "TKDB Check Pauli Princip
628
629 for ( G4int j = 0 ; j < i ; j++ )
630 {
631 phase[ j ] = 0.0;
632 //std::cout << "TKDB Check Pauli Prin
633 G4double expa = 0.0;
634 if ( participants[j]->GetDefinition()
635 {
636
637 expa = - meanfield->GetRR2(i,j) *
638
639 if ( expa > epsx )
640 {
641 G4ThreeVector p_i = participant
642 G4ThreeVector pj = participants
643 G4double dist2_p = p_i.diff2( p
644
645 dist2_p = dist2_p*cph;
646 expa = expa - dist2_p;
647
648 if ( expa > epsx )
649 {
650
651 phase[j] = G4Exp ( expa );
652
653 if ( phase[j] * cpc > 0.2 )
654 {
655 /*
656 G4cout << "TKDB Check Pauli Principle
657 G4cout << "TKDB Check Pauli Principle
658 G4cout << "TKDB Check Pauli Principle
659 */
660 isThisOK = false;
661 break;
662 }
663 if ( ( phase_g[j] + phase[j]
664 {
665 /*
666 G4cout << "TKDB Check Pauli Principle
667 G4cout << "TKDB Check Pauli Principle
668 G4cout << "TKDB Check Pauli Principle
669 G4cout << "TKDB Check Pauli Principle
670 */
671 isThisOK = false;
672 break;
673 }
674
675 phase[i] += phase[j];
676 if ( phase[i] * cpc > 0.3 )
677 {
678 /*
679 G4cout << "TKDB Check Pauli Principle
680 G4cout << "TKDB Check Pauli Principle
681 G4cout << "TKDB Check Pauli Principle
682 */
683 isThisOK = false;
684 break;
685 }
686
687 //std::cout << "TKDB Check P
688
689 }
690 else
691 {
692 //std::cout << "TKDB Check P
693 }
694
695 }
696 else
697 {
698 //std::cout << "TKDB Check Paul
699 }
700
701 }
702 else
703 {
704 //std::cout << "TKDB Check Pauli P
705 }
706
707 }
708
709 if ( isThisOK == true )
710 {
711
712 phase_g[i] = phase[i];
713
714 for ( int j = 0 ; j < i ; j++ )
715 {
716 phase_g[j] += phase[j];
717 }
718
719 result = true;
720 return result;
721 }
722
723 }
724
725 return result;
726
727 }
728
729
730
731 void G4QMDGroundStateNucleus::killCMMotionAndA
732 {
733
734 // CalEnergyAndAngularMomentumInCM();
735
736 //std::vector< G4ThreeVector > p ( GetMassN
737 //std::vector< G4ThreeVector > r ( GetMassN
738
739 // Move to cm system
740
741 G4ThreeVector pcm_tmp ( 0.0 );
742 G4ThreeVector rcm_tmp ( 0.0 );
743 G4double sumMass = 0.0;
744
745 for ( G4int i = 0 ; i < GetMassNumber() ; i
746 {
747 pcm_tmp += participants[i]->GetMomentum(
748 rcm_tmp += participants[i]->GetPosition(
749 sumMass += participants[i]->GetMass();
750 }
751
752 pcm_tmp = pcm_tmp/GetMassNumber();
753 rcm_tmp = rcm_tmp/sumMass;
754
755 for ( G4int i = 0 ; i < GetMassNumber() ; i
756 {
757 participants[i]->SetMomentum( participan
758 participants[i]->SetPosition( participan
759 }
760
761 // kill the angular momentum
762
763 G4ThreeVector ll ( 0.0 );
764 for ( G4int i = 0 ; i < GetMassNumber() ; i
765 {
766 ll += participants[i]->GetPosition().cro
767 }
768
769 G4double rr[3][3];
770 G4double ss[3][3];
771 for ( G4int i = 0 ; i < 3 ; i++ )
772 {
773 for ( G4int j = 0 ; j < 3 ; j++ )
774 {
775 rr [i][j] = 0.0;
776
777 if ( i == j )
778 {
779 ss [i][j] = 1.0;
780 }
781 else
782 {
783 ss [i][j] = 0.0;
784 }
785 }
786 }
787
788 for ( G4int i = 0 ; i < GetMassNumber() ; i
789 {
790 G4ThreeVector r = participants[i]->GetPo
791 rr[0][0] += r.y() * r.y() + r.z() * r.z(
792 rr[1][0] += - r.y() * r.x();
793 rr[2][0] += - r.z() * r.x();
794 rr[0][1] += - r.x() * r.y();
795 rr[1][1] += r.z() * r.z() + r.x() * r.x(
796 rr[2][1] += - r.z() * r.y();
797 rr[2][0] += - r.x() * r.z();
798 rr[2][1] += - r.y() * r.z();
799 rr[2][2] += r.x() * r.x() + r.y() * r.y(
800 }
801
802 for ( G4int i = 0 ; i < 3 ; i++ )
803 {
804 G4double x = rr [i][i];
805 for ( G4int j = 0 ; j < 3 ; j++ )
806 {
807 rr[i][j] = rr[i][j] / x;
808 ss[i][j] = ss[i][j] / x;
809 }
810 for ( G4int j = 0 ; j < 3 ; j++ )
811 {
812 if ( j != i )
813 {
814 G4double y = rr [j][i];
815 for ( G4int k = 0 ; k < 3 ; k++ )
816 {
817 rr[j][k] += -y * rr[i][k];
818 ss[j][k] += -y * ss[i][k];
819 }
820 }
821 }
822 }
823
824 G4double opl[3];
825 G4double rll[3];
826
827 rll[0] = ll.x();
828 rll[1] = ll.y();
829 rll[2] = ll.z();
830
831 for ( G4int i = 0 ; i < 3 ; i++ )
832 {
833 opl[i] = 0.0;
834
835 for ( G4int j = 0 ; j < 3 ; j++ )
836 {
837 opl[i] += ss[i][j]*rll[j];
838 }
839 }
840
841 for ( G4int i = 0 ; i < GetMassNumber() ; i
842 {
843 G4ThreeVector p_i = participants[i]->Get
844 G4ThreeVector ri = participants[i]->GetP
845 G4ThreeVector opl_v ( opl[0] , opl[1] ,
846
847 p_i += ri.cross(opl_v);
848
849 participants[i]->SetMomentum( p_i );
850 }
851
852 }
853