#include <G4QMDMeanField.hh>
Public Member Functions | |
G4QMDMeanField () | |
~G4QMDMeanField () | |
void | SetSystem (G4QMDSystem *aSystem) |
void | SetNucleus (G4QMDNucleus *aSystem) |
G4QMDSystem * | GetSystem () |
void | Cal2BodyQuantities () |
void | Cal2BodyQuantities (G4int) |
void | CalGraduate () |
G4bool | IsPauliBlocked (G4int) |
G4double | GetTotalPotential () |
G4double | GetPotential (G4int) |
void | DoPropagation (G4double) |
std::vector< G4QMDNucleus * > | DoClusterJudgment () |
G4double | GetRR2 (G4int i, G4int j) |
G4double | GetRHA (G4int i, G4int j) |
G4double | GetRHE (G4int i, G4int j) |
G4ThreeVector | GetFFr (G4int i) |
G4ThreeVector | GetFFp (G4int i) |
std::vector< G4double > | GetLocalDensity () |
std::vector< G4double > | GetDepthOfPotential () |
void | Update () |
Definition at line 44 of file G4QMDMeanField.hh.
G4QMDMeanField::G4QMDMeanField | ( | ) |
Definition at line 41 of file G4QMDMeanField.cc.
References G4QMDParameters::Get_c0(), G4QMDParameters::Get_c3(), G4QMDParameters::Get_cl(), G4QMDParameters::Get_cpc(), G4QMDParameters::Get_cph(), G4QMDParameters::Get_cpw(), G4QMDParameters::Get_cs(), G4QMDParameters::Get_gamm(), G4QMDParameters::Get_hbc(), G4QMDParameters::Get_rho0(), G4QMDParameters::Get_wl(), G4QMDParameters::GetInstance(), and G4INCL::Math::pi.
00042 : rclds ( 4.0 ) // distance for cluster judgement 00043 , epsx ( -20.0 ) // gauss term 00044 , epscl ( 0.0001 ) // coulomb term 00045 , irelcr ( 1 ) 00046 { 00047 00048 G4QMDParameters* parameters = G4QMDParameters::GetInstance(); 00049 wl = parameters->Get_wl(); 00050 cl = parameters->Get_cl(); 00051 rho0 = parameters->Get_rho0(); 00052 hbc = parameters->Get_hbc(); 00053 gamm = parameters->Get_gamm(); 00054 00055 cpw = parameters->Get_cpw(); 00056 cph = parameters->Get_cph(); 00057 cpc = parameters->Get_cpc(); 00058 00059 c0 = parameters->Get_c0(); 00060 c3 = parameters->Get_c3(); 00061 cs = parameters->Get_cs(); 00062 00063 // distance 00064 c0w = 1.0/4.0/wl; 00065 //c3w = 1.0/4.0/wl; //no need 00066 c0sw = std::sqrt( c0w ); 00067 clw = 2.0 / std::sqrt ( 4.0 * pi * wl ); 00068 00069 // graduate 00070 c0g = - c0 / ( 2.0 * wl ); 00071 c3g = - c3 / ( 4.0 * wl ) * gamm; 00072 csg = - cs / ( 2.0 * wl ); 00073 pag = gamm - 1; 00074 00075 }
G4QMDMeanField::~G4QMDMeanField | ( | ) |
void G4QMDMeanField::Cal2BodyQuantities | ( | G4int | ) |
Definition at line 248 of file G4QMDMeanField.cc.
References G4QMDParticipant::Get4Momentum(), G4QMDParticipant::GetBaryonNumber(), G4QMDParticipant::GetChargeInUnitOfEplus(), G4QMDSystem::GetParticipant(), G4QMDParticipant::GetPosition(), and G4QMDSystem::GetTotalNumberOfParticipant().
00249 { 00250 00251 //std::cout << "Cal2BodyQuantities " << i << std::endl; 00252 00253 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 00254 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum(); 00255 00256 00257 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; j ++ ) 00258 { 00259 if ( j == i ) continue; 00260 00261 G4ThreeVector rj = system->GetParticipant( j )->GetPosition(); 00262 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum(); 00263 00264 G4ThreeVector rij = ri - rj; 00265 G4ThreeVector pij = (p4i - p4j).v(); 00266 G4LorentzVector p4ij = p4i - p4j; 00267 G4ThreeVector bij = ( p4i + p4j ).boostVector(); 00268 G4double gammaij = ( p4i + p4j ).gamma(); 00269 00270 G4double eij = ( p4i + p4j ).e(); 00271 00272 G4double rbrb = rij*bij; 00273 // G4double bij2 = bij*bij; 00274 G4double rij2 = rij*rij; 00275 G4double pij2 = pij*pij; 00276 00277 rbrb = irelcr * rbrb; 00278 G4double gamma2_ij = gammaij*gammaij; 00279 00280 /* 00281 G4double rbrb = 0.0; 00282 G4double beta2_ij = 0.0; 00283 G4double rij2 = 0.0; 00284 G4double pij2 = 0.0; 00285 00286 // 00287 G4LorentzVector p4ip4j = p4i + p4j; 00288 G4double eij = p4ip4j.e(); 00289 00290 G4ThreeVector r = ri - rj; 00291 G4LorentzVector p4 = p4i - p4j; 00292 00293 rbrb = r.x()*p4ip4j.x()/eij 00294 + r.y()*p4ip4j.y()/eij 00295 + r.z()*p4ip4j.z()/eij; 00296 00297 beta2_ij = ( p4ip4j.x()*p4ip4j.x() + p4ip4j.y()*p4ip4j.y() + p4ip4j.z()*p4ip4j.z() ) / ( eij*eij ); 00298 rij2 = r*r; 00299 pij2 = p4.v()*p4.v(); 00300 00301 rbrb = irelcr * rbrb; 00302 00303 G4double gamma2_ij = 1 / ( 1 - beta2_ij ); 00304 */ 00305 00306 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb; 00307 rr2[j][i] = rr2[i][j]; 00308 00309 rbij[i][j] = gamma2_ij * rbrb; 00310 rbij[j][i] = - rbij[i][j]; 00311 00312 pp2[i][j] = pij2 00313 + irelcr * ( - std::pow ( p4i.e() - p4j.e() , 2 ) 00314 + gamma2_ij * std::pow ( ( ( p4i.m2() - p4j.m2() ) / eij ) , 2 ) ); 00315 00316 pp2[j][i] = pp2[i][j]; 00317 00318 // Gauss term 00319 00320 G4double expa1 = - rr2[i][j] * c0w; 00321 00322 G4double rh1; 00323 if ( expa1 > epsx ) 00324 { 00325 rh1 = std::exp( expa1 ); 00326 } 00327 else 00328 { 00329 rh1 = 0.0; 00330 } 00331 00332 G4int ibry = system->GetParticipant(i)->GetBaryonNumber(); 00333 G4int jbry = system->GetParticipant(j)->GetBaryonNumber(); 00334 00335 00336 rha[i][j] = ibry*jbry*rh1; 00337 rha[j][i] = rha[i][j]; 00338 00339 // Coulomb terms 00340 00341 G4double rrs2 = rr2[i][j] + epscl; 00342 G4double rrs = std::sqrt ( rrs2 ); 00343 00344 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 00345 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 00346 00347 G4double erf = 0.0; 00348 // T. K. add this protection. 5.8 is good enough for double 00349 if ( rrs*c0sw < 5.8 ) 00350 erf = CLHEP::HepStat::erf ( rrs*c0sw ); 00351 else 00352 erf = 1.0; 00353 00354 G4double erfij = erf/rrs; 00355 00356 00357 rhe[i][j] = icharge*jcharge * erfij; 00358 00359 rhe[j][i] = rhe[i][j]; 00360 00361 // G4double clw; 00362 00363 rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2; 00364 00365 rhc[j][i] = rhc[i][j]; 00366 00367 } 00368 00369 }
void G4QMDMeanField::Cal2BodyQuantities | ( | ) |
Definition at line 148 of file G4QMDMeanField.cc.
References G4QMDParticipant::Get4Momentum(), G4QMDParticipant::GetBaryonNumber(), G4QMDParticipant::GetChargeInUnitOfEplus(), G4QMDSystem::GetParticipant(), G4QMDParticipant::GetPosition(), and G4QMDSystem::GetTotalNumberOfParticipant().
Referenced by G4QMDCollision::CalFinalStateOfTheBinaryCollision(), G4QMDCollision::CalFinalStateOfTheBinaryCollisionJQMD(), G4QMDCollision::CalKinematicsOfBinaryCollisions(), DoClusterJudgment(), DoPropagation(), and SetSystem().
00149 { 00150 00151 if ( system->GetTotalNumberOfParticipant() < 2 ) return; 00152 00153 for ( G4int j = 1 ; j < system->GetTotalNumberOfParticipant() ; j++ ) 00154 { 00155 00156 G4ThreeVector rj = system->GetParticipant( j )->GetPosition(); 00157 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum(); 00158 00159 for ( G4int i = 0 ; i < j ; i++ ) 00160 { 00161 00162 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 00163 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum(); 00164 00165 G4ThreeVector rij = ri - rj; 00166 G4ThreeVector pij = (p4i - p4j).v(); 00167 G4LorentzVector p4ij = p4i - p4j; 00168 G4ThreeVector bij = ( p4i + p4j ).boostVector(); 00169 G4double gammaij = ( p4i + p4j ).gamma(); 00170 00171 G4double eij = ( p4i + p4j ).e(); 00172 00173 G4double rbrb = rij*bij; 00174 // G4double bij2 = bij*bij; 00175 G4double rij2 = rij*rij; 00176 G4double pij2 = pij*pij; 00177 00178 rbrb = irelcr * rbrb; 00179 G4double gamma2_ij = gammaij*gammaij; 00180 00181 00182 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb; 00183 rr2[j][i] = rr2[i][j]; 00184 00185 rbij[i][j] = gamma2_ij * rbrb; 00186 rbij[j][i] = - rbij[i][j]; 00187 00188 pp2[i][j] = pij2 00189 + irelcr * ( - std::pow ( p4i.e() - p4j.e() , 2 ) 00190 + gamma2_ij * std::pow ( ( ( p4i.m2() - p4j.m2() ) / eij ) , 2 ) ); 00191 00192 00193 pp2[j][i] = pp2[i][j]; 00194 00195 // Gauss term 00196 00197 G4double expa1 = - rr2[i][j] * c0w; 00198 00199 G4double rh1; 00200 if ( expa1 > epsx ) 00201 { 00202 rh1 = std::exp( expa1 ); 00203 } 00204 else 00205 { 00206 rh1 = 0.0; 00207 } 00208 00209 G4int ibry = system->GetParticipant(i)->GetBaryonNumber(); 00210 G4int jbry = system->GetParticipant(j)->GetBaryonNumber(); 00211 00212 00213 rha[i][j] = ibry*jbry*rh1; 00214 rha[j][i] = rha[i][j]; 00215 00216 // Coulomb terms 00217 00218 G4double rrs2 = rr2[i][j] + epscl; 00219 G4double rrs = std::sqrt ( rrs2 ); 00220 00221 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 00222 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 00223 00224 G4double erf = 0.0; 00225 // T. K. add this protection. 5.8 is good enough for double 00226 if ( rrs*c0sw < 5.8 ) 00227 erf = CLHEP::HepStat::erf ( rrs*c0sw ); 00228 else 00229 erf = 1.0; 00230 00231 G4double erfij = erf/rrs; 00232 00233 00234 rhe[i][j] = icharge*jcharge * erfij; 00235 00236 rhe[j][i] = rhe[i][j]; 00237 00238 rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2; 00239 00240 rhc[j][i] = rhc[i][j]; 00241 00242 } // i 00243 } // j 00244 }
void G4QMDMeanField::CalGraduate | ( | ) |
Definition at line 373 of file G4QMDMeanField.cc.
References G4QMDParticipant::Get4Momentum(), G4QMDParticipant::GetChargeInUnitOfEplus(), G4QMDParticipant::GetNuc(), G4QMDSystem::GetParticipant(), G4QMDParticipant::GetPosition(), GetPotential(), and G4QMDSystem::GetTotalNumberOfParticipant().
Referenced by DoPropagation().
00374 { 00375 00376 ffr.resize( system->GetTotalNumberOfParticipant() ); 00377 ffp.resize( system->GetTotalNumberOfParticipant() ); 00378 rh3d.resize( system->GetTotalNumberOfParticipant() ); 00379 00380 for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; i ++ ) 00381 { 00382 G4double rho3 = 0.0; 00383 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; j ++ ) 00384 { 00385 rho3 += rha[j][i]; 00386 } 00387 rh3d[i] = std::pow ( rho3 , pag ); 00388 } 00389 00390 00391 for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; i ++ ) 00392 { 00393 00394 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 00395 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum(); 00396 00397 G4ThreeVector betai = p4i.v()/p4i.e(); 00398 00399 // R-JQMD 00400 G4double Vi = GetPotential( i ); 00401 G4double p_zero = std::sqrt( p4i.e()*p4i.e() + 2*p4i.m()*Vi); 00402 G4ThreeVector betai_R = p4i.v()/p_zero; 00403 G4double mi_R = p4i.m()/p_zero; 00404 // 00405 ffr[i] = betai_R; 00406 ffp[i] = G4ThreeVector( 0.0 ); 00407 00408 if ( false ) 00409 { 00410 ffr[i] = betai; 00411 mi_R = 1.0; 00412 } 00413 00414 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; j ++ ) 00415 { 00416 00417 G4ThreeVector rj = system->GetParticipant( j )->GetPosition(); 00418 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum(); 00419 00420 G4double eij = p4i.e() + p4j.e(); 00421 00422 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 00423 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 00424 00425 G4int inuc = system->GetParticipant(i)->GetNuc(); 00426 G4int jnuc = system->GetParticipant(j)->GetNuc(); 00427 00428 G4double ccpp = c0g * rha[j][i] 00429 + c3g * rha[j][i] * ( rh3d[j] + rh3d[i] ) 00430 + csg * rha[j][i] * jnuc * inuc 00431 * ( 1. - 2. * std::abs( jcharge - icharge ) ) 00432 + cl * rhc[j][i]; 00433 ccpp *= mi_R; 00434 00435 /* 00436 std::cout << c0g << " " << c3g << " " << csg << " " << cl << std::endl; 00437 std::cout << "ccpp " << i << " " << j << " " << ccpp << std::endl; 00438 std::cout << "rha[j][i] " << rha[j][i] << std::endl; 00439 std::cout << "rh3d " << rh3d[j] << " " << rh3d[i] << std::endl; 00440 std::cout << "rhc[j][i] " << rhc[j][i] << std::endl; 00441 */ 00442 00443 G4double grbb = - rbij[j][i]; 00444 G4double ccrr = grbb * ccpp / eij; 00445 00446 /* 00447 std::cout << "ccrr " << ccrr << std::endl; 00448 std::cout << "grbb " << grbb << std::endl; 00449 */ 00450 00451 00452 G4ThreeVector rij = ri - rj; 00453 G4ThreeVector betaij = ( p4i + p4j ).v()/eij; 00454 00455 G4ThreeVector cij = betaij - betai; 00456 00457 ffr[i] = ffr[i] + 2*ccrr* ( rij + grbb*cij ); 00458 00459 ffp[i] = ffp[i] - 2*ccpp* ( rij + grbb*betaij ); 00460 00461 } 00462 } 00463 00464 //std::cout << "gradu 0 " << ffr[0] << " " << ffp[0] << std::endl; 00465 //std::cout << "gradu 1 " << ffr[1] << " " << ffp[1] << std::endl; 00466 00467 }
std::vector< G4QMDNucleus * > G4QMDMeanField::DoClusterJudgment | ( | ) |
Definition at line 680 of file G4QMDMeanField.cc.
References Cal2BodyQuantities(), G4QMDParticipant::GetBaryonNumber(), G4QMDSystem::GetParticipant(), G4QMDSystem::GetTotalNumberOfParticipant(), CLHEP::detail::n, G4INCL::Math::pi, G4QMDSystem::SetParticipant(), and G4QMDSystem::SubtractSystem().
Referenced by G4QMDReaction::ApplyYourself().
00681 { 00682 00683 //std::cout << "MeanField DoClusterJudgemnt" << std::endl; 00684 00685 Cal2BodyQuantities(); 00686 00687 G4double cpf2 = std::pow ( 1.5 * pi*pi * std::pow ( 4.0 * pi * wl , -1.5 ) 00688 , 00689 2./3. ) 00690 * hbc * hbc; 00691 G4double rcc2 = rclds*rclds; 00692 00693 G4int n = system->GetTotalNumberOfParticipant(); 00694 std::vector < G4double > rhoa; 00695 rhoa.resize ( n ); 00696 00697 for ( G4int i = 0 ; i < n ; i++ ) 00698 { 00699 rhoa[i] = 0.0; 00700 00701 if ( system->GetParticipant( i )->GetBaryonNumber() == 1 ) 00702 { 00703 for ( G4int j = 0 ; j < n ; j++ ) 00704 { 00705 if ( system->GetParticipant( j )->GetBaryonNumber() == 1 ) 00706 rhoa[i] += rha[i][j]; 00707 } 00708 } 00709 00710 rhoa[i] = std::pow ( rhoa[i] + 1 , 1.0/3.0 ); 00711 00712 } 00713 00714 // identification of the cluster 00715 00716 std::map < G4int , std::vector < G4int > > cluster_map; 00717 std::vector < G4bool > is_already_belong_some_cluster; 00718 00719 // cluster_id participant_id 00720 std::multimap < G4int , G4int > comb_map; 00721 std::multimap < G4int , G4int > assign_map; 00722 assign_map.clear(); 00723 00724 std::vector < G4int > mascl; 00725 std::vector < G4int > num; 00726 mascl.resize ( n ); 00727 num.resize ( n ); 00728 is_already_belong_some_cluster.resize ( n ); 00729 00730 std::vector < G4int > is_assigned_to ( n , -1 ); 00731 std::multimap < G4int , G4int > clusters; 00732 00733 for ( G4int i = 0 ; i < n ; i++ ) 00734 { 00735 mascl[i] = 1; 00736 num[i] = 1; 00737 00738 is_already_belong_some_cluster[i] = false; 00739 } 00740 00741 00742 G4int nclst = 1; 00743 G4int ichek = 1; 00744 00745 00746 G4int id = 0; 00747 G4int cluster_id = -1; 00748 for ( G4int i = 0 ; i < n-1 ; i++ ) 00749 { 00750 00751 G4bool hasThisCompany = false; 00752 // Check only for bryons? 00753 // std::cout << "Check Baryon " << i << std::endl; 00754 00755 if ( system->GetParticipant( i )->GetBaryonNumber() == 1 ) 00756 { 00757 00758 // if ( is_already_belong_some_cluster[i] != true ) 00759 // { 00760 //G4int j1 = ichek + 1; 00761 G4int j1 = i + 1; 00762 for ( G4int j = j1 ; j < n ; j++ ) 00763 { 00764 00765 std::vector < G4int > cluster_participants; 00766 if ( system->GetParticipant( j )->GetBaryonNumber() == 1 ) 00767 { 00768 G4double rdist2 = rr2[ i ][ j ]; 00769 G4double pdist2 = pp2[ i ][ j ]; 00770 //G4double rdist2 = rr2[ num[i] ][ num[j] ]; 00771 //G4double pdist2 = pp2[ num[i] ][ num[j] ]; 00772 G4double pcc2 = cpf2 00773 * ( rhoa[ i ] + rhoa[ j ] ) 00774 * ( rhoa[ i ] + rhoa[ j ] ); 00775 00776 // Check phase space: close enough? 00777 if ( rdist2 < rcc2 && pdist2 < pcc2 ) 00778 { 00779 00780 /* 00781 std::cout << "G4QMDRESULT " 00782 << i << " " << j << " " << id << " " 00783 << is_assigned_to [ i ] << " " << is_assigned_to [ j ] 00784 << std::endl; 00785 */ 00786 00787 if ( is_assigned_to [ j ] == -1 ) 00788 { 00789 if ( is_assigned_to [ i ] == -1 ) 00790 { 00791 if ( clusters.size() != 0 ) 00792 { 00793 id = clusters.rbegin()->first + 1; 00794 //std::cout << "id is increare " << id << std::endl; 00795 } 00796 else 00797 { 00798 id = 0; 00799 } 00800 clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , i ) ); 00801 is_assigned_to [ i ] = id; 00802 clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , j ) ); 00803 is_assigned_to [ j ] = id; 00804 } 00805 else 00806 { 00807 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ i ] , j ) ); 00808 is_assigned_to [ j ] = is_assigned_to [ i ]; 00809 } 00810 } 00811 else 00812 { 00813 // j is already belong to some cluester 00814 if ( is_assigned_to [ i ] == -1 ) 00815 { 00816 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , i ) ); 00817 is_assigned_to [ i ] = is_assigned_to [ j ]; 00818 } 00819 else 00820 { 00821 // i has companion 00822 if ( is_assigned_to [ i ] != is_assigned_to [ j ] ) 00823 { 00824 // move companions to the cluster 00825 // 00826 //std::cout << "combine " << is_assigned_to [ i ] << " to " << is_assigned_to [ j ] << std::endl; 00827 std::multimap< G4int , G4int > clusters_tmp; 00828 G4int target_cluster_id; 00829 if ( is_assigned_to [ i ] > is_assigned_to [ j ] ) 00830 target_cluster_id = is_assigned_to [ i ]; 00831 else 00832 target_cluster_id = is_assigned_to [ j ]; 00833 00834 for ( std::multimap< G4int , G4int >::iterator it 00835 = clusters.begin() ; it != clusters.end() ; it++ ) 00836 { 00837 00838 //std::cout << it->first << " " << it->second << " " << target_cluster_id << std::endl; 00839 if ( it->first == target_cluster_id ) 00840 { 00841 //std::cout << "move " << it->first << " " << it->second << std::endl; 00842 is_assigned_to [ it->second ] = is_assigned_to [ j ]; 00843 clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , it->second ) ); 00844 } 00845 else 00846 { 00847 clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( it->first , it->second ) ); 00848 } 00849 } 00850 00851 clusters = clusters_tmp; 00852 //id = clusters.rbegin()->first; 00853 //id = target_cluster_id; 00854 //std::cout << "id " << id << std::endl; 00855 } 00856 } 00857 } 00858 00859 //std::cout << "combination " << i << " " << j << std::endl; 00860 comb_map.insert( std::multimap<G4int,G4int>::value_type ( i , j ) ); 00861 cluster_participants.push_back ( j ); 00862 00863 00864 00865 if ( assign_map.find( cluster_id ) == assign_map.end() ) 00866 { 00867 is_already_belong_some_cluster[i] = true; 00868 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , i ) ); 00869 hasThisCompany = true; 00870 } 00871 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , j ) ); 00872 is_already_belong_some_cluster[j] = true; 00873 00874 } 00875 00876 if ( ichek == i ) 00877 { 00878 nclst++; 00879 ichek++; 00880 } 00881 } 00882 00883 if ( cluster_participants.size() > 0 ) 00884 { 00885 // cluster , participant 00886 cluster_map.insert ( std::pair < G4int , std::vector < G4int > > ( i , cluster_participants ) ); 00887 } 00888 } 00889 // } 00890 } 00891 if ( hasThisCompany == true ) cluster_id++; 00892 } 00893 00894 //std::cout << " id " << id << std::endl; 00895 00896 // sort 00897 // Heavy cluster comes first 00898 // size cluster_id 00899 std::multimap< G4int , G4int > sorted_cluster_map; 00900 for ( G4int i = 0 ; i <= id ; i++ ) // << "<=" because id is highest cluster nubmer. 00901 { 00902 00903 //std::cout << i << " cluster has " << clusters.count( i ) << " nucleons." << std::endl; 00904 sorted_cluster_map.insert ( std::multimap<G4int,G4int>::value_type ( (G4int) clusters.count( i ) , i ) ); 00905 00906 } 00907 00908 00909 // create nucleus from devided clusters 00910 std::vector < G4QMDNucleus* > result; 00911 for ( std::multimap < G4int , G4int >::reverse_iterator it 00912 = sorted_cluster_map.rbegin() ; it != sorted_cluster_map.rend() ; it ++) 00913 { 00914 00915 //G4cout << "Add Participants to cluseter " << it->second << G4endl; 00916 00917 if ( it->first != 0 ) 00918 { 00919 G4QMDNucleus* nucleus = new G4QMDNucleus(); 00920 for ( std::multimap < G4int , G4int >::iterator itt 00921 = clusters.begin() ; itt != clusters.end() ; itt ++) 00922 { 00923 00924 if ( it->second == itt->first ) 00925 { 00926 nucleus->SetParticipant( system->GetParticipant ( itt->second ) ); 00927 //G4cout << "Add Participants " << itt->second << " " << system->GetParticipant ( itt->second )->GetPosition() << G4endl; 00928 } 00929 00930 } 00931 result.push_back( nucleus ); 00932 } 00933 00934 } 00935 00936 // delete participants from current system 00937 00938 for ( std::vector < G4QMDNucleus* > ::iterator it 00939 = result.begin() ; it != result.end() ; it++ ) 00940 { 00941 system->SubtractSystem ( *it ); 00942 } 00943 00944 return result; 00945 00946 }
void G4QMDMeanField::DoPropagation | ( | G4double | ) |
Definition at line 618 of file G4QMDMeanField.cc.
References Cal2BodyQuantities(), CalGraduate(), G4QMDParticipant::GetMomentum(), G4QMDSystem::GetParticipant(), G4QMDParticipant::GetPosition(), G4QMDSystem::GetTotalNumberOfParticipant(), CLHEP::detail::n, G4QMDParticipant::SetMomentum(), and G4QMDParticipant::SetPosition().
Referenced by G4QMDReaction::ApplyYourself().
00619 { 00620 00621 G4double cc2 = 1.0; 00622 G4double cc1 = 1.0 - cc2; 00623 G4double cc3 = 1.0 / 2.0 / cc2; 00624 00625 G4double dt3 = dt * cc3; 00626 G4double dt1 = dt * ( cc1 - cc3 ); 00627 G4double dt2 = dt * cc2; 00628 00629 CalGraduate(); 00630 00631 G4int n = system->GetTotalNumberOfParticipant(); 00632 00633 // 1st Step 00634 00635 std::vector< G4ThreeVector > f0r, f0p; 00636 f0r.resize( n ); 00637 f0p.resize( n ); 00638 00639 for ( G4int i = 0 ; i < n ; i++ ) 00640 { 00641 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 00642 G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum(); 00643 00644 ri += dt3* ffr[i]; 00645 p3i += dt3* ffp[i]; 00646 00647 f0r[i] = ffr[i]; 00648 f0p[i] = ffp[i]; 00649 00650 system->GetParticipant( i )->SetPosition( ri ); 00651 system->GetParticipant( i )->SetMomentum( p3i ); 00652 00653 // we do not need set total momentum by ourselvs 00654 } 00655 00656 // 2nd Step 00657 Cal2BodyQuantities(); 00658 CalGraduate(); 00659 00660 for ( G4int i = 0 ; i < n ; i++ ) 00661 { 00662 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 00663 G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum(); 00664 00665 ri += dt1* f0r[i] + dt2* ffr[i]; 00666 p3i += dt1* f0p[i] + dt2* ffp[i]; 00667 00668 system->GetParticipant( i )->SetPosition( ri ); 00669 system->GetParticipant( i )->SetMomentum( p3i ); 00670 00671 // we do not need set total momentum by ourselvs 00672 } 00673 00674 Cal2BodyQuantities(); 00675 00676 }
std::vector< G4double > G4QMDMeanField::GetDepthOfPotential | ( | ) |
G4ThreeVector G4QMDMeanField::GetFFp | ( | G4int | i | ) | [inline] |
G4ThreeVector G4QMDMeanField::GetFFr | ( | G4int | i | ) | [inline] |
std::vector< G4double > G4QMDMeanField::GetLocalDensity | ( | ) |
Definition at line 471 of file G4QMDMeanField.cc.
References G4QMDParticipant::GetChargeInUnitOfEplus(), G4QMDParticipant::GetNuc(), G4QMDSystem::GetParticipant(), G4QMDSystem::GetTotalNumberOfParticipant(), and CLHEP::detail::n.
Referenced by CalGraduate().
00472 { 00473 G4int n = system->GetTotalNumberOfParticipant(); 00474 00475 G4double rhoa = 0.0; 00476 G4double rho3 = 0.0; 00477 G4double rhos = 0.0; 00478 G4double rhoc = 0.0; 00479 00480 00481 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 00482 G4int inuc = system->GetParticipant(i)->GetNuc(); 00483 00484 for ( G4int j = 0 ; j < n ; j ++ ) 00485 { 00486 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 00487 G4int jnuc = system->GetParticipant(j)->GetNuc(); 00488 00489 rhoa += rha[j][i]; 00490 rhoc += rhe[j][i]; 00491 rhos += rha[j][i] * jnuc * inuc 00492 * ( 1 - 2 * std::abs ( jcharge - icharge ) ); 00493 } 00494 00495 rho3 = std::pow ( rhoa , gamm ); 00496 00497 G4double potential = c0 * rhoa 00498 + c3 * rho3 00499 + cs * rhos 00500 + cl * rhoc; 00501 00502 return potential; 00503 }
Definition at line 68 of file G4QMDMeanField.hh.
Referenced by G4QMDCollision::CalKinematicsOfBinaryCollisions().
G4QMDSystem* G4QMDMeanField::GetSystem | ( | ) | [inline] |
G4double G4QMDMeanField::GetTotalPotential | ( | ) |
Definition at line 507 of file G4QMDMeanField.cc.
References G4QMDParticipant::GetChargeInUnitOfEplus(), G4QMDParticipant::GetNuc(), G4QMDSystem::GetParticipant(), G4QMDSystem::GetTotalNumberOfParticipant(), and CLHEP::detail::n.
Referenced by G4QMDReaction::ApplyYourself(), G4QMDCollision::CalFinalStateOfTheBinaryCollision(), G4QMDCollision::CalFinalStateOfTheBinaryCollisionJQMD(), G4QMDCollision::CalKinematicsOfBinaryCollisions(), and SetNucleus().
00508 { 00509 00510 G4int n = system->GetTotalNumberOfParticipant(); 00511 00512 std::vector < G4double > rhoa ( n , 0.0 ); 00513 std::vector < G4double > rho3 ( n , 0.0 ); 00514 std::vector < G4double > rhos ( n , 0.0 ); 00515 std::vector < G4double > rhoc ( n , 0.0 ); 00516 00517 00518 for ( G4int i = 0 ; i < n ; i ++ ) 00519 { 00520 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 00521 G4int inuc = system->GetParticipant(i)->GetNuc(); 00522 00523 for ( G4int j = 0 ; j < n ; j ++ ) 00524 { 00525 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 00526 G4int jnuc = system->GetParticipant(j)->GetNuc(); 00527 00528 rhoa[i] += rha[j][i]; 00529 rhoc[i] += rhe[j][i]; 00530 rhos[i] += rha[j][i] * jnuc * inuc 00531 * ( 1 - 2 * std::abs ( jcharge - icharge ) ); 00532 } 00533 00534 rho3[i] = std::pow ( rhoa[i] , gamm ); 00535 } 00536 00537 G4double potential = c0 * std::accumulate( rhoa.begin() , rhoa.end() , 0.0 ) 00538 + c3 * std::accumulate( rho3.begin() , rho3.end() , 0.0 ) 00539 + cs * std::accumulate( rhos.begin() , rhos.end() , 0.0 ) 00540 + cl * std::accumulate( rhoc.begin() , rhoc.end() , 0.0 ); 00541 00542 return potential; 00543 00544 }
Definition at line 602 of file G4QMDMeanField.cc.
References G4UniformRand, G4QMDParticipant::GetNuc(), and G4QMDSystem::GetParticipant().
Referenced by G4QMDCollision::CalFinalStateOfTheBinaryCollision(), and G4QMDCollision::CalKinematicsOfBinaryCollisions().
00603 { 00604 G4bool result = false; 00605 00606 if ( system->GetParticipant( i )->GetNuc() == 1 ) 00607 { 00608 G4double pf = calPauliBlockingFactor( i ); 00609 G4double rand = G4UniformRand(); 00610 if ( pf > rand ) result = true; 00611 } 00612 00613 return result; 00614 }
void G4QMDMeanField::SetNucleus | ( | G4QMDNucleus * | aSystem | ) |
Definition at line 132 of file G4QMDMeanField.cc.
References G4QMDNucleus::CalEnergyAndAngularMomentumInCM(), GetTotalPotential(), SetSystem(), and G4QMDNucleus::SetTotalPotential().
Referenced by G4QMDReaction::ApplyYourself().
00133 { 00134 00135 //std::cout << "QMDMeanField SetNucleus" << std::endl; 00136 00137 SetSystem( aNucleus ); 00138 00139 G4double totalPotential = GetTotalPotential(); 00140 aNucleus->SetTotalPotential( totalPotential ); 00141 00142 aNucleus->CalEnergyAndAngularMomentumInCM(); 00143 00144 }
void G4QMDMeanField::SetSystem | ( | G4QMDSystem * | aSystem | ) |
Definition at line 86 of file G4QMDMeanField.cc.
References Cal2BodyQuantities(), G4QMDSystem::GetTotalNumberOfParticipant(), and CLHEP::detail::n.
Referenced by G4QMDReaction::ApplyYourself(), G4QMDGroundStateNucleus::G4QMDGroundStateNucleus(), SetNucleus(), and Update().
00087 { 00088 00089 //std::cout << "QMDMeanField SetSystem" << std::endl; 00090 00091 system = aSystem; 00092 00093 G4int n = system->GetTotalNumberOfParticipant(); 00094 00095 pp2.clear(); 00096 rr2.clear(); 00097 rbij.clear(); 00098 rha.clear(); 00099 rhe.clear(); 00100 rhc.clear(); 00101 00102 rr2.resize( n ); 00103 pp2.resize( n ); 00104 rbij.resize( n ); 00105 rha.resize( n ); 00106 rhe.resize( n ); 00107 rhc.resize( n ); 00108 00109 for ( int i = 0 ; i < n ; i++ ) 00110 { 00111 rr2[i].resize( n ); 00112 pp2[i].resize( n ); 00113 rbij[i].resize( n ); 00114 rha[i].resize( n ); 00115 rhe[i].resize( n ); 00116 rhc[i].resize( n ); 00117 } 00118 00119 00120 ffr.clear(); 00121 ffp.clear(); 00122 rh3d.clear(); 00123 00124 ffr.resize( n ); 00125 ffp.resize( n ); 00126 rh3d.resize( n ); 00127 00128 Cal2BodyQuantities(); 00129 00130 }
void G4QMDMeanField::Update | ( | ) |
Definition at line 950 of file G4QMDMeanField.cc.
References SetSystem().
Referenced by G4QMDCollision::CalFinalStateOfTheBinaryCollision(), and G4QMDCollision::CalKinematicsOfBinaryCollisions().
00951 { 00952 SetSystem( system ); 00953 }