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 #include "G4Proton.hh" 31 #include "G4Proton.hh" 32 #include "G4Neutron.hh" 32 #include "G4Neutron.hh" 33 #include "G4Exp.hh" 33 #include "G4Exp.hh" 34 #include "G4Pow.hh" 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 72 // Following 10 lines should be here, right before the line 90. 74 // Otherwise, mass number cannot be conserv 73 // Otherwise, mass number cannot be conserved if the projectile or 75 // the target are nucleons. 74 // the target are nucleons. 76 //Nucleon primary or target case; 75 //Nucleon primary or target case; 77 if ( z == 1 && a == 1 ) { // Hydrogen Cas 76 if ( z == 1 && a == 1 ) { // Hydrogen Case or proton primary 78 SetParticipant( new G4QMDParticipant( G4 77 SetParticipant( new G4QMDParticipant( G4Proton::Proton() , G4ThreeVector( 0.0 ) , G4ThreeVector( 0.0 ) ) ); 79 // ebini = 0.0; << 78 ebini = 0.0; 80 return; 79 return; 81 } 80 } 82 else if ( z == 0 && a == 1 ) { // Neutron p 81 else if ( z == 0 && a == 1 ) { // Neutron primary 83 SetParticipant( new G4QMDParticipant( G4 82 SetParticipant( new G4QMDParticipant( G4Neutron::Neutron() , G4ThreeVector( 0.0 ) , G4ThreeVector( 0.0 ) ) ); 84 // ebini = 0.0; << 83 ebini = 0.0; 85 return; 84 return; 86 } 85 } 87 86 88 87 89 //edepth = 0.0; 88 //edepth = 0.0; 90 89 91 for ( G4int i = 0 ; i < a ; ++i ) << 90 for ( int i = 0 ; i < a ; i++ ) 92 { 91 { 93 92 94 G4ParticleDefinition* pd; 93 G4ParticleDefinition* pd; 95 94 96 if ( i < z ) 95 if ( i < z ) 97 { 96 { 98 pd = G4Proton::Proton(); 97 pd = G4Proton::Proton(); 99 } 98 } 100 else 99 else 101 { 100 { 102 pd = G4Neutron::Neutron(); 101 pd = G4Neutron::Neutron(); 103 } 102 } 104 103 105 G4ThreeVector p( 0.0 ); 104 G4ThreeVector p( 0.0 ); 106 G4ThreeVector r( 0.0 ); 105 G4ThreeVector r( 0.0 ); 107 G4QMDParticipant* aParticipant = new G4Q 106 G4QMDParticipant* aParticipant = new G4QMDParticipant( pd , p , r ); 108 SetParticipant( aParticipant ); 107 SetParticipant( aParticipant ); 109 108 110 } 109 } 111 110 112 G4double radious = r00 * G4Pow::GetInstance << 111 G4double radious = r00 * G4Pow::GetInstance()->A13( double ( GetMassNumber() ) ); 113 112 114 rt00 = radious - r01; 113 rt00 = radious - r01; 115 radm = radious - rada * ( gamm - 1.0 ) + ra 114 radm = radious - rada * ( gamm - 1.0 ) + radb; 116 rmax = 1.0 / ( 1.0 + G4Exp ( -rt00/saa ) ); 115 rmax = 1.0 / ( 1.0 + G4Exp ( -rt00/saa ) ); 117 116 118 //maxTrial = 1000; 117 //maxTrial = 1000; 119 118 120 119 121 meanfield = new G4QMDMeanField(); 120 meanfield = new G4QMDMeanField(); 122 meanfield->SetSystem( this ); 121 meanfield->SetSystem( this ); 123 122 124 //std::cout << "G4QMDGroundStateNucleus( G4 123 //std::cout << "G4QMDGroundStateNucleus( G4int z , G4int a ) packNucleons Begin ( z , a ) ( " << z << ", " << a << " )" << std::endl; 125 packNucleons(); 124 packNucleons(); 126 //std::cout << "G4QMDGroundStateNucleus( G4 125 //std::cout << "G4QMDGroundStateNucleus( G4int z , G4int a ) packNucleons End" << std::endl; 127 126 128 delete meanfield; 127 delete meanfield; 129 128 130 } 129 } 131 130 132 131 133 132 134 void G4QMDGroundStateNucleus::packNucleons() 133 void G4QMDGroundStateNucleus::packNucleons() 135 { 134 { 136 135 137 //std::cout << "G4QMDGroundStateNucleus::pa 136 //std::cout << "G4QMDGroundStateNucleus::packNucleons" << std::endl; 138 137 139 ebini = - G4NucleiProperties::GetBindingEne 138 ebini = - G4NucleiProperties::GetBindingEnergy( GetMassNumber() , GetAtomicNumber() ) / GetMassNumber(); 140 139 141 G4double ebin00 = ebini * 0.001; 140 G4double ebin00 = ebini * 0.001; 142 141 143 G4double ebin0 = 0.0; 142 G4double ebin0 = 0.0; 144 G4double ebin1 = 0.0; 143 G4double ebin1 = 0.0; 145 144 146 if ( GetMassNumber() != 4 ) 145 if ( GetMassNumber() != 4 ) 147 { 146 { 148 ebin0 = ( ebini - 0.5 ) * 0.001; 147 ebin0 = ( ebini - 0.5 ) * 0.001; 149 ebin1 = ( ebini + 0.5 ) * 0.001; 148 ebin1 = ( ebini + 0.5 ) * 0.001; 150 } 149 } 151 else 150 else 152 { 151 { 153 ebin0 = ( ebini - 1.5 ) * 0.001; 152 ebin0 = ( ebini - 1.5 ) * 0.001; 154 ebin1 = ( ebini + 1.5 ) * 0.001; 153 ebin1 = ( ebini + 1.5 ) * 0.001; 155 } 154 } 156 155 157 G4int n0Try = 0; 156 G4int n0Try = 0; 158 G4bool isThisOK = false; 157 G4bool isThisOK = false; 159 while ( n0Try < maxTrial ) // Loop checking 158 while ( n0Try < maxTrial ) // Loop checking, 11.03.2015, T. Koi 160 { 159 { 161 n0Try++; 160 n0Try++; 162 //std::cout << "TKDB packNucleons n0Try 161 //std::cout << "TKDB packNucleons n0Try " << n0Try << std::endl; 163 162 164 // Sampling Position 163 // Sampling Position 165 164 166 //std::cout << "TKDB Sampling Position " 165 //std::cout << "TKDB Sampling Position " << std::endl; 167 166 168 G4bool areThesePsOK = false; 167 G4bool areThesePsOK = false; 169 G4int npTry = 0; 168 G4int npTry = 0; 170 while ( npTry < maxTrial ) // Loop check 169 while ( npTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi 171 { 170 { 172 //std::cout << "TKDB Sampling Position n 171 //std::cout << "TKDB Sampling Position npTry " << npTry << std::endl; 173 npTry++; 172 npTry++; 174 G4int i = 0; 173 G4int i = 0; 175 if ( samplingPosition( i ) ) 174 if ( samplingPosition( i ) ) 176 { 175 { 177 //std::cout << "packNucleons sampl 176 //std::cout << "packNucleons samplingPosition 0 succeed " << std::endl; 178 for ( i = 1 ; i < GetMassNumber() 177 for ( i = 1 ; i < GetMassNumber() ; i++ ) 179 { 178 { 180 //std::cout << "packNucleons sa 179 //std::cout << "packNucleons samplingPosition " << i << " trying " << std::endl; 181 if ( !( samplingPosition( i ) ) 180 if ( !( samplingPosition( i ) ) ) 182 { 181 { 183 //std::cout << "packNucleons 182 //std::cout << "packNucleons samplingPosition " << i << " failed" << std::endl; 184 break; 183 break; 185 } 184 } 186 } 185 } 187 if ( i == GetMassNumber() ) 186 if ( i == GetMassNumber() ) 188 { 187 { 189 //std::cout << "packNucleons sa 188 //std::cout << "packNucleons samplingPosition all scucceed " << std::endl; 190 areThesePsOK = true; 189 areThesePsOK = true; 191 break; 190 break; 192 } 191 } 193 } 192 } 194 } 193 } 195 if ( areThesePsOK == false ) continue; 194 if ( areThesePsOK == false ) continue; 196 195 197 //std::cout << "TKDB Sampling Position E 196 //std::cout << "TKDB Sampling Position End" << std::endl; 198 197 199 // Calculate Two-body quantities 198 // Calculate Two-body quantities 200 199 201 meanfield->Cal2BodyQuantities(); 200 meanfield->Cal2BodyQuantities(); 202 std::vector< G4double > rho_a ( GetMassN 201 std::vector< G4double > rho_a ( GetMassNumber() , 0.0 ); 203 std::vector< G4double > rho_s ( GetMassN 202 std::vector< G4double > rho_s ( GetMassNumber() , 0.0 ); 204 std::vector< G4double > rho_c ( GetMassN 203 std::vector< G4double > rho_c ( GetMassNumber() , 0.0 ); 205 204 206 rho_l.resize ( GetMassNumber() , 0.0 ); 205 rho_l.resize ( GetMassNumber() , 0.0 ); 207 d_pot.resize ( GetMassNumber() , 0.0 ); 206 d_pot.resize ( GetMassNumber() , 0.0 ); 208 207 209 for ( G4int i = 0 ; i < GetMassNumber() 208 for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) 210 { 209 { 211 for ( G4int j = 0 ; j < GetMassNumber 210 for ( G4int j = 0 ; j < GetMassNumber() ; j++ ) 212 { 211 { 213 212 214 rho_a[ i ] += meanfield->GetRHA( j 213 rho_a[ i ] += meanfield->GetRHA( j , i ); 215 G4int k = 0; 214 G4int k = 0; 216 if ( participants[j]->GetDefinitio 215 if ( participants[j]->GetDefinition() != participants[i]->GetDefinition() ) 217 { 216 { 218 k = 1; 217 k = 1; 219 } 218 } 220 219 221 rho_s[ i ] += meanfield->GetRHA( j 220 rho_s[ i ] += meanfield->GetRHA( j , i )*( 1.0 - 2.0 * k ); // OK? 222 221 223 rho_c[ i ] += meanfield->GetRHE( j 222 rho_c[ i ] += meanfield->GetRHE( j , i ); 224 } 223 } 225 224 226 } 225 } 227 226 228 for ( G4int i = 0 ; i < GetMassNumber() 227 for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) 229 { 228 { 230 rho_l[i] = cdp * ( rho_a[i] + 1.0 ); 229 rho_l[i] = cdp * ( rho_a[i] + 1.0 ); 231 d_pot[i] = c0p * rho_a[i] 230 d_pot[i] = c0p * rho_a[i] 232 + c3p * G4Pow::GetInstance() 231 + c3p * G4Pow::GetInstance()->powA ( rho_a[i] , gamm ) 233 + csp * rho_s[i] 232 + csp * rho_s[i] 234 + clp * rho_c[i]; 233 + clp * rho_c[i]; 235 234 236 //std::cout << "d_pot[i] " << i << " 235 //std::cout << "d_pot[i] " << i << " " << d_pot[i] << std::endl; 237 } 236 } 238 237 239 238 240 // Sampling Momentum 239 // Sampling Momentum 241 240 242 //std::cout << "TKDB Sampling Momentum " 241 //std::cout << "TKDB Sampling Momentum " << std::endl; 243 242 244 phase_g.clear(); 243 phase_g.clear(); 245 phase_g.resize( GetMassNumber() , 0.0 ); 244 phase_g.resize( GetMassNumber() , 0.0 ); 246 245 247 //std::cout << "TKDB Sampling Momentum 1 246 //std::cout << "TKDB Sampling Momentum 1st " << std::endl; 248 G4bool isThis1stMOK = false; 247 G4bool isThis1stMOK = false; 249 G4int nmTry = 0; 248 G4int nmTry = 0; 250 while ( nmTry < maxTrial ) // Loop check 249 while ( nmTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi 251 { 250 { 252 nmTry++; 251 nmTry++; 253 G4int i = 0; 252 G4int i = 0; 254 if ( samplingMomentum( i ) ) 253 if ( samplingMomentum( i ) ) 255 { 254 { 256 isThis1stMOK = true; 255 isThis1stMOK = true; 257 break; 256 break; 258 } 257 } 259 } 258 } 260 if ( isThis1stMOK == false ) continue; 259 if ( isThis1stMOK == false ) continue; 261 260 262 //std::cout << "TKDB Sampling Momentum 2 261 //std::cout << "TKDB Sampling Momentum 2nd so on" << std::endl; 263 262 264 G4bool areTheseMsOK = true; 263 G4bool areTheseMsOK = true; 265 nmTry = 0; 264 nmTry = 0; 266 while ( nmTry < maxTrial ) // Loop check 265 while ( nmTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi 267 { 266 { 268 nmTry++; 267 nmTry++; 269 G4int i = 0; 268 G4int i = 0; 270 for ( i = 1 ; i < GetMassNumber() 269 for ( i = 1 ; i < GetMassNumber() ; i++ ) 271 { 270 { 272 //std::cout << "TKDB packNucleo 271 //std::cout << "TKDB packNucleons samplingMomentum try " << i << std::endl; 273 if ( !( samplingMomentum( i ) ) 272 if ( !( samplingMomentum( i ) ) ) 274 { 273 { 275 //std::cout << "TKDB packNuc 274 //std::cout << "TKDB packNucleons samplingMomentum " << i << " failed" << std::endl; 276 areTheseMsOK = false; 275 areTheseMsOK = false; 277 break; 276 break; 278 } 277 } 279 } 278 } 280 if ( i == GetMassNumber() ) 279 if ( i == GetMassNumber() ) 281 { 280 { 282 areTheseMsOK = true; 281 areTheseMsOK = true; 283 } 282 } 284 283 285 break; 284 break; 286 } 285 } 287 if ( areTheseMsOK == false ) continue; 286 if ( areTheseMsOK == false ) continue; 288 287 289 // Kill Angluar Momentum 288 // Kill Angluar Momentum 290 289 291 //std::cout << "TKDB Sampling Kill Anglu 290 //std::cout << "TKDB Sampling Kill Angluar Momentum " << std::endl; 292 291 293 killCMMotionAndAngularM(); 292 killCMMotionAndAngularM(); 294 293 295 294 296 // Check Binding Energy 295 // Check Binding Energy 297 296 298 //std::cout << "packNucleons Check Bindi 297 //std::cout << "packNucleons Check Binding Energy Begin " << std::endl; 299 298 300 G4double ekinal = 0.0; 299 G4double ekinal = 0.0; 301 for ( int i = 0 ; i < GetMassNumber() ; 300 for ( int i = 0 ; i < GetMassNumber() ; i++ ) 302 { 301 { 303 ekinal += participants[i]->GetKinetic 302 ekinal += participants[i]->GetKineticEnergy(); 304 } 303 } 305 304 306 meanfield->Cal2BodyQuantities(); 305 meanfield->Cal2BodyQuantities(); 307 306 308 G4double totalPotentialE = meanfield->Ge 307 G4double totalPotentialE = meanfield->GetTotalPotential(); 309 G4double ebinal = ( totalPotentialE + ek 308 G4double ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() ); 310 309 311 //std::cout << "packNucleons totalPotent 310 //std::cout << "packNucleons totalPotentialE " << totalPotentialE << std::endl; 312 //std::cout << "packNucleons ebinal " << 311 //std::cout << "packNucleons ebinal " << ebinal << std::endl; 313 //std::cout << "packNucleons ekinal " << 312 //std::cout << "packNucleons ekinal " << ekinal << std::endl; 314 313 315 if ( ebinal < ebin0 || ebinal > ebin1 ) 314 if ( ebinal < ebin0 || ebinal > ebin1 ) 316 { 315 { 317 //std::cout << "packNucleons ebin0 " 316 //std::cout << "packNucleons ebin0 " << ebin0 << std::endl; 318 //std::cout << "packNucleons ebin1 " 317 //std::cout << "packNucleons ebin1 " << ebin1 << std::endl; 319 //std::cout << "packNucleons ebinal " 318 //std::cout << "packNucleons ebinal " << ebinal << std::endl; 320 //std::cout << "packNucleons Check Bindi 319 //std::cout << "packNucleons Check Binding Energy Failed " << std::endl; 321 continue; 320 continue; 322 } 321 } 323 322 324 //std::cout << "packNucleons Check Bindi 323 //std::cout << "packNucleons Check Binding Energy End = OK " << std::endl; 325 324 326 325 327 // Energy Adujstment 326 // Energy Adujstment 328 327 329 G4double dtc = 1.0; 328 G4double dtc = 1.0; 330 G4double frg = -0.1; 329 G4double frg = -0.1; 331 G4double rdf0 = 0.5; 330 G4double rdf0 = 0.5; 332 331 333 G4double edif0 = ebinal - ebin00; 332 G4double edif0 = ebinal - ebin00; 334 333 335 G4double cfrc = 0.0; 334 G4double cfrc = 0.0; 336 if ( 0 < edif0 ) 335 if ( 0 < edif0 ) 337 { 336 { 338 cfrc = frg; 337 cfrc = frg; 339 } 338 } 340 else 339 else 341 { 340 { 342 cfrc = -frg; 341 cfrc = -frg; 343 } 342 } 344 343 345 G4int ifrc = 1; 344 G4int ifrc = 1; 346 345 347 G4int neaTry = 0; 346 G4int neaTry = 0; 348 347 349 G4bool isThisEAOK = false; 348 G4bool isThisEAOK = false; 350 while ( neaTry < maxTrial ) // Loop che 349 while ( neaTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi 351 { 350 { 352 neaTry++; 351 neaTry++; 353 352 354 G4double edif = ebinal - ebin00; 353 G4double edif = ebinal - ebin00; 355 354 356 //std::cout << "TKDB edif " << edif < 355 //std::cout << "TKDB edif " << edif << std::endl; 357 if ( std::abs ( edif ) < epse ) 356 if ( std::abs ( edif ) < epse ) 358 { 357 { 359 358 360 isThisEAOK = true; 359 isThisEAOK = true; 361 //std::cout << "isThisEAOK " << is 360 //std::cout << "isThisEAOK " << isThisEAOK << std::endl; 362 break; 361 break; 363 } 362 } 364 363 365 G4int jfrc = 0; 364 G4int jfrc = 0; 366 if ( edif < 0.0 ) 365 if ( edif < 0.0 ) 367 { 366 { 368 jfrc = 1; 367 jfrc = 1; 369 } 368 } 370 else 369 else 371 { 370 { 372 jfrc = -1; 371 jfrc = -1; 373 } 372 } 374 373 375 if ( jfrc != ifrc ) 374 if ( jfrc != ifrc ) 376 { 375 { 377 cfrc = -rdf0 * cfrc; 376 cfrc = -rdf0 * cfrc; 378 dtc = rdf0 * dtc; 377 dtc = rdf0 * dtc; 379 } 378 } 380 379 381 if ( jfrc == ifrc && std::abs( edif0 380 if ( jfrc == ifrc && std::abs( edif0 ) < std::abs( edif ) ) 382 { 381 { 383 cfrc = -rdf0 * cfrc; 382 cfrc = -rdf0 * cfrc; 384 dtc = rdf0 * dtc; 383 dtc = rdf0 * dtc; 385 } 384 } 386 385 387 ifrc = jfrc; 386 ifrc = jfrc; 388 edif0 = edif; 387 edif0 = edif; 389 388 390 meanfield->CalGraduate(); 389 meanfield->CalGraduate(); 391 390 392 for ( int i = 0 ; i < GetMassNumber() 391 for ( int i = 0 ; i < GetMassNumber() ; i++ ) 393 { 392 { 394 G4ThreeVector ri = participants[i] 393 G4ThreeVector ri = participants[i]->GetPosition(); 395 G4ThreeVector p_i = participants[i 394 G4ThreeVector p_i = participants[i]->GetMomentum(); 396 395 397 ri += dtc * ( meanfield->GetFFr(i) 396 ri += dtc * ( meanfield->GetFFr(i) - cfrc * ( meanfield->GetFFp(i) ) ); 398 p_i += dtc * ( meanfield->GetFFp(i 397 p_i += dtc * ( meanfield->GetFFp(i) + cfrc * ( meanfield->GetFFr(i) ) ); 399 398 400 participants[i]->SetPosition( ri ) 399 participants[i]->SetPosition( ri ); 401 participants[i]->SetMomentum( p_i 400 participants[i]->SetMomentum( p_i ); 402 } 401 } 403 402 404 ekinal = 0.0; 403 ekinal = 0.0; 405 404 406 for ( int i = 0 ; i < GetMassNumber() 405 for ( int i = 0 ; i < GetMassNumber() ; i++ ) 407 { 406 { 408 ekinal += participants[i]->GetKine 407 ekinal += participants[i]->GetKineticEnergy(); 409 } 408 } 410 409 411 meanfield->Cal2BodyQuantities(); 410 meanfield->Cal2BodyQuantities(); 412 totalPotentialE = meanfield->GetTotal 411 totalPotentialE = meanfield->GetTotalPotential(); 413 412 414 ebinal = ( totalPotentialE + ekinal ) 413 ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() ); 415 414 416 } 415 } 417 //std::cout << "isThisEAOK " << isThisEA 416 //std::cout << "isThisEAOK " << isThisEAOK << std::endl; 418 if ( isThisEAOK == false ) continue; 417 if ( isThisEAOK == false ) continue; 419 418 420 isThisOK = true; 419 isThisOK = true; 421 //std::cout << "isThisOK " << isThisOK < 420 //std::cout << "isThisOK " << isThisOK << std::endl; 422 break; 421 break; 423 422 424 } 423 } 425 424 426 if ( isThisOK == false ) 425 if ( isThisOK == false ) 427 { 426 { 428 G4cout << "GroundStateNucleus state cann 427 G4cout << "GroundStateNucleus state cannot be created. Try again with another parameters." << G4endl; 429 } 428 } 430 429 431 //std::cout << "packNucleons End" << std::e 430 //std::cout << "packNucleons End" << std::endl; 432 return; 431 return; 433 } 432 } 434 433 435 434 436 G4bool G4QMDGroundStateNucleus::samplingPositi 435 G4bool G4QMDGroundStateNucleus::samplingPosition( G4int i ) 437 { 436 { 438 437 439 G4bool result = false; 438 G4bool result = false; 440 439 441 G4int nTry = 0; 440 G4int nTry = 0; 442 441 443 while ( nTry < maxTrial ) // Loop checking 442 while ( nTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi 444 { 443 { 445 444 446 //std::cout << "samplingPosition i th pa 445 //std::cout << "samplingPosition i th particle, nTtry " << i << " " << nTry << std::endl; 447 446 448 G4double rwod = -1.0; 447 G4double rwod = -1.0; 449 G4double rrr = 0.0; 448 G4double rrr = 0.0; 450 449 451 G4double rx = 0.0; 450 G4double rx = 0.0; 452 G4double ry = 0.0; 451 G4double ry = 0.0; 453 G4double rz = 0.0; 452 G4double rz = 0.0; 454 453 455 G4int icounter = 0; 454 G4int icounter = 0; 456 G4int icounter_max = 1024; 455 G4int icounter_max = 1024; 457 while ( G4UniformRand() * rmax > rwod ) 456 while ( G4UniformRand() * rmax > rwod ) // Loop checking, 11.03.2015, T. Koi 458 { 457 { 459 icounter++; 458 icounter++; 460 if ( icounter > icounter_max ) { 459 if ( icounter > icounter_max ) { 461 G4cout << "Loop-counter exceeded the thr 460 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl; 462 break; 461 break; 463 } 462 } 464 463 465 G4double rsqr = 10.0; 464 G4double rsqr = 10.0; 466 G4int jcounter = 0; 465 G4int jcounter = 0; 467 G4int jcounter_max = 1024; 466 G4int jcounter_max = 1024; 468 while ( rsqr > 1.0 ) // Loop checking 467 while ( rsqr > 1.0 ) // Loop checking, 11.03.2015, T. Koi 469 { 468 { 470 jcounter++; 469 jcounter++; 471 if ( jcounter > jcounter_max ) { 470 if ( jcounter > jcounter_max ) { 472 G4cout << "Loop-counter exceeded the 471 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl; 473 break; 472 break; 474 } 473 } 475 rx = 1.0 - 2.0 * G4UniformRand(); 474 rx = 1.0 - 2.0 * G4UniformRand(); 476 ry = 1.0 - 2.0 * G4UniformRand(); 475 ry = 1.0 - 2.0 * G4UniformRand(); 477 rz = 1.0 - 2.0 * G4UniformRand(); 476 rz = 1.0 - 2.0 * G4UniformRand(); 478 rsqr = rx*rx + ry*ry + rz*rz; 477 rsqr = rx*rx + ry*ry + rz*rz; 479 } 478 } 480 rrr = radm * std::sqrt ( rsqr ); 479 rrr = radm * std::sqrt ( rsqr ); 481 rwod = 1.0 / ( 1.0 + G4Exp ( ( rrr - 480 rwod = 1.0 / ( 1.0 + G4Exp ( ( rrr - rt00 ) / saa ) ); 482 481 483 } 482 } 484 483 485 participants[i]->SetPosition( G4ThreeVec 484 participants[i]->SetPosition( G4ThreeVector( rx , ry , rz )*radm ); 486 485 487 if ( i == 0 ) 486 if ( i == 0 ) 488 { 487 { 489 result = true; 488 result = true; 490 return result; 489 return result; 491 } 490 } 492 491 493 // i > 1 ( Second Particle or later ) 492 // i > 1 ( Second Particle or later ) 494 // Check Distance to others 493 // Check Distance to others 495 494 496 G4bool isThisOK = true; 495 G4bool isThisOK = true; 497 for ( G4int j = 0 ; j < i ; j++ ) 496 for ( G4int j = 0 ; j < i ; j++ ) 498 { 497 { 499 498 500 G4double r2 = participants[j]->GetPo 499 G4double r2 = participants[j]->GetPosition().diff2( participants[i]->GetPosition() ); 501 G4double dmin2 = 0.0; 500 G4double dmin2 = 0.0; 502 501 503 if ( participants[j]->GetDefinition() 502 if ( participants[j]->GetDefinition() == participants[i]->GetDefinition() ) 504 { 503 { 505 dmin2 = dsam2; 504 dmin2 = dsam2; 506 } 505 } 507 else 506 else 508 { 507 { 509 dmin2 = ddif2; 508 dmin2 = ddif2; 510 } 509 } 511 510 512 //std::cout << "distance between j an 511 //std::cout << "distance between j and i " << j << " " << i << " " << r2 << " " << dmin2 << std::endl; 513 if ( r2 < dmin2 ) 512 if ( r2 < dmin2 ) 514 { 513 { 515 isThisOK = false; 514 isThisOK = false; 516 break; 515 break; 517 } 516 } 518 517 519 } 518 } 520 519 521 if ( isThisOK == true ) 520 if ( isThisOK == true ) 522 { 521 { 523 result = true; 522 result = true; 524 return result; 523 return result; 525 } 524 } 526 525 527 nTry++; 526 nTry++; 528 527 529 } 528 } 530 529 531 // Here return "false" 530 // Here return "false" 532 return result; 531 return result; 533 } 532 } 534 533 535 534 536 535 537 G4bool G4QMDGroundStateNucleus::samplingMoment 536 G4bool G4QMDGroundStateNucleus::samplingMomentum( G4int i ) 538 { 537 { 539 538 540 //std::cout << "TKDB samplingMomentum for " 539 //std::cout << "TKDB samplingMomentum for " << i << std::endl; 541 540 542 G4bool result = false; 541 G4bool result = false; 543 542 544 G4double pfm = hbc * G4Pow::GetInstance()-> 543 G4double pfm = hbc * G4Pow::GetInstance()->A13 ( ( 3.0 / 2.0 * pi*pi * rho_l[i] ) ); 545 544 546 if ( 10 < GetMassNumber() && -5.5 < ebini 545 if ( 10 < GetMassNumber() && -5.5 < ebini ) 547 { 546 { 548 pfm = pfm * ( 1.0 + 0.2 * std::sqrt( std 547 pfm = pfm * ( 1.0 + 0.2 * std::sqrt( std::abs( 8.0 + ebini ) / 8.0 ) ); 549 } 548 } 550 549 551 //std::cout << "TKDB samplingMomentum pfm " 550 //std::cout << "TKDB samplingMomentum pfm " << pfm << std::endl; 552 551 553 std::vector< G4double > phase; 552 std::vector< G4double > phase; 554 phase.resize( i+1 ); // i start from 0 553 phase.resize( i+1 ); // i start from 0 555 554 556 G4int ntry = 0; 555 G4int ntry = 0; 557 // 710 556 // 710 558 while ( ntry < maxTrial ) // Loop checking 557 while ( ntry < maxTrial ) // Loop checking, 11.03.2015, T. Koi 559 { 558 { 560 559 561 //std::cout << " TKDB ntry " << ntry << 560 //std::cout << " TKDB ntry " << ntry << std::endl; 562 ntry++; 561 ntry++; 563 562 564 G4double ke = DBL_MAX; 563 G4double ke = DBL_MAX; 565 564 566 G4int tkdb_i =0; 565 G4int tkdb_i =0; 567 // 700 566 // 700 568 G4int icounter = 0; 567 G4int icounter = 0; 569 G4int icounter_max = 1024; 568 G4int icounter_max = 1024; 570 while ( ke + d_pot [i] > edepth ) // Loo 569 while ( ke + d_pot [i] > edepth ) // Loop checking, 11.03.2015, T. Koi 571 { 570 { 572 icounter++; 571 icounter++; 573 if ( icounter > icounter_max ) { 572 if ( icounter > icounter_max ) { 574 G4cout << "Loop-counter exceeded the thr 573 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl; 575 break; 574 break; 576 } 575 } 577 576 578 G4double psqr = 10.0; 577 G4double psqr = 10.0; 579 G4double px = 0.0; 578 G4double px = 0.0; 580 G4double py = 0.0; 579 G4double py = 0.0; 581 G4double pz = 0.0; 580 G4double pz = 0.0; 582 581 583 G4int jcounter = 0; 582 G4int jcounter = 0; 584 G4int jcounter_max = 1024; 583 G4int jcounter_max = 1024; 585 while ( psqr > 1.0 ) // Loop checking 584 while ( psqr > 1.0 ) // Loop checking, 11.03.2015, T. Koi 586 { 585 { 587 jcounter++; 586 jcounter++; 588 if ( jcounter > jcounter_max ) { 587 if ( jcounter > jcounter_max ) { 589 G4cout << "Loop-counter exceeded the 588 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl; 590 break; 589 break; 591 } 590 } 592 px = 1.0 - 2.0*G4UniformRand(); 591 px = 1.0 - 2.0*G4UniformRand(); 593 py = 1.0 - 2.0*G4UniformRand(); 592 py = 1.0 - 2.0*G4UniformRand(); 594 pz = 1.0 - 2.0*G4UniformRand(); 593 pz = 1.0 - 2.0*G4UniformRand(); 595 594 596 psqr = px*px + py*py + pz*pz; 595 psqr = px*px + py*py + pz*pz; 597 } 596 } 598 597 599 G4ThreeVector p ( px , py , pz ); 598 G4ThreeVector p ( px , py , pz ); 600 p = pfm * p; 599 p = pfm * p; 601 participants[i]->SetMomentum( p ); 600 participants[i]->SetMomentum( p ); 602 G4LorentzVector p4 = participants[i]- 601 G4LorentzVector p4 = participants[i]->Get4Momentum(); 603 //ke = p4.e() - p4.restMass(); 602 //ke = p4.e() - p4.restMass(); 604 ke = participants[i]->GetKineticEnerg 603 ke = participants[i]->GetKineticEnergy(); 605 604 606 605 607 tkdb_i++; 606 tkdb_i++; 608 if ( tkdb_i > maxTrial ) return resul 607 if ( tkdb_i > maxTrial ) return result; // return false 609 608 610 } 609 } 611 610 612 //std::cout << "TKDB ke d_pot[i] " << ke 611 //std::cout << "TKDB ke d_pot[i] " << ke << " " << d_pot[i] << std::endl; 613 612 614 613 615 if ( i == 0 ) 614 if ( i == 0 ) 616 { 615 { 617 result = true; 616 result = true; 618 return result; 617 return result; 619 } 618 } 620 619 621 G4bool isThisOK = true; 620 G4bool isThisOK = true; 622 621 623 // Check Pauli principle 622 // Check Pauli principle 624 623 625 phase[ i ] = 0.0; 624 phase[ i ] = 0.0; 626 625 627 //std::cout << "TKDB Check Pauli Princip 626 //std::cout << "TKDB Check Pauli Principle " << i << std::endl; 628 627 629 for ( G4int j = 0 ; j < i ; j++ ) 628 for ( G4int j = 0 ; j < i ; j++ ) 630 { 629 { 631 phase[ j ] = 0.0; 630 phase[ j ] = 0.0; 632 //std::cout << "TKDB Check Pauli Prin 631 //std::cout << "TKDB Check Pauli Principle i , j " << i << " , " << j << std::endl; 633 G4double expa = 0.0; 632 G4double expa = 0.0; 634 if ( participants[j]->GetDefinition() 633 if ( participants[j]->GetDefinition() == participants[i]->GetDefinition() ) 635 { 634 { 636 635 637 expa = - meanfield->GetRR2(i,j) * 636 expa = - meanfield->GetRR2(i,j) * cpw; 638 637 639 if ( expa > epsx ) 638 if ( expa > epsx ) 640 { 639 { 641 G4ThreeVector p_i = participant 640 G4ThreeVector p_i = participants[i]->GetMomentum(); 642 G4ThreeVector pj = participants 641 G4ThreeVector pj = participants[j]->GetMomentum(); 643 G4double dist2_p = p_i.diff2( p 642 G4double dist2_p = p_i.diff2( pj ); 644 643 645 dist2_p = dist2_p*cph; 644 dist2_p = dist2_p*cph; 646 expa = expa - dist2_p; 645 expa = expa - dist2_p; 647 646 648 if ( expa > epsx ) 647 if ( expa > epsx ) 649 { 648 { 650 649 651 phase[j] = G4Exp ( expa ); 650 phase[j] = G4Exp ( expa ); 652 651 653 if ( phase[j] * cpc > 0.2 ) 652 if ( phase[j] * cpc > 0.2 ) 654 { 653 { 655 /* 654 /* 656 G4cout << "TKDB Check Pauli Principle 655 G4cout << "TKDB Check Pauli Principle A i , j " << i << " , " << j << G4endl; 657 G4cout << "TKDB Check Pauli Principle 656 G4cout << "TKDB Check Pauli Principle phase[j] " << phase[j] << G4endl; 658 G4cout << "TKDB Check Pauli Principle 657 G4cout << "TKDB Check Pauli Principle phase[j]*cpc > 0.2 " << phase[j]*cpc << G4endl; 659 */ 658 */ 660 isThisOK = false; 659 isThisOK = false; 661 break; 660 break; 662 } 661 } 663 if ( ( phase_g[j] + phase[j] 662 if ( ( phase_g[j] + phase[j] ) * cpc > 0.5 ) 664 { 663 { 665 /* 664 /* 666 G4cout << "TKDB Check Pauli Principle 665 G4cout << "TKDB Check Pauli Principle B i , j " << i << " , " << j << G4endl; 667 G4cout << "TKDB Check Pauli Principle 666 G4cout << "TKDB Check Pauli Principle B phase_g[j] " << phase_g[j] << G4endl; 668 G4cout << "TKDB Check Pauli Principle 667 G4cout << "TKDB Check Pauli Principle B phase[j] " << phase[j] << G4endl; 669 G4cout << "TKDB Check Pauli Principle 668 G4cout << "TKDB Check Pauli Principle B phase_g[j] + phase[j] ) * cpc > 0.5 " << ( phase_g[j] + phase[j] ) * cpc << G4endl; 670 */ 669 */ 671 isThisOK = false; 670 isThisOK = false; 672 break; 671 break; 673 } 672 } 674 673 675 phase[i] += phase[j]; 674 phase[i] += phase[j]; 676 if ( phase[i] * cpc > 0.3 ) 675 if ( phase[i] * cpc > 0.3 ) 677 { 676 { 678 /* 677 /* 679 G4cout << "TKDB Check Pauli Principle 678 G4cout << "TKDB Check Pauli Principle C i , j " << i << " , " << j << G4endl; 680 G4cout << "TKDB Check Pauli Principle 679 G4cout << "TKDB Check Pauli Principle C phase[i] " << phase[i] << G4endl; 681 G4cout << "TKDB Check Pauli Principle 680 G4cout << "TKDB Check Pauli Principle C phase[i] * cpc > 0.3 " << phase[i] * cpc << G4endl; 682 */ 681 */ 683 isThisOK = false; 682 isThisOK = false; 684 break; 683 break; 685 } 684 } 686 685 687 //std::cout << "TKDB Check P 686 //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl; 688 687 689 } 688 } 690 else 689 else 691 { 690 { 692 //std::cout << "TKDB Check P 691 //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl; 693 } 692 } 694 693 695 } 694 } 696 else 695 else 697 { 696 { 698 //std::cout << "TKDB Check Paul 697 //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl; 699 } 698 } 700 699 701 } 700 } 702 else 701 else 703 { 702 { 704 //std::cout << "TKDB Check Pauli P 703 //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl; 705 } 704 } 706 705 707 } 706 } 708 707 709 if ( isThisOK == true ) 708 if ( isThisOK == true ) 710 { 709 { 711 710 712 phase_g[i] = phase[i]; 711 phase_g[i] = phase[i]; 713 712 714 for ( int j = 0 ; j < i ; j++ ) 713 for ( int j = 0 ; j < i ; j++ ) 715 { 714 { 716 phase_g[j] += phase[j]; 715 phase_g[j] += phase[j]; 717 } 716 } 718 717 719 result = true; 718 result = true; 720 return result; 719 return result; 721 } 720 } 722 721 723 } 722 } 724 723 725 return result; 724 return result; 726 725 727 } 726 } 728 727 729 728 730 729 731 void G4QMDGroundStateNucleus::killCMMotionAndA 730 void G4QMDGroundStateNucleus::killCMMotionAndAngularM() 732 { 731 { 733 732 734 // CalEnergyAndAngularMomentumInCM(); 733 // CalEnergyAndAngularMomentumInCM(); 735 734 736 //std::vector< G4ThreeVector > p ( GetMassN 735 //std::vector< G4ThreeVector > p ( GetMassNumber() , 0.0 ); 737 //std::vector< G4ThreeVector > r ( GetMassN 736 //std::vector< G4ThreeVector > r ( GetMassNumber() , 0.0 ); 738 737 739 // Move to cm system 738 // Move to cm system 740 739 741 G4ThreeVector pcm_tmp ( 0.0 ); 740 G4ThreeVector pcm_tmp ( 0.0 ); 742 G4ThreeVector rcm_tmp ( 0.0 ); 741 G4ThreeVector rcm_tmp ( 0.0 ); 743 G4double sumMass = 0.0; 742 G4double sumMass = 0.0; 744 743 745 for ( G4int i = 0 ; i < GetMassNumber() ; i 744 for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) 746 { 745 { 747 pcm_tmp += participants[i]->GetMomentum( 746 pcm_tmp += participants[i]->GetMomentum(); 748 rcm_tmp += participants[i]->GetPosition( 747 rcm_tmp += participants[i]->GetPosition() * participants[i]->GetMass(); 749 sumMass += participants[i]->GetMass(); 748 sumMass += participants[i]->GetMass(); 750 } 749 } 751 750 752 pcm_tmp = pcm_tmp/GetMassNumber(); 751 pcm_tmp = pcm_tmp/GetMassNumber(); 753 rcm_tmp = rcm_tmp/sumMass; 752 rcm_tmp = rcm_tmp/sumMass; 754 753 755 for ( G4int i = 0 ; i < GetMassNumber() ; i 754 for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) 756 { 755 { 757 participants[i]->SetMomentum( participan 756 participants[i]->SetMomentum( participants[i]->GetMomentum() - pcm_tmp ); 758 participants[i]->SetPosition( participan 757 participants[i]->SetPosition( participants[i]->GetPosition() - rcm_tmp ); 759 } 758 } 760 759 761 // kill the angular momentum 760 // kill the angular momentum 762 761 763 G4ThreeVector ll ( 0.0 ); 762 G4ThreeVector ll ( 0.0 ); 764 for ( G4int i = 0 ; i < GetMassNumber() ; i 763 for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) 765 { 764 { 766 ll += participants[i]->GetPosition().cro 765 ll += participants[i]->GetPosition().cross ( participants[i]->GetMomentum() ); 767 } 766 } 768 767 769 G4double rr[3][3]; 768 G4double rr[3][3]; 770 G4double ss[3][3]; 769 G4double ss[3][3]; 771 for ( G4int i = 0 ; i < 3 ; i++ ) 770 for ( G4int i = 0 ; i < 3 ; i++ ) 772 { 771 { 773 for ( G4int j = 0 ; j < 3 ; j++ ) 772 for ( G4int j = 0 ; j < 3 ; j++ ) 774 { 773 { 775 rr [i][j] = 0.0; 774 rr [i][j] = 0.0; 776 775 777 if ( i == j ) 776 if ( i == j ) 778 { 777 { 779 ss [i][j] = 1.0; 778 ss [i][j] = 1.0; 780 } 779 } 781 else 780 else 782 { 781 { 783 ss [i][j] = 0.0; 782 ss [i][j] = 0.0; 784 } 783 } 785 } 784 } 786 } 785 } 787 786 788 for ( G4int i = 0 ; i < GetMassNumber() ; i 787 for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) 789 { 788 { 790 G4ThreeVector r = participants[i]->GetPo 789 G4ThreeVector r = participants[i]->GetPosition(); 791 rr[0][0] += r.y() * r.y() + r.z() * r.z( 790 rr[0][0] += r.y() * r.y() + r.z() * r.z(); 792 rr[1][0] += - r.y() * r.x(); 791 rr[1][0] += - r.y() * r.x(); 793 rr[2][0] += - r.z() * r.x(); 792 rr[2][0] += - r.z() * r.x(); 794 rr[0][1] += - r.x() * r.y(); 793 rr[0][1] += - r.x() * r.y(); 795 rr[1][1] += r.z() * r.z() + r.x() * r.x( 794 rr[1][1] += r.z() * r.z() + r.x() * r.x(); 796 rr[2][1] += - r.z() * r.y(); 795 rr[2][1] += - r.z() * r.y(); 797 rr[2][0] += - r.x() * r.z(); 796 rr[2][0] += - r.x() * r.z(); 798 rr[2][1] += - r.y() * r.z(); 797 rr[2][1] += - r.y() * r.z(); 799 rr[2][2] += r.x() * r.x() + r.y() * r.y( 798 rr[2][2] += r.x() * r.x() + r.y() * r.y(); 800 } 799 } 801 800 802 for ( G4int i = 0 ; i < 3 ; i++ ) 801 for ( G4int i = 0 ; i < 3 ; i++ ) 803 { 802 { 804 G4double x = rr [i][i]; 803 G4double x = rr [i][i]; 805 for ( G4int j = 0 ; j < 3 ; j++ ) 804 for ( G4int j = 0 ; j < 3 ; j++ ) 806 { 805 { 807 rr[i][j] = rr[i][j] / x; 806 rr[i][j] = rr[i][j] / x; 808 ss[i][j] = ss[i][j] / x; 807 ss[i][j] = ss[i][j] / x; 809 } 808 } 810 for ( G4int j = 0 ; j < 3 ; j++ ) 809 for ( G4int j = 0 ; j < 3 ; j++ ) 811 { 810 { 812 if ( j != i ) 811 if ( j != i ) 813 { 812 { 814 G4double y = rr [j][i]; 813 G4double y = rr [j][i]; 815 for ( G4int k = 0 ; k < 3 ; k++ ) 814 for ( G4int k = 0 ; k < 3 ; k++ ) 816 { 815 { 817 rr[j][k] += -y * rr[i][k]; 816 rr[j][k] += -y * rr[i][k]; 818 ss[j][k] += -y * ss[i][k]; 817 ss[j][k] += -y * ss[i][k]; 819 } 818 } 820 } 819 } 821 } 820 } 822 } 821 } 823 822 824 G4double opl[3]; 823 G4double opl[3]; 825 G4double rll[3]; 824 G4double rll[3]; 826 825 827 rll[0] = ll.x(); 826 rll[0] = ll.x(); 828 rll[1] = ll.y(); 827 rll[1] = ll.y(); 829 rll[2] = ll.z(); 828 rll[2] = ll.z(); 830 829 831 for ( G4int i = 0 ; i < 3 ; i++ ) 830 for ( G4int i = 0 ; i < 3 ; i++ ) 832 { 831 { 833 opl[i] = 0.0; 832 opl[i] = 0.0; 834 833 835 for ( G4int j = 0 ; j < 3 ; j++ ) 834 for ( G4int j = 0 ; j < 3 ; j++ ) 836 { 835 { 837 opl[i] += ss[i][j]*rll[j]; 836 opl[i] += ss[i][j]*rll[j]; 838 } 837 } 839 } 838 } 840 839 841 for ( G4int i = 0 ; i < GetMassNumber() ; i 840 for ( G4int i = 0 ; i < GetMassNumber() ; i++ ) 842 { 841 { 843 G4ThreeVector p_i = participants[i]->Get 842 G4ThreeVector p_i = participants[i]->GetMomentum() ; 844 G4ThreeVector ri = participants[i]->GetP 843 G4ThreeVector ri = participants[i]->GetPosition() ; 845 G4ThreeVector opl_v ( opl[0] , opl[1] , 844 G4ThreeVector opl_v ( opl[0] , opl[1] , opl[2] ); 846 845 847 p_i += ri.cross(opl_v); 846 p_i += ri.cross(opl_v); 848 847 849 participants[i]->SetMomentum( p_i ); 848 participants[i]->SetMomentum( p_i ); 850 } 849 } 851 850 852 } 851 } 853 852