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