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