Geant4 Cross Reference |
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 // 081120 Add Update by T. Koi 27 // 28 // 230307 Skyrme-QMD parameters added by Y-H. Sato and A. Haga 29 // 230307 "CalDensityProfile" and "CalChargeDensityProfile" functions added by Y-H. Sato and A. Haga 30 // 230307 "GetSingleEnergy" and "GetTotalEnergy" functions added by Y-H. Sato and A. Haga 31 32 #include <map> 33 #include <algorithm> 34 #include <numeric> 35 36 #include <cmath> 37 #include <CLHEP/Random/Stat.h> 38 39 #include "G4LightIonQMDMeanField.hh" 40 #include "G4LightIonQMDParameters.hh" 41 #include "G4Exp.hh" 42 #include "G4Pow.hh" 43 #include "G4PhysicalConstants.hh" 44 #include "Randomize.hh" 45 46 G4LightIonQMDMeanField::G4LightIonQMDMeanField() 47 { 48 G4LightIonQMDParameters* parameters = G4LightIonQMDParameters::GetInstance(); 49 wl = parameters->Get_wl(); 50 cl = parameters->Get_cl(); 51 rho0 = parameters->Get_rho0(); 52 hbc = parameters->Get_hbc(); 53 gamm = parameters->Get_gamm(); 54 eta = parameters->Get_eta(); // Skyrme-QMD 55 kappas = parameters->Get_kappas(); // Skyrme-QMD 56 57 cpw = parameters->Get_cpw(); 58 cph = parameters->Get_cph(); 59 cpc = parameters->Get_cpc(); 60 61 c0 = parameters->Get_c0(); 62 c3 = parameters->Get_c3(); 63 cs = parameters->Get_cs(); 64 g0 = parameters->Get_g0(); // Skyrme-QMD 65 g0iso = parameters->Get_g0iso(); // Skyrme-QMD 66 gtau0 = parameters->Get_gtau0(); // Skyrme-QMD 67 68 // distance 69 c0w = 1.0/4.0/wl; 70 c0sw = std::sqrt( c0w ); 71 clw = 2.0 / std::sqrt ( 4.0 * pi * wl ); 72 73 // graduate 74 c0g = - c0 / ( 2.0 * wl ); 75 c3g = - c3 / ( 4.0 * wl ) * gamm; 76 csg = - cs / ( 2.0 * wl ); 77 pag = gamm - 1; 78 pag_tau = eta - 1; // Skyrme-QMD 79 cg0 = - g0 / ( 2.0 * wl ); // Skyrme-QMD 80 cgtau0 = - gtau0 / ( 4.0 * wl ) * eta; // Skyrme-QMD 81 82 system = nullptr; // will be set through SetSystem method 83 } 84 85 void G4LightIonQMDMeanField::SetSystem ( G4QMDSystem* aSystem ) 86 { 87 system = aSystem; 88 89 G4int n = system->GetTotalNumberOfParticipant(); 90 91 pp2.clear(); 92 rr2.clear(); 93 rbij.clear(); 94 rha.clear(); 95 rhe.clear(); 96 rhc.clear(); 97 98 rr2.resize( n ); 99 pp2.resize( n ); 100 rbij.resize( n ); 101 rha.resize( n ); 102 rhe.resize( n ); 103 rhc.resize( n ); 104 105 for ( G4int i = 0 ; i < n ; ++i ) 106 { 107 rr2[i].resize( n ); 108 pp2[i].resize( n ); 109 rbij[i].resize( n ); 110 rha[i].resize( n ); 111 rhe[i].resize( n ); 112 rhc[i].resize( n ); 113 } 114 115 ffr.clear(); 116 ffp.clear(); 117 rh3d.clear(); 118 rh3d_tau.clear(); // Skyrme-QMD 119 120 ffr.resize( n ); 121 ffp.resize( n ); 122 rh3d.resize( n ); 123 rh3d_tau.resize( n ); // Skyrme-QMD 124 125 Cal2BodyQuantities(); 126 } 127 128 void G4LightIonQMDMeanField::SetNucleus ( G4LightIonQMDNucleus* aNucleus ) 129 { 130 SetSystem( aNucleus ); 131 132 G4double totalPotential = GetTotalPotential(); 133 aNucleus->SetTotalPotential( totalPotential ); 134 aNucleus->CalEnergyAndAngularMomentumInCM(); 135 } 136 137 void G4LightIonQMDMeanField::Cal2BodyQuantities() 138 { 139 if ( system->GetTotalNumberOfParticipant() < 2 ) { return; } 140 141 for ( G4int j = 1 ; j < system->GetTotalNumberOfParticipant() ; ++j ) 142 { 143 G4ThreeVector rj = system->GetParticipant( j )->GetPosition(); 144 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum(); 145 146 for ( G4int i = 0 ; i < j ; ++i ) 147 { 148 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 149 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum(); 150 151 G4ThreeVector rij = ri - rj; 152 G4ThreeVector pij = (p4i - p4j).v(); 153 G4LorentzVector p4ij = p4i - p4j; 154 G4ThreeVector bij = ( p4i + p4j ).boostVector(); 155 G4double gammaij = ( p4i + p4j ).gamma(); 156 157 G4double eij = ( p4i + p4j ).e(); 158 159 G4double rbrb = rij*bij; 160 G4double rij2 = rij*rij; 161 G4double pij2 = pij*pij; 162 163 rbrb = irelcr * rbrb; 164 G4double gamma2_ij = gammaij*gammaij; 165 166 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb; 167 rr2[j][i] = rr2[i][j]; 168 169 rbij[i][j] = gamma2_ij * rbrb; 170 rbij[j][i] = - rbij[i][j]; 171 172 pp2[i][j] = pij2 173 + irelcr * ( - G4Pow::GetInstance()->powN ( p4i.e() - p4j.e() , 2 ) 174 + gamma2_ij * G4Pow::GetInstance()->powN ( ( ( p4i.m2() - p4j.m2() ) / eij ) , 2 ) ); 175 176 177 pp2[j][i] = pp2[i][j]; 178 179 // Gauss term 180 181 G4double expa1 = - rr2[i][j] * c0w; 182 183 G4double rh1; 184 if ( expa1 > epsx ) 185 { 186 rh1 = G4Exp( expa1 ); 187 } 188 else 189 { 190 rh1 = 0.0; 191 } 192 193 G4int ibry = system->GetParticipant(i)->GetBaryonNumber(); 194 G4int jbry = system->GetParticipant(j)->GetBaryonNumber(); 195 196 rha[i][j] = ibry*jbry*rh1; 197 rha[j][i] = rha[i][j]; 198 199 // Coulomb terms 200 201 G4double rrs2 = rr2[i][j] + epscl; 202 G4double rrs = std::sqrt ( rrs2 ); 203 204 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 205 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 206 207 G4double xerf = 0.0; 208 // T. K. add this protection. 5.8 is good enough for double 209 if ( rrs*c0sw < 5.8 ) 210 { 211 #if defined WIN32-VC 212 xerf = CLHEP::HepStat::erf ( rrs*c0sw ); 213 #else 214 xerf = std::erf ( rrs*c0sw ); 215 #endif 216 } 217 else 218 { 219 xerf = 1.0; 220 } 221 222 G4double erfij = xerf/rrs; 223 224 rhe[i][j] = icharge*jcharge * erfij; 225 rhe[j][i] = rhe[i][j]; 226 rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2; 227 rhc[j][i] = rhc[i][j]; 228 } // i 229 } // j 230 } 231 232 void G4LightIonQMDMeanField::Cal2BodyQuantities( G4int i ) 233 { 234 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 235 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum(); 236 237 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j ) 238 { 239 if ( j == i ) { continue; } 240 241 G4ThreeVector rj = system->GetParticipant( j )->GetPosition(); 242 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum(); 243 244 G4ThreeVector rij = ri - rj; 245 G4ThreeVector pij = (p4i - p4j).v(); 246 G4LorentzVector p4ij = p4i - p4j; 247 G4ThreeVector bij = ( p4i + p4j ).boostVector(); 248 G4double gammaij = ( p4i + p4j ).gamma(); 249 250 G4double eij = ( p4i + p4j ).e(); 251 252 G4double rbrb = rij*bij; 253 G4double rij2 = rij*rij; 254 G4double pij2 = pij*pij; 255 256 rbrb = irelcr * rbrb; 257 G4double gamma2_ij = gammaij*gammaij; 258 259 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb; 260 rr2[j][i] = rr2[i][j]; 261 262 rbij[i][j] = gamma2_ij * rbrb; 263 rbij[j][i] = - rbij[i][j]; 264 265 pp2[i][j] = pij2 266 + irelcr * ( - G4Pow::GetInstance()->powN ( p4i.e() - p4j.e() , 2 ) 267 + gamma2_ij * G4Pow::GetInstance()->powN ( ( ( p4i.m2() - p4j.m2() ) / eij ) , 2 ) ); 268 269 pp2[j][i] = pp2[i][j]; 270 271 // Gauss term 272 273 G4double expa1 = - rr2[i][j] * c0w; 274 275 G4double rh1; 276 if ( expa1 > epsx ) 277 { 278 rh1 = G4Exp( expa1 ); 279 } 280 else 281 { 282 rh1 = 0.0; 283 } 284 285 G4int ibry = system->GetParticipant(i)->GetBaryonNumber(); 286 G4int jbry = system->GetParticipant(j)->GetBaryonNumber(); 287 288 rha[i][j] = ibry*jbry*rh1; 289 rha[j][i] = rha[i][j]; 290 291 // Coulomb terms 292 293 G4double rrs2 = rr2[i][j] + epscl; 294 G4double rrs = std::sqrt ( rrs2 ); 295 296 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 297 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 298 299 G4double xerf = 0.0; 300 // T. K. add this protection. 5.8 is good enough for double 301 if ( rrs*c0sw < 5.8 ) 302 { 303 #if defined WIN32-VC 304 xerf = CLHEP::HepStat::erf ( rrs*c0sw ); 305 #else 306 xerf = std::erf ( rrs*c0sw ); 307 #endif 308 } 309 else 310 { 311 xerf = 1.0; 312 } 313 314 G4double erfij = xerf/rrs; 315 316 rhe[i][j] = icharge*jcharge * erfij; 317 rhe[j][i] = rhe[i][j]; 318 rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2; 319 rhc[j][i] = rhc[i][j]; 320 } 321 } 322 323 void G4LightIonQMDMeanField::CalGraduate() 324 { 325 ffr.resize( system->GetTotalNumberOfParticipant() ); 326 ffp.resize( system->GetTotalNumberOfParticipant() ); 327 rh3d.resize( system->GetTotalNumberOfParticipant() ); 328 rh3d_tau.resize( system->GetTotalNumberOfParticipant() ); // Skyrme-QMD 329 330 for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; ++i ) 331 { 332 G4double rho3 = 0.0; 333 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j ) 334 { 335 rho3 += rha[j][i]; 336 } 337 rh3d[i] = G4Pow::GetInstance()->powA ( rho3 , pag ); 338 rh3d_tau[i] = G4Pow::GetInstance()->powA ( rho3 , pag_tau ); // Skyrme-QMD 339 } 340 341 for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; ++i ) 342 { 343 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 344 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum(); 345 346 G4ThreeVector betai = p4i.v()/p4i.e(); 347 348 // R-JQMD 349 G4double Vi = GetPotential( i ); 350 G4double p_zero = std::sqrt( p4i.e()*p4i.e() + 2*p4i.m()*Vi); 351 G4ThreeVector betai_R = p4i.v()/p_zero; 352 G4double mi_R = p4i.m()/p_zero; 353 354 ffr[i] = betai_R; 355 ffp[i] = G4ThreeVector( 0.0 ); 356 357 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j ) 358 { 359 G4ThreeVector rj = system->GetParticipant( j )->GetPosition(); 360 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum(); 361 362 G4double eij = p4i.e() + p4j.e(); 363 364 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 365 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 366 367 G4int inuc = system->GetParticipant(i)->GetNuc(); 368 G4int jnuc = system->GetParticipant(j)->GetNuc(); 369 370 G4double fsij = 3.0/(2*wl) - rr2[j][i]/(2*wl)/(2*wl); // Add for Skyrme-QMD 371 372 G4double ccpp = c0g * rha[j][i] 373 + c3g * rha[j][i] * ( rh3d[j] + rh3d[i] ) 374 + cg0 * rha[j][i]/wl 375 + cg0 * rha[j][i] * fsij 376 + cgtau0 * rha[j][i] * ( rh3d_tau[j] + rh3d_tau[i] ) 377 + csg * rha[j][i] * jnuc * inuc 378 * ( 1. - 2. * std::abs( jcharge - icharge ) ) 379 * (1. - kappas * fsij + kappas / wl) 380 + cl * rhc[j][i]; 381 382 ccpp *= mi_R; 383 384 G4double grbb = - rbij[j][i]; 385 G4double ccrr = grbb * ccpp / eij; 386 387 G4ThreeVector rij = ri - rj; 388 G4ThreeVector betaij = ( p4i + p4j ).v()/eij; 389 G4ThreeVector cij = betaij - betai; 390 391 ffr[i] = ffr[i] + 2*ccrr* ( rij + grbb*cij ); 392 ffp[i] = ffp[i] - 2*ccpp* ( rij + grbb*betaij ); 393 } 394 } 395 } 396 397 G4double G4LightIonQMDMeanField::GetPotential( G4int i ) 398 { 399 G4int n = system->GetTotalNumberOfParticipant(); 400 401 G4double rhoa = 0.0; 402 G4double rho3 = 0.0; 403 G4double fsij_rhoa = 0.0; // Skyrme-QMD 404 //G4double fsij_rhos = 0.0; // Skyrme-QMD 405 G4double rho3_tau = 0.0; // Skyrme-QMD 406 G4double rhos = 0.0; 407 G4double rhoc = 0.0; 408 409 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 410 G4int inuc = system->GetParticipant(i)->GetNuc(); 411 412 for ( G4int j = 0 ; j < n ; ++j ) 413 { 414 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 415 G4int jnuc = system->GetParticipant(j)->GetNuc(); 416 G4double fsij = 3.0/(2*wl) - rr2[j][i]/(2*wl)/(2*wl); // Add for Skyrme-QMD 417 418 rhoa += rha[j][i]; 419 fsij_rhoa += fsij * rha[j][i]; // Skyrme-QMD 420 rhoc += rhe[j][i]; 421 rhos += rha[j][i] * jnuc * inuc 422 * ( 1. - 2. * std::abs( jcharge - icharge ) ) // Skyrme-QMD 423 * (1. - kappas * fsij); // Skyrme-QMD 424 } 425 426 rho3 = G4Pow::GetInstance()->powA ( rhoa , gamm ); 427 rho3_tau = G4Pow::GetInstance()->powA ( rhoa , eta ); 428 429 G4double potential = c0 * rhoa 430 + c3 * rho3 431 + g0 * fsij_rhoa // Skyrme-QMD 432 //+ g0iso * fsij_rhos // Skyrme-QMD 433 + gtau0 * rho3_tau // Skyrme-QMD 434 + cs * rhos 435 + cl * rhoc; 436 return potential; 437 } 438 439 G4double G4LightIonQMDMeanField::GetTotalPotential() 440 { 441 G4int n = system->GetTotalNumberOfParticipant(); 442 443 std::vector < G4double > rhoa ( n , 0.0 ); 444 std::vector < G4double > rho3 ( n , 0.0 ); 445 std::vector < G4double > rho3_tau ( n , 0.0 ); // Skyrme-QMD 446 //std::vector < G4double > fsij_rhos ( n , 0.0 ); // Skyrme-QMD 447 std::vector < G4double > fsij_rhoa ( n , 0.0 ); // Skyrme-QMD 448 std::vector < G4double > rhos ( n , 0.0 ); 449 std::vector < G4double > rhoc ( n , 0.0 ); 450 451 for ( G4int i = 0 ; i < n ; ++i ) 452 { 453 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 454 G4int inuc = system->GetParticipant(i)->GetNuc(); 455 456 for ( G4int j = 0 ; j < n ; ++j ) 457 { 458 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 459 G4int jnuc = system->GetParticipant(j)->GetNuc(); 460 G4double fsij = 3.0/(2*wl) - rr2[j][i]/(2*wl)/(2*wl); // Add for Skyrme-QMD 461 462 rhoa[i] += rha[j][i]; 463 fsij_rhoa[i] += fsij * rha[j][i]; // Skyrme-QMD 464 rhoc[i] += rhe[j][i]; 465 rhos[i] += rha[j][i] * jnuc * inuc 466 //* ( 1 - 2 * std::abs ( jcharge - icharge ) ); 467 * ( 1. - 2. * std::abs( jcharge - icharge ) ) // Skyrme-QMD 468 * (1. - kappas * fsij); // Skyrme-QMD 469 //fsij_rhos[i] += fsij * rha[j][i] * jnuc * inuc 470 //* ( 1. - 2. * std::abs( jcharge - icharge ) ) // Skyrme-QMD 471 //* (1. - kappas * fsij); // Skyrme-QMD 472 } 473 474 rho3[i] = G4Pow::GetInstance()->powA ( rhoa[i] , gamm ); 475 rho3_tau[i] = G4Pow::GetInstance()->powA ( rhoa[i] , eta ); 476 } 477 478 G4double potential = c0 * std::accumulate( rhoa.cbegin() , rhoa.cend() , 0.0 ) 479 + c3 * std::accumulate( rho3.cbegin() , rho3.cend() , 0.0 ) 480 + g0 * std::accumulate( fsij_rhoa.cbegin() , fsij_rhoa.cend() , 0.0 ) 481 //+ g0iso * std::accumulate( fsij_rhos.cbegin() , fsij_rhos.cend() , 0.0 ) 482 + gtau0 * std::accumulate( rho3_tau.cbegin() , rho3_tau.cend() , 0.0 ) 483 + cs * std::accumulate( rhos.cbegin() , rhos.cend() , 0.0 ) 484 + cl * std::accumulate( rhoc.cbegin() , rhoc.cend() , 0.0 ); 485 486 return potential; 487 } 488 489 G4double G4LightIonQMDMeanField::GetSingleEnergy( G4int j ) 490 { 491 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum(); 492 G4double emass = p4j.m(); 493 G4double ekinal2 = p4j.e()*p4j.e(); 494 G4double esingle = std::sqrt(ekinal2 + 2*emass*GetPotential(j)); 495 return esingle; 496 } 497 498 G4double G4LightIonQMDMeanField::GetTotalEnergy() 499 { 500 501 G4int n = system->GetTotalNumberOfParticipant(); 502 G4double etotal = 0.0; 503 for ( int j = 0 ; j < n ; j++ ) 504 { 505 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum(); 506 G4double emass = p4j.m(); 507 G4double ekinal2 = p4j.e()*p4j.e(); 508 etotal += std::sqrt(ekinal2 + 2*emass*GetPotential(j)); 509 } 510 return etotal; 511 512 } 513 514 G4double G4LightIonQMDMeanField::calPauliBlockingFactor( G4int i ) 515 { 516 // i is supposed beyond total number of Participant() 517 518 G4double pf = 0.0; 519 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 520 521 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j ) 522 { 523 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 524 G4int jnuc = system->GetParticipant(j)->GetNuc(); 525 526 if ( jcharge == icharge && jnuc == 1 ) 527 { 528 G4double expa = -rr2[i][j]*cpw; 529 if ( expa > epsx ) 530 { 531 expa = expa - pp2[i][j]*cph; 532 if ( expa > epsx ) 533 { 534 pf = pf + G4Exp ( expa ); 535 } 536 } 537 } 538 } 539 540 return ( pf - 1.0 ) * cpc; 541 } 542 543 G4bool G4LightIonQMDMeanField::IsPauliBlocked( G4int i ) 544 { 545 G4bool result = false; 546 547 if ( system->GetParticipant( i )->GetNuc() == 1 ) 548 { 549 G4double pf = calPauliBlockingFactor( i ); 550 G4double rand = G4UniformRand(); 551 if ( pf > rand ) { result = true; } 552 } 553 554 return result; 555 } 556 557 void G4LightIonQMDMeanField::DoPropagation( G4double dt ) 558 { 559 G4double cc2 = 1.0; 560 G4double cc1 = 1.0 - cc2; 561 G4double cc3 = 1.0 / 2.0 / cc2; 562 563 G4double dt3 = dt * cc3; 564 G4double dt1 = dt * ( cc1 - cc3 ); 565 G4double dt2 = dt * cc2; 566 567 CalGraduate(); 568 569 G4int n = system->GetTotalNumberOfParticipant(); 570 571 // 1st Step 572 573 std::vector< G4ThreeVector > f0r, f0p; 574 f0r.resize( n ); 575 f0p.resize( n ); 576 577 for ( G4int i = 0 ; i < n ; ++i ) 578 { 579 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 580 G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum(); 581 582 ri += dt3* ffr[i]; 583 p3i += dt3* ffp[i]; 584 585 f0r[i] = ffr[i]; 586 f0p[i] = ffp[i]; 587 588 system->GetParticipant( i )->SetPosition( ri ); 589 system->GetParticipant( i )->SetMomentum( p3i ); 590 591 // we do not need set total momentum by ourselvs 592 } 593 594 // 2nd Step 595 596 Cal2BodyQuantities(); 597 CalGraduate(); 598 599 for ( G4int i = 0 ; i < n ; ++i ) 600 { 601 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 602 G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum(); 603 604 ri += dt1* f0r[i] + dt2* ffr[i]; 605 p3i += dt1* f0p[i] + dt2* ffp[i]; 606 607 system->GetParticipant( i )->SetPosition( ri ); 608 system->GetParticipant( i )->SetMomentum( p3i ); 609 610 // we do not need set total momentum by ourselvs 611 } 612 613 Cal2BodyQuantities(); 614 } 615 616 std::vector< G4LightIonQMDNucleus* > G4LightIonQMDMeanField::DoClusterJudgment() 617 { 618 Cal2BodyQuantities(); 619 620 G4double cpf2 = G4Pow::GetInstance()->A23 ( 1.5 * pi*pi * G4Pow::GetInstance()->powA ( 4.0 * pi * wl , -1.5 ) ) * hbc * hbc; 621 G4double rcc2 = rclds*rclds; 622 623 G4int n = system->GetTotalNumberOfParticipant(); 624 std::vector < G4double > rhoa; 625 rhoa.resize ( n ); 626 627 for ( G4int i = 0 ; i < n ; ++i ) 628 { 629 rhoa[i] = 0.0; 630 631 if ( system->GetParticipant( i )->GetBaryonNumber() == 1 ) 632 { 633 for ( G4int j = 0 ; j < n ; ++j ) 634 { 635 if ( system->GetParticipant( j )->GetBaryonNumber() == 1 ) 636 rhoa[i] += rha[i][j]; 637 } 638 } 639 640 rhoa[i] = G4Pow::GetInstance()->A13 ( rhoa[i] + 1 ); 641 } 642 643 // identification of the cluster 644 std::vector < G4bool > is_already_belong_some_cluster; 645 646 // cluster_id participant_id 647 std::multimap < G4int , G4int > comb_map; 648 std::multimap < G4int , G4int > assign_map; 649 assign_map.clear(); 650 651 std::vector < G4int > mascl; 652 std::vector < G4int > num; 653 mascl.resize ( n ); 654 num.resize ( n ); 655 is_already_belong_some_cluster.resize ( n ); 656 657 std::vector < G4int > is_assigned_to ( n , -1 ); 658 std::multimap < G4int , G4int > clusters; 659 660 for ( G4int i = 0 ; i < n ; ++i ) 661 { 662 mascl[i] = 1; 663 num[i] = 1; 664 is_already_belong_some_cluster[i] = false; 665 } 666 667 G4int ichek = 1; 668 G4int id = 0; 669 G4int cluster_id = -1; 670 for ( G4int i = 0 ; i < n-1 ; ++i ) 671 { 672 G4bool hasThisCompany = false; 673 674 if ( system->GetParticipant( i )->GetBaryonNumber() == 1 ) 675 { 676 G4int j1 = i + 1; 677 for ( G4int j = j1 ; j < n ; ++j ) 678 { 679 std::vector < G4int > cluster_participants; 680 if ( system->GetParticipant( j )->GetBaryonNumber() == 1 ) 681 { 682 G4double rdist2 = rr2[ i ][ j ]; 683 G4double pdist2 = pp2[ i ][ j ]; 684 G4double pcc2 = cpf2 685 * ( rhoa[ i ] + rhoa[ j ] ) 686 * ( rhoa[ i ] + rhoa[ j ] ); 687 688 // Check phase space: close enough? 689 if ( rdist2 < rcc2 && pdist2 < pcc2 ) 690 { 691 if ( is_assigned_to [ j ] == -1 ) 692 { 693 if ( is_assigned_to [ i ] == -1 ) 694 { 695 if ( clusters.size() != 0 ) 696 { 697 id = clusters.rbegin()->first + 1; 698 } 699 else 700 { 701 id = 0; 702 } 703 clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , i ) ); 704 is_assigned_to [ i ] = id; 705 clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , j ) ); 706 is_assigned_to [ j ] = id; 707 } 708 else 709 { 710 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ i ] , j ) ); 711 is_assigned_to [ j ] = is_assigned_to [ i ]; 712 } 713 } 714 else 715 { 716 // j is already belong to some cluster 717 if ( is_assigned_to [ i ] == -1 ) 718 { 719 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , i ) ); 720 is_assigned_to [ i ] = is_assigned_to [ j ]; 721 } 722 else 723 { 724 // i has companion 725 if ( is_assigned_to [ i ] != is_assigned_to [ j ] ) 726 { 727 // move companions to the cluster 728 std::multimap< G4int , G4int > clusters_tmp; 729 G4int target_cluster_id; 730 if ( is_assigned_to [ i ] > is_assigned_to [ j ] ) 731 { 732 target_cluster_id = is_assigned_to [ i ]; 733 } 734 else 735 { 736 target_cluster_id = is_assigned_to [ j ]; 737 } 738 for ( auto it = clusters.cbegin() ; it != clusters.cend() ; ++it ) 739 { 740 if ( it->first == target_cluster_id ) 741 { 742 is_assigned_to [ it->second ] = is_assigned_to [ j ]; 743 clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , it->second ) ); 744 } 745 else 746 { 747 clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( it->first , it->second ) ); 748 } 749 } 750 clusters = std::move(clusters_tmp); 751 } 752 } 753 } 754 755 comb_map.insert( std::multimap<G4int,G4int>::value_type ( i , j ) ); 756 cluster_participants.push_back ( j ); 757 758 if ( assign_map.find( cluster_id ) == assign_map.end() ) 759 { 760 is_already_belong_some_cluster[i] = true; 761 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , i ) ); 762 hasThisCompany = true; 763 } 764 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , j ) ); 765 is_already_belong_some_cluster[j] = true; 766 } 767 768 if ( ichek == i ) 769 { 770 ++ichek; 771 } 772 } 773 } 774 } 775 if ( hasThisCompany == true ) { ++cluster_id; } 776 } 777 778 // sort 779 // Heavy cluster comes first 780 // size cluster_id 781 std::multimap< G4int , G4int > sorted_cluster_map; 782 for ( G4int i = 0 ; i <= id ; ++i ) // << "<=" because id is highest cluster nubmer. 783 { 784 sorted_cluster_map.insert ( std::multimap<G4int,G4int>::value_type ( (G4int) clusters.count( i ) , i ) ); 785 } 786 787 // create nucleus from divided clusters 788 std::vector < G4LightIonQMDNucleus* > result; 789 for ( auto it = sorted_cluster_map.crbegin(); it != sorted_cluster_map.crend(); ++it ) 790 { 791 if ( it->first != 0 ) 792 { 793 G4LightIonQMDNucleus* nucleus = new G4LightIonQMDNucleus(); 794 for ( auto itt = clusters.cbegin(); itt != clusters.cend(); ++itt ) 795 { 796 if ( it->second == itt->first ) 797 { 798 nucleus->SetParticipant( system->GetParticipant ( itt->second ) ); 799 } 800 } 801 result.push_back( nucleus ); 802 } 803 } 804 805 // delete participants from current system 806 for ( auto it = result.cbegin(); it != result.cend(); ++it ) 807 { 808 system->SubtractSystem ( *it ); 809 } 810 811 return result; 812 } 813 814 void G4LightIonQMDMeanField::Update() 815 { 816 SetSystem( system ); 817 } 818