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

Last change on this file since 830 was 819, checked in by garnier, 17 years ago

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