00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "G4QMDGroundStateNucleus.hh"
00029
00030 #include "G4NucleiProperties.hh"
00031
00032 #include "G4Proton.hh"
00033 #include "G4Neutron.hh"
00034
00035 #include "G4PhysicalConstants.hh"
00036 #include "Randomize.hh"
00037
00038 G4QMDGroundStateNucleus::G4QMDGroundStateNucleus( G4int z , G4int a )
00039 : r00 ( 1.124 )
00040 , r01 ( 0.5 )
00041 , saa ( 0.2 )
00042 , rada ( 0.9 )
00043 , radb ( 0.3 )
00044 , dsam ( 1.5 )
00045 , ddif ( 1.0 )
00046 , epse ( 0.000001 )
00047 {
00048
00049
00050
00051 if ( z == 1 && a == 1 )
00052 {
00053 SetParticipant( new G4QMDParticipant( G4Proton::Proton() , G4ThreeVector( 0.0 ) , G4ThreeVector( 0.0 ) ) );
00054 return;
00055 }
00056 else if ( z == 0 && a == 1 )
00057 {
00058 SetParticipant( new G4QMDParticipant( G4Neutron::Neutron() , G4ThreeVector( 0.0 ) , G4ThreeVector( 0.0 ) ) );
00059 return;
00060 }
00061
00062 dsam2 = dsam*dsam;
00063 ddif2 = ddif*ddif;
00064
00065 G4QMDParameters* parameters = G4QMDParameters::GetInstance();
00066
00067 hbc = parameters->Get_hbc();
00068 gamm = parameters->Get_gamm();
00069 cpw = parameters->Get_cpw();
00070 cph = parameters->Get_cph();
00071 epsx = parameters->Get_epsx();
00072 cpc = parameters->Get_cpc();
00073
00074 cdp = parameters->Get_cdp();
00075 c0p = parameters->Get_c0p();
00076 c3p = parameters->Get_c3p();
00077 csp = parameters->Get_csp();
00078 clp = parameters->Get_clp();
00079
00080 edepth = 0.0;
00081
00082 for ( int i = 0 ; i < a ; i++ )
00083 {
00084
00085 G4ParticleDefinition* pd;
00086
00087 if ( i < z )
00088 {
00089 pd = G4Proton::Proton();
00090 }
00091 else
00092 {
00093 pd = G4Neutron::Neutron();
00094 }
00095
00096 G4ThreeVector p( 0.0 );
00097 G4ThreeVector r( 0.0 );
00098 G4QMDParticipant* aParticipant = new G4QMDParticipant( pd , p , r );
00099 SetParticipant( aParticipant );
00100
00101 }
00102
00103 G4double radious = r00 * std::pow ( double ( GetMassNumber() ) , 1.0/3.0 );
00104
00105 rt00 = radious - r01;
00106 radm = radious - rada * ( gamm - 1.0 ) + radb;
00107 rmax = 1.0 / ( 1.0 + std::exp ( -rt00/saa ) );
00108
00109 maxTrial = 1000;
00110 meanfield = new G4QMDMeanField();
00111 meanfield->SetSystem( this );
00112
00113
00114 packNucleons();
00115
00116
00117 delete meanfield;
00118
00119 }
00120
00121
00122
00123 void G4QMDGroundStateNucleus::packNucleons()
00124 {
00125
00126
00127
00128 ebini = - G4NucleiProperties::GetBindingEnergy( GetMassNumber() , GetAtomicNumber() ) / GetMassNumber();
00129
00130 G4double ebin00 = ebini * 0.001;
00131
00132 G4double ebin0 = 0.0;
00133 G4double ebin1 = 0.0;
00134
00135 if ( GetMassNumber() != 4 )
00136 {
00137 ebin0 = ( ebini - 0.5 ) * 0.001;
00138 ebin1 = ( ebini + 0.5 ) * 0.001;
00139 }
00140 else
00141 {
00142 ebin0 = ( ebini - 1.5 ) * 0.001;
00143 ebin1 = ( ebini + 1.5 ) * 0.001;
00144 }
00145
00146 {
00147 G4int n0Try = 0;
00148 G4bool isThisOK = false;
00149 while ( n0Try < maxTrial )
00150 {
00151 n0Try++;
00152
00153
00154
00155
00156
00157
00158 G4bool areThesePsOK = false;
00159 G4int npTry = 0;
00160 while ( npTry < maxTrial )
00161 {
00162
00163 npTry++;
00164 G4int i = 0;
00165 if ( samplingPosition( i ) )
00166 {
00167
00168 for ( i = 1 ; i < GetMassNumber() ; i++ )
00169 {
00170
00171 if ( !( samplingPosition( i ) ) )
00172 {
00173
00174 break;
00175 }
00176 }
00177 if ( i == GetMassNumber() )
00178 {
00179
00180 areThesePsOK = true;
00181 break;
00182 }
00183 }
00184 }
00185 if ( areThesePsOK == false ) continue;
00186
00187
00188
00189
00190
00191 meanfield->Cal2BodyQuantities();
00192 std::vector< G4double > rho_a ( GetMassNumber() , 0.0 );
00193 std::vector< G4double > rho_s ( GetMassNumber() , 0.0 );
00194 std::vector< G4double > rho_c ( GetMassNumber() , 0.0 );
00195
00196 rho_l.resize ( GetMassNumber() , 0.0 );
00197 d_pot.resize ( GetMassNumber() , 0.0 );
00198
00199 for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
00200 {
00201 for ( G4int j = 0 ; j < GetMassNumber() ; j++ )
00202 {
00203
00204 rho_a[ i ] += meanfield->GetRHA( j , i );
00205 G4int k = 0;
00206 if ( participants[j]->GetDefinition() != participants[i]->GetDefinition() )
00207 {
00208 k = 1;
00209 }
00210
00211 rho_s[ i ] += meanfield->GetRHA( j , i )*( 1.0 - 2.0 * k );
00212
00213 rho_c[ i ] += meanfield->GetRHE( j , i );
00214 }
00215
00216 }
00217
00218 for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
00219 {
00220 rho_l[i] = cdp * ( rho_a[i] + 1.0 );
00221 d_pot[i] = c0p * rho_a[i]
00222 + c3p * std::pow ( rho_a[i] , gamm )
00223 + csp * rho_s[i]
00224 + clp * rho_c[i];
00225
00226
00227 }
00228
00229
00230
00231
00232
00233
00234 phase_g.clear();
00235 phase_g.resize( GetMassNumber() , 0.0 );
00236
00237
00238 G4bool isThis1stMOK = false;
00239 G4int nmTry = 0;
00240 while ( nmTry < maxTrial )
00241 {
00242 nmTry++;
00243 G4int i = 0;
00244 if ( samplingMomentum( i ) )
00245 {
00246 isThis1stMOK = true;
00247 break;
00248 }
00249 }
00250 if ( isThis1stMOK == false ) continue;
00251
00252
00253
00254 G4bool areTheseMsOK = true;
00255 nmTry = 0;
00256 while ( nmTry < maxTrial )
00257 {
00258 nmTry++;
00259 G4int i = 0;
00260 for ( i = 1 ; i < GetMassNumber() ; i++ )
00261 {
00262
00263 if ( !( samplingMomentum( i ) ) )
00264 {
00265
00266 areTheseMsOK = false;
00267 break;
00268 }
00269 }
00270 if ( i == GetMassNumber() )
00271 {
00272 areTheseMsOK = true;
00273 }
00274
00275 break;
00276 }
00277 if ( areTheseMsOK == false ) continue;
00278
00279
00280
00281
00282
00283 killCMMotionAndAngularM();
00284
00285
00286
00287
00288
00289
00290 G4double ekinal = 0.0;
00291 for ( int i = 0 ; i < GetMassNumber() ; i++ )
00292 {
00293 ekinal += participants[i]->GetKineticEnergy();
00294 }
00295
00296 meanfield->Cal2BodyQuantities();
00297
00298 G4double totalPotentialE = meanfield->GetTotalPotential();
00299 G4double ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() );
00300
00301
00302
00303
00304
00305 if ( ebinal < ebin0 || ebinal > ebin1 )
00306 {
00307
00308
00309
00310
00311 continue;
00312 }
00313
00314
00315
00316
00317
00318
00319 G4double dtc = 1.0;
00320 G4double frg = -0.1;
00321 G4double rdf0 = 0.5;
00322
00323 G4double edif0 = ebinal - ebin00;
00324
00325 G4double cfrc = 0.0;
00326 if ( 0 < edif0 )
00327 {
00328 cfrc = frg;
00329 }
00330 else
00331 {
00332 cfrc = -frg;
00333 }
00334
00335 G4int ifrc = 1;
00336
00337 G4int neaTry = 0;
00338
00339 G4bool isThisEAOK = false;
00340 while ( neaTry < maxTrial )
00341 {
00342 neaTry++;
00343
00344 G4double edif = ebinal - ebin00;
00345
00346
00347 if ( std::abs ( edif ) < epse )
00348 {
00349
00350 isThisEAOK = true;
00351
00352 break;
00353 }
00354
00355 G4int jfrc = 0;
00356 if ( edif < 0.0 )
00357 {
00358 jfrc = 1;
00359 }
00360 else
00361 {
00362 jfrc = -1;
00363 }
00364
00365 if ( jfrc != ifrc )
00366 {
00367 cfrc = -rdf0 * cfrc;
00368 dtc = rdf0 * dtc;
00369 }
00370
00371 if ( jfrc == ifrc && std::abs( edif0 ) < std::abs( edif ) )
00372 {
00373 cfrc = -rdf0 * cfrc;
00374 dtc = rdf0 * dtc;
00375 }
00376
00377 ifrc = jfrc;
00378 edif0 = edif;
00379
00380 meanfield->CalGraduate();
00381
00382 for ( int i = 0 ; i < GetMassNumber() ; i++ )
00383 {
00384 G4ThreeVector ri = participants[i]->GetPosition();
00385 G4ThreeVector p_i = participants[i]->GetMomentum();
00386
00387 ri += dtc * ( meanfield->GetFFr(i) - cfrc * ( meanfield->GetFFp(i) ) );
00388 p_i += dtc * ( meanfield->GetFFp(i) + cfrc * ( meanfield->GetFFr(i) ) );
00389
00390 participants[i]->SetPosition( ri );
00391 participants[i]->SetMomentum( p_i );
00392 }
00393
00394 ekinal = 0.0;
00395
00396 for ( int i = 0 ; i < GetMassNumber() ; i++ )
00397 {
00398 ekinal += participants[i]->GetKineticEnergy();
00399 }
00400
00401 meanfield->Cal2BodyQuantities();
00402 totalPotentialE = meanfield->GetTotalPotential();
00403
00404 ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() );
00405
00406 }
00407
00408 if ( isThisEAOK == false ) continue;
00409
00410 isThisOK = true;
00411
00412 break;
00413
00414 }
00415
00416 if ( isThisOK == false )
00417 {
00418 std::cout << "GroundStateNucleus state cannot be created. Try again with another parameters." << std::endl;
00419 }
00420
00421
00422 return;
00423
00424 }
00425
00426
00427
00428
00429
00430
00431
00432
00433 G4int n0Try = 0;
00434
00435 while ( n0Try < maxTrial )
00436 {
00437 if ( samplingPosition( 0 ) )
00438 {
00439 G4int i = 0;
00440 for ( i = 1 ; i < GetMassNumber() ; i++ )
00441 {
00442 if ( !( samplingPosition( i ) ) )
00443 {
00444 break;
00445 }
00446 }
00447 if ( i == GetMassNumber() ) break;
00448 }
00449 n0Try++;
00450 }
00451
00452 if ( n0Try > maxTrial )
00453 {
00454 std::cout << "GroundStateNucleus state cannot be created. Try again with another parameters." << std::endl;
00455 return;
00456 }
00457
00458 meanfield->Cal2BodyQuantities();
00459 std::vector< G4double > rho_a ( GetMassNumber() , 0.0 );
00460 std::vector< G4double > rho_s ( GetMassNumber() , 0.0 );
00461 std::vector< G4double > rho_c ( GetMassNumber() , 0.0 );
00462
00463 rho_l.resize ( GetMassNumber() , 0.0 );
00464 d_pot.resize ( GetMassNumber() , 0.0 );
00465
00466 for ( int i = 0 ; i < GetMassNumber() ; i++ )
00467 {
00468 for ( int j = 0 ; j < GetMassNumber() ; j++ )
00469 {
00470
00471 rho_a[ i ] += meanfield->GetRHA( j , i );
00472 G4int k = 0;
00473 if ( participants[i]->GetDefinition() != participants[j]->GetDefinition() )
00474 {
00475 k = 1;
00476 }
00477
00478 rho_s[ i ] += meanfield->GetRHA( j , i )*( 1.0 - 2.0 * k );
00479
00480 rho_c[ i ] += meanfield->GetRHE( j , i );
00481 }
00482 }
00483
00484 for ( int i = 0 ; i < GetMassNumber() ; i++ )
00485 {
00486 rho_l[i] = cdp * ( rho_a[i] + 1.0 );
00487 d_pot[i] = c0p * rho_a[i]
00488 + c3p * std::pow ( rho_a[i] , gamm )
00489 + csp * rho_s[i]
00490 + clp * rho_c[i];
00491 }
00492
00493
00494
00495
00496
00497
00498
00499
00500 phase_g.resize( GetMassNumber() , 0.0 );
00501
00502
00503 samplingMomentum( 0 );
00504
00505 G4int n1Try = 0;
00506
00507
00508 while ( n1Try < maxTrial )
00509 {
00510 if ( samplingPosition( 0 ) )
00511 {
00512 G4int i = 0;
00513 G4bool isThisOK = true;
00514 for ( i = 1 ; i < GetMassNumber() ; i++ )
00515 {
00516 if ( !( samplingMomentum( i ) ) )
00517 {
00518 isThisOK = false;
00519 break;
00520 }
00521 }
00522 if ( isThisOK == true ) break;
00523
00524 }
00525 n1Try++;
00526 }
00527
00528 if ( n1Try > maxTrial )
00529 {
00530 std::cout << "GroundStateNucleus state cannot be created. Try again with another parameters." << std::endl;
00531 return;
00532 }
00533
00534
00535
00536
00537 killCMMotionAndAngularM();
00538
00539
00540
00541 G4double ekinal = 0.0;
00542 for ( int i = 0 ; i < GetMassNumber() ; i++ )
00543 {
00544 ekinal += participants[i]->GetKineticEnergy();
00545 }
00546
00547 meanfield->Cal2BodyQuantities();
00548 G4double totalPotentialE = meanfield->GetTotalPotential();
00549 G4double ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() );
00550
00551 if ( ebinal < ebin0 || ebinal > ebin1 )
00552 {
00553
00554 }
00555
00556
00557
00558
00559 G4double dtc = 1.0;
00560 G4double frg = -0.1;
00561 G4double rdf0 = 0.5;
00562
00563 G4double edif0 = ebinal - ebin00;
00564
00565 G4double cfrc = 0.0;
00566 if ( 0 < edif0 )
00567 {
00568 cfrc = frg;
00569 }
00570 else
00571 {
00572 cfrc = -frg;
00573 }
00574
00575 G4int ifrc = 1;
00576
00577 G4int ntryACH = 0;
00578
00579 G4bool isThisOK = false;
00580 while ( ntryACH < maxTrial )
00581 {
00582
00583 G4double edif = ebinal - ebin00;
00584 if ( std::abs ( edif ) < epse )
00585 {
00586 isThisOK = true;
00587 break;
00588 }
00589
00590 G4int jfrc = 0;
00591 if ( edif < 0.0 )
00592 {
00593 jfrc = 1;
00594 }
00595 else
00596 {
00597 jfrc = -1;
00598 }
00599
00600 if ( jfrc != ifrc )
00601 {
00602 cfrc = -rdf0 * cfrc;
00603 dtc = rdf0 * dtc;
00604 }
00605
00606 if ( jfrc == ifrc && std::abs( edif0 ) < std::abs( edif ) )
00607 {
00608 cfrc = -rdf0 * cfrc;
00609 dtc = rdf0 * dtc;
00610 }
00611
00612 ifrc = jfrc;
00613 edif0 = edif;
00614
00615 meanfield->CalGraduate();
00616
00617 for ( int i = 0 ; i < GetMassNumber() ; i++ )
00618 {
00619 G4ThreeVector ri = participants[i]->GetPosition();
00620 G4ThreeVector p_i = participants[i]->GetMomentum();
00621
00622 ri += dtc * ( meanfield->GetFFr(i) - cfrc * ( meanfield->GetFFp(i) ) );
00623 p_i += dtc * ( meanfield->GetFFp(i) + cfrc * ( meanfield->GetFFr(i) ) );
00624
00625 participants[i]->SetPosition( ri );
00626 participants[i]->SetMomentum( p_i );
00627 }
00628
00629 ekinal = 0.0;
00630
00631 for ( int i = 0 ; i < GetMassNumber() ; i++ )
00632 {
00633 ekinal += participants[i]->GetKineticEnergy();
00634 }
00635
00636 meanfield->Cal2BodyQuantities();
00637 totalPotentialE = meanfield->GetTotalPotential();
00638
00639 ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() );
00640
00641
00642 ntryACH++;
00643 }
00644
00645 if ( isThisOK )
00646 {
00647 return;
00648 }
00649
00650 }
00651
00652
00653 G4bool G4QMDGroundStateNucleus::samplingPosition( G4int i )
00654 {
00655
00656 G4bool result = false;
00657
00658 G4int nTry = 0;
00659
00660 while ( nTry < maxTrial )
00661 {
00662
00663
00664
00665 G4double rwod = -1.0;
00666 G4double rrr = 0.0;
00667
00668 G4double rx = 0.0;
00669 G4double ry = 0.0;
00670 G4double rz = 0.0;
00671
00672 while ( G4UniformRand() * rmax > rwod )
00673 {
00674
00675 G4double rsqr = 10.0;
00676 while ( rsqr > 1.0 )
00677 {
00678 rx = 1.0 - 2.0 * G4UniformRand();
00679 ry = 1.0 - 2.0 * G4UniformRand();
00680 rz = 1.0 - 2.0 * G4UniformRand();
00681 rsqr = rx*rx + ry*ry + rz*rz;
00682 }
00683 rrr = radm * std::sqrt ( rsqr );
00684 rwod = 1.0 / ( 1.0 + std::exp ( ( rrr - rt00 ) / saa ) );
00685
00686 }
00687
00688 participants[i]->SetPosition( G4ThreeVector( rx , ry , rz )*radm );
00689
00690 if ( i == 0 )
00691 {
00692 result = true;
00693 return result;
00694 }
00695
00696
00697
00698
00699 G4bool isThisOK = true;
00700 for ( G4int j = 0 ; j < i ; j++ )
00701 {
00702
00703 G4double r2 = participants[j]->GetPosition().diff2( participants[i]->GetPosition() );
00704 G4double dmin2 = 0.0;
00705
00706 if ( participants[j]->GetDefinition() == participants[i]->GetDefinition() )
00707 {
00708 dmin2 = dsam2;
00709 }
00710 else
00711 {
00712 dmin2 = ddif2;
00713 }
00714
00715
00716 if ( r2 < dmin2 )
00717 {
00718 isThisOK = false;
00719 break;
00720 }
00721
00722 }
00723
00724 if ( isThisOK == true )
00725 {
00726 result = true;
00727 return result;
00728 }
00729
00730 nTry++;
00731
00732 }
00733
00734
00735 return result;
00736 }
00737
00738
00739
00740 G4bool G4QMDGroundStateNucleus::samplingMomentum( G4int i )
00741 {
00742
00743
00744
00745 G4bool result = false;
00746
00747 G4double pfm = hbc * std::pow ( ( 3.0 / 2.0 * pi*pi * rho_l[i] ) , 1./3. );
00748
00749 if ( 10 < GetMassNumber() && -5.5 < ebini )
00750 {
00751 pfm = pfm * ( 1.0 + 0.2 * std::sqrt( std::abs( 8.0 + ebini ) / 8.0 ) );
00752 }
00753
00754
00755
00756 std::vector< G4double > phase;
00757 phase.resize( i+1 );
00758
00759 G4int ntry = 0;
00760
00761 while ( ntry < maxTrial )
00762 {
00763
00764
00765 ntry++;
00766
00767 G4double ke = DBL_MAX;
00768
00769 G4int tkdb_i =0;
00770
00771 while ( ke + d_pot [i] > edepth )
00772 {
00773
00774 G4double psqr = 10.0;
00775 G4double px = 0.0;
00776 G4double py = 0.0;
00777 G4double pz = 0.0;
00778
00779 while ( psqr > 1.0 )
00780 {
00781 px = 1.0 - 2.0*G4UniformRand();
00782 py = 1.0 - 2.0*G4UniformRand();
00783 pz = 1.0 - 2.0*G4UniformRand();
00784
00785 psqr = px*px + py*py + pz*pz;
00786 }
00787
00788 G4ThreeVector p ( px , py , pz );
00789 p = pfm * p;
00790 participants[i]->SetMomentum( p );
00791 G4LorentzVector p4 = participants[i]->Get4Momentum();
00792
00793 ke = participants[i]->GetKineticEnergy();
00794
00795
00796 tkdb_i++;
00797 if ( tkdb_i > maxTrial ) return result;
00798
00799 }
00800
00801
00802
00803
00804 if ( i == 0 )
00805 {
00806 result = true;
00807 return result;
00808 }
00809
00810 G4bool isThisOK = true;
00811
00812
00813
00814 phase[ i ] = 0.0;
00815
00816
00817
00818 for ( G4int j = 0 ; j < i ; j++ )
00819 {
00820 phase[ j ] = 0.0;
00821
00822 G4double expa = 0.0;
00823 if ( participants[j]->GetDefinition() == participants[i]->GetDefinition() )
00824 {
00825
00826 expa = - meanfield->GetRR2(i,j) * cpw;
00827
00828 if ( expa > epsx )
00829 {
00830 G4ThreeVector p_i = participants[i]->GetMomentum();
00831 G4ThreeVector pj = participants[j]->GetMomentum();
00832 G4double dist2_p = p_i.diff2( pj );
00833
00834 dist2_p = dist2_p*cph;
00835 expa = expa - dist2_p;
00836
00837 if ( expa > epsx )
00838 {
00839
00840 phase[j] = std::exp ( expa );
00841
00842 if ( phase[j] * cpc > 0.2 )
00843 {
00844
00845
00846
00847
00848
00849 isThisOK = false;
00850 break;
00851 }
00852 if ( ( phase_g[j] + phase[j] ) * cpc > 0.5 )
00853 {
00854
00855
00856
00857
00858
00859
00860 isThisOK = false;
00861 break;
00862 }
00863
00864 phase[i] += phase[j];
00865 if ( phase[i] * cpc > 0.3 )
00866 {
00867
00868
00869
00870
00871
00872 isThisOK = false;
00873 break;
00874 }
00875
00876
00877
00878 }
00879 else
00880 {
00881
00882 }
00883
00884 }
00885 else
00886 {
00887
00888 }
00889
00890 }
00891 else
00892 {
00893
00894 }
00895
00896 }
00897
00898 if ( isThisOK == true )
00899 {
00900
00901 phase_g[i] = phase[i];
00902
00903 for ( int j = 0 ; j < i ; j++ )
00904 {
00905 phase_g[j] += phase[j];
00906 }
00907
00908 result = true;
00909 return result;
00910 }
00911
00912 }
00913
00914 return result;
00915
00916 }
00917
00918
00919
00920 void G4QMDGroundStateNucleus::killCMMotionAndAngularM()
00921 {
00922
00923
00924
00925
00926
00927
00928
00929
00930 G4ThreeVector pcm_tmp ( 0.0 );
00931 G4ThreeVector rcm_tmp ( 0.0 );
00932 G4double sumMass = 0.0;
00933
00934 for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
00935 {
00936 pcm_tmp += participants[i]->GetMomentum();
00937 rcm_tmp += participants[i]->GetPosition() * participants[i]->GetMass();
00938 sumMass += participants[i]->GetMass();
00939 }
00940
00941 pcm_tmp = pcm_tmp/GetMassNumber();
00942 rcm_tmp = rcm_tmp/sumMass;
00943
00944 for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
00945 {
00946 participants[i]->SetMomentum( participants[i]->GetMomentum() - pcm_tmp );
00947 participants[i]->SetPosition( participants[i]->GetPosition() - rcm_tmp );
00948 }
00949
00950
00951
00952 G4ThreeVector ll ( 0.0 );
00953 for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
00954 {
00955 ll += participants[i]->GetPosition().cross ( participants[i]->GetMomentum() );
00956 }
00957
00958 G4double rr[3][3];
00959 G4double ss[3][3];
00960 for ( G4int i = 0 ; i < 3 ; i++ )
00961 {
00962 for ( G4int j = 0 ; j < 3 ; j++ )
00963 {
00964 rr [i][j] = 0.0;
00965
00966 if ( i == j )
00967 {
00968 ss [i][j] = 1.0;
00969 }
00970 else
00971 {
00972 ss [i][j] = 0.0;
00973 }
00974 }
00975 }
00976
00977 for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
00978 {
00979 G4ThreeVector r = participants[i]->GetPosition();
00980 rr[0][0] += r.y() * r.y() + r.z() * r.z();
00981 rr[1][0] += - r.y() * r.x();
00982 rr[2][0] += - r.z() * r.x();
00983 rr[0][1] += - r.x() * r.y();
00984 rr[1][1] += r.z() * r.z() + r.x() * r.x();
00985 rr[2][1] += - r.z() * r.y();
00986 rr[2][0] += - r.x() * r.z();
00987 rr[2][1] += - r.y() * r.z();
00988 rr[2][2] += r.x() * r.x() + r.y() * r.y();
00989 }
00990
00991 for ( G4int i = 0 ; i < 3 ; i++ )
00992 {
00993 G4double x = rr [i][i];
00994 for ( G4int j = 0 ; j < 3 ; j++ )
00995 {
00996 rr[i][j] = rr[i][j] / x;
00997 ss[i][j] = ss[i][j] / x;
00998 }
00999 for ( G4int j = 0 ; j < 3 ; j++ )
01000 {
01001 if ( j != i )
01002 {
01003 G4double y = rr [j][i];
01004 for ( G4int k = 0 ; k < 3 ; k++ )
01005 {
01006 rr[j][k] += -y * rr[i][k];
01007 ss[j][k] += -y * ss[i][k];
01008 }
01009 }
01010 }
01011 }
01012
01013 G4double opl[3];
01014 G4double rll[3];
01015
01016 rll[0] = ll.x();
01017 rll[1] = ll.y();
01018 rll[2] = ll.z();
01019
01020 for ( G4int i = 0 ; i < 3 ; i++ )
01021 {
01022 opl[i] = 0.0;
01023
01024 for ( G4int j = 0 ; j < 3 ; j++ )
01025 {
01026 opl[i] += ss[i][j]*rll[j];
01027 }
01028 }
01029
01030 for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
01031 {
01032 G4ThreeVector p_i = participants[i]->GetMomentum() ;
01033 G4ThreeVector ri = participants[i]->GetPosition() ;
01034 G4ThreeVector opl_v ( opl[0] , opl[1] , opl[2] );
01035
01036 p_i += ri.cross(opl_v);
01037
01038 participants[i]->SetMomentum( p_i );
01039 }
01040
01041 }