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

Last change on this file since 1036 was 962, checked in by garnier, 17 years ago

update processes
File size: 26.2 KB
Line 
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
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
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
118void 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
648G4bool 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
735G4bool 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
915void 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}
Note: See TracBrowser for help on using the repository browser.