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