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