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