| 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 // 081120 Add Update by T. Koi 26 // 081120 Add Update by T. Koi
27 // 27 //
28 // 230307 Skyrme-QMD parameters added by Y-H. 28 // 230307 Skyrme-QMD parameters added by Y-H. Sato and A. Haga
29 // 230307 "CalDensityProfile" and "CalChargeDe 29 // 230307 "CalDensityProfile" and "CalChargeDensityProfile" functions added by Y-H. Sato and A. Haga
30 // 230307 "GetSingleEnergy" and "GetTotalEnerg 30 // 230307 "GetSingleEnergy" and "GetTotalEnergy" functions added by Y-H. Sato and A. Haga
31 31
32 #include <map> 32 #include <map>
33 #include <algorithm> 33 #include <algorithm>
34 #include <numeric> 34 #include <numeric>
35 35
36 #include <cmath> 36 #include <cmath>
37 #include <CLHEP/Random/Stat.h> 37 #include <CLHEP/Random/Stat.h>
38 38
39 #include "G4LightIonQMDMeanField.hh" 39 #include "G4LightIonQMDMeanField.hh"
40 #include "G4LightIonQMDParameters.hh" 40 #include "G4LightIonQMDParameters.hh"
41 #include "G4Exp.hh" 41 #include "G4Exp.hh"
42 #include "G4Pow.hh" 42 #include "G4Pow.hh"
43 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh"
44 #include "Randomize.hh" 44 #include "Randomize.hh"
45 45
46 G4LightIonQMDMeanField::G4LightIonQMDMeanField 46 G4LightIonQMDMeanField::G4LightIonQMDMeanField()
47 { 47 {
48 G4LightIonQMDParameters* parameters = G4Lig 48 G4LightIonQMDParameters* parameters = G4LightIonQMDParameters::GetInstance();
49 wl = parameters->Get_wl(); 49 wl = parameters->Get_wl();
50 cl = parameters->Get_cl(); 50 cl = parameters->Get_cl();
51 rho0 = parameters->Get_rho0(); 51 rho0 = parameters->Get_rho0();
52 hbc = parameters->Get_hbc(); 52 hbc = parameters->Get_hbc();
53 gamm = parameters->Get_gamm(); 53 gamm = parameters->Get_gamm();
54 eta = parameters->Get_eta(); // Skyrme-QMD 54 eta = parameters->Get_eta(); // Skyrme-QMD
55 kappas = parameters->Get_kappas(); // Skyrm 55 kappas = parameters->Get_kappas(); // Skyrme-QMD
56 56
57 cpw = parameters->Get_cpw(); 57 cpw = parameters->Get_cpw();
58 cph = parameters->Get_cph(); 58 cph = parameters->Get_cph();
59 cpc = parameters->Get_cpc(); 59 cpc = parameters->Get_cpc();
60 60
61 c0 = parameters->Get_c0(); 61 c0 = parameters->Get_c0();
62 c3 = parameters->Get_c3(); 62 c3 = parameters->Get_c3();
63 cs = parameters->Get_cs(); 63 cs = parameters->Get_cs();
64 g0 = parameters->Get_g0(); // Skyrme-QMD 64 g0 = parameters->Get_g0(); // Skyrme-QMD
65 g0iso = parameters->Get_g0iso(); // Skyrme- 65 g0iso = parameters->Get_g0iso(); // Skyrme-QMD
66 gtau0 = parameters->Get_gtau0(); // Skyrme- 66 gtau0 = parameters->Get_gtau0(); // Skyrme-QMD
67 67
68 // distance 68 // distance
69 c0w = 1.0/4.0/wl; 69 c0w = 1.0/4.0/wl;
70 c0sw = std::sqrt( c0w ); 70 c0sw = std::sqrt( c0w );
71 clw = 2.0 / std::sqrt ( 4.0 * pi * wl ); 71 clw = 2.0 / std::sqrt ( 4.0 * pi * wl );
72 72
73 // graduate 73 // graduate
74 c0g = - c0 / ( 2.0 * wl ); 74 c0g = - c0 / ( 2.0 * wl );
75 c3g = - c3 / ( 4.0 * wl ) * gamm; 75 c3g = - c3 / ( 4.0 * wl ) * gamm;
76 csg = - cs / ( 2.0 * wl ); 76 csg = - cs / ( 2.0 * wl );
77 pag = gamm - 1; 77 pag = gamm - 1;
78 pag_tau = eta - 1; // Skyrme-QMD 78 pag_tau = eta - 1; // Skyrme-QMD
79 cg0 = - g0 / ( 2.0 * wl ); // Skyrme-QMD 79 cg0 = - g0 / ( 2.0 * wl ); // Skyrme-QMD
80 cgtau0 = - gtau0 / ( 4.0 * wl ) * eta; // 80 cgtau0 = - gtau0 / ( 4.0 * wl ) * eta; // Skyrme-QMD
81 81
82 system = nullptr; // will be set through Se 82 system = nullptr; // will be set through SetSystem method
83 } 83 }
84 84
85 void G4LightIonQMDMeanField::SetSystem ( G4QMD 85 void G4LightIonQMDMeanField::SetSystem ( G4QMDSystem* aSystem )
86 { 86 {
87 system = aSystem; 87 system = aSystem;
88 88
89 G4int n = system->GetTotalNumberOfParticipa 89 G4int n = system->GetTotalNumberOfParticipant();
90 90
91 pp2.clear(); 91 pp2.clear();
92 rr2.clear(); 92 rr2.clear();
93 rbij.clear(); 93 rbij.clear();
94 rha.clear(); 94 rha.clear();
95 rhe.clear(); 95 rhe.clear();
96 rhc.clear(); 96 rhc.clear();
97 97
98 rr2.resize( n ); 98 rr2.resize( n );
99 pp2.resize( n ); 99 pp2.resize( n );
100 rbij.resize( n ); 100 rbij.resize( n );
101 rha.resize( n ); 101 rha.resize( n );
102 rhe.resize( n ); 102 rhe.resize( n );
103 rhc.resize( n ); 103 rhc.resize( n );
104 104
105 for ( G4int i = 0 ; i < n ; ++i ) 105 for ( G4int i = 0 ; i < n ; ++i )
106 { 106 {
107 rr2[i].resize( n ); 107 rr2[i].resize( n );
108 pp2[i].resize( n ); 108 pp2[i].resize( n );
109 rbij[i].resize( n ); 109 rbij[i].resize( n );
110 rha[i].resize( n ); 110 rha[i].resize( n );
111 rhe[i].resize( n ); 111 rhe[i].resize( n );
112 rhc[i].resize( n ); 112 rhc[i].resize( n );
113 } 113 }
114 114
115 ffr.clear(); 115 ffr.clear();
116 ffp.clear(); 116 ffp.clear();
117 rh3d.clear(); 117 rh3d.clear();
118 rh3d_tau.clear(); // Skyrme-QMD 118 rh3d_tau.clear(); // Skyrme-QMD
119 119
120 ffr.resize( n ); 120 ffr.resize( n );
121 ffp.resize( n ); 121 ffp.resize( n );
122 rh3d.resize( n ); 122 rh3d.resize( n );
123 rh3d_tau.resize( n ); // Skyrme-QMD 123 rh3d_tau.resize( n ); // Skyrme-QMD
124 124
125 Cal2BodyQuantities(); 125 Cal2BodyQuantities();
126 } 126 }
127 127
128 void G4LightIonQMDMeanField::SetNucleus ( G4Li 128 void G4LightIonQMDMeanField::SetNucleus ( G4LightIonQMDNucleus* aNucleus )
129 { 129 {
130 SetSystem( aNucleus ); 130 SetSystem( aNucleus );
131 131
132 G4double totalPotential = GetTotalPotential 132 G4double totalPotential = GetTotalPotential();
133 aNucleus->SetTotalPotential( totalPotential 133 aNucleus->SetTotalPotential( totalPotential );
134 aNucleus->CalEnergyAndAngularMomentumInCM() 134 aNucleus->CalEnergyAndAngularMomentumInCM();
135 } 135 }
136 136
137 void G4LightIonQMDMeanField::Cal2BodyQuantitie 137 void G4LightIonQMDMeanField::Cal2BodyQuantities()
138 { 138 {
139 if ( system->GetTotalNumberOfParticipant() 139 if ( system->GetTotalNumberOfParticipant() < 2 ) { return; }
140 140
141 for ( G4int j = 1 ; j < system->GetTotalNum 141 for ( G4int j = 1 ; j < system->GetTotalNumberOfParticipant() ; ++j )
142 { 142 {
143 G4ThreeVector rj = system->GetParticipan 143 G4ThreeVector rj = system->GetParticipant( j )->GetPosition();
144 G4LorentzVector p4j = system->GetPartici 144 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
145 145
146 for ( G4int i = 0 ; i < j ; ++i ) 146 for ( G4int i = 0 ; i < j ; ++i )
147 { 147 {
148 G4ThreeVector ri = system->GetPartici 148 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
149 G4LorentzVector p4i = system->GetPart 149 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum();
150 150
151 G4ThreeVector rij = ri - rj; 151 G4ThreeVector rij = ri - rj;
152 G4ThreeVector pij = (p4i - p4j).v(); 152 G4ThreeVector pij = (p4i - p4j).v();
153 G4LorentzVector p4ij = p4i - p4j; 153 G4LorentzVector p4ij = p4i - p4j;
154 G4ThreeVector bij = ( p4i + p4j ).boo 154 G4ThreeVector bij = ( p4i + p4j ).boostVector();
155 G4double gammaij = ( p4i + p4j ).gamm 155 G4double gammaij = ( p4i + p4j ).gamma();
156 156
157 G4double eij = ( p4i + p4j ).e(); 157 G4double eij = ( p4i + p4j ).e();
158 158
159 G4double rbrb = rij*bij; 159 G4double rbrb = rij*bij;
160 G4double rij2 = rij*rij; 160 G4double rij2 = rij*rij;
161 G4double pij2 = pij*pij; 161 G4double pij2 = pij*pij;
162 162
163 rbrb = irelcr * rbrb; 163 rbrb = irelcr * rbrb;
164 G4double gamma2_ij = gammaij*gammaij 164 G4double gamma2_ij = gammaij*gammaij;
165 165
166 rr2[i][j] = rij2 + gamma2_ij * rbrb*r 166 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb;
167 rr2[j][i] = rr2[i][j]; 167 rr2[j][i] = rr2[i][j];
168 168
169 rbij[i][j] = gamma2_ij * rbrb; 169 rbij[i][j] = gamma2_ij * rbrb;
170 rbij[j][i] = - rbij[i][j]; 170 rbij[j][i] = - rbij[i][j];
171 171
172 pp2[i][j] = pij2 172 pp2[i][j] = pij2
173 + irelcr * ( - G4Pow::GetIn 173 + irelcr * ( - G4Pow::GetInstance()->powN ( p4i.e() - p4j.e() , 2 )
174 + gamma2_ij * G4Pow::GetIns 174 + gamma2_ij * G4Pow::GetInstance()->powN ( ( ( p4i.m2() - p4j.m2() ) / eij ) , 2 ) );
175 175
176 176
177 pp2[j][i] = pp2[i][j]; 177 pp2[j][i] = pp2[i][j];
178 178
179 // Gauss term 179 // Gauss term
180 180
181 G4double expa1 = - rr2[i][j] * c0w; 181 G4double expa1 = - rr2[i][j] * c0w;
182 182
183 G4double rh1; 183 G4double rh1;
184 if ( expa1 > epsx ) 184 if ( expa1 > epsx )
185 { 185 {
186 rh1 = G4Exp( expa1 ); 186 rh1 = G4Exp( expa1 );
187 } 187 }
188 else 188 else
189 { 189 {
190 rh1 = 0.0; 190 rh1 = 0.0;
191 } 191 }
192 192
193 G4int ibry = system->GetParticipant(i 193 G4int ibry = system->GetParticipant(i)->GetBaryonNumber();
194 G4int jbry = system->GetParticipant(j 194 G4int jbry = system->GetParticipant(j)->GetBaryonNumber();
195 195
196 rha[i][j] = ibry*jbry*rh1; 196 rha[i][j] = ibry*jbry*rh1;
197 rha[j][i] = rha[i][j]; 197 rha[j][i] = rha[i][j];
198 198
199 // Coulomb terms 199 // Coulomb terms
200 200
201 G4double rrs2 = rr2[i][j] + epscl; 201 G4double rrs2 = rr2[i][j] + epscl;
202 G4double rrs = std::sqrt ( rrs2 ); 202 G4double rrs = std::sqrt ( rrs2 );
203 203
204 G4int icharge = system->GetParticipan 204 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
205 G4int jcharge = system->GetParticipan 205 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
206 206
207 G4double xerf = 0.0; 207 G4double xerf = 0.0;
208 // T. K. add this protection. 5.8 is 208 // T. K. add this protection. 5.8 is good enough for double
209 if ( rrs*c0sw < 5.8 ) 209 if ( rrs*c0sw < 5.8 )
210 { 210 {
211 #if defined WIN32-VC 211 #if defined WIN32-VC
212 xerf = CLHEP::HepStat::erf ( rrs*c 212 xerf = CLHEP::HepStat::erf ( rrs*c0sw );
213 #else 213 #else
214 xerf = std::erf ( rrs*c0sw ); 214 xerf = std::erf ( rrs*c0sw );
215 #endif 215 #endif
216 } 216 }
217 else 217 else
218 { 218 {
219 xerf = 1.0; 219 xerf = 1.0;
220 } 220 }
221 221
222 G4double erfij = xerf/rrs; 222 G4double erfij = xerf/rrs;
223 223
224 rhe[i][j] = icharge*jcharge * erfij; 224 rhe[i][j] = icharge*jcharge * erfij;
225 rhe[j][i] = rhe[i][j]; 225 rhe[j][i] = rhe[i][j];
226 rhc[i][j] = icharge*jcharge * ( - erf 226 rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2;
227 rhc[j][i] = rhc[i][j]; 227 rhc[j][i] = rhc[i][j];
228 } // i 228 } // i
229 } // j 229 } // j
230 } 230 }
231 231
232 void G4LightIonQMDMeanField::Cal2BodyQuantitie 232 void G4LightIonQMDMeanField::Cal2BodyQuantities( G4int i )
233 { 233 {
234 G4ThreeVector ri = system->GetParticipant( 234 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
235 G4LorentzVector p4i = system->GetParticipan 235 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum();
236 236
237 for ( G4int j = 0 ; j < system->GetTotalNum 237 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j )
238 { 238 {
239 if ( j == i ) { continue; } 239 if ( j == i ) { continue; }
240 240
241 G4ThreeVector rj = system->GetParticipan 241 G4ThreeVector rj = system->GetParticipant( j )->GetPosition();
242 G4LorentzVector p4j = system->GetPartici 242 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
243 243
244 G4ThreeVector rij = ri - rj; 244 G4ThreeVector rij = ri - rj;
245 G4ThreeVector pij = (p4i - p4j).v(); 245 G4ThreeVector pij = (p4i - p4j).v();
246 G4LorentzVector p4ij = p4i - p4j; 246 G4LorentzVector p4ij = p4i - p4j;
247 G4ThreeVector bij = ( p4i + p4j ).boostV 247 G4ThreeVector bij = ( p4i + p4j ).boostVector();
248 G4double gammaij = ( p4i + p4j ).gamma() 248 G4double gammaij = ( p4i + p4j ).gamma();
249 249
250 G4double eij = ( p4i + p4j ).e(); 250 G4double eij = ( p4i + p4j ).e();
251 251
252 G4double rbrb = rij*bij; 252 G4double rbrb = rij*bij;
253 G4double rij2 = rij*rij; 253 G4double rij2 = rij*rij;
254 G4double pij2 = pij*pij; 254 G4double pij2 = pij*pij;
255 255
256 rbrb = irelcr * rbrb; 256 rbrb = irelcr * rbrb;
257 G4double gamma2_ij = gammaij*gammaij; 257 G4double gamma2_ij = gammaij*gammaij;
258 258
259 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb 259 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb;
260 rr2[j][i] = rr2[i][j]; 260 rr2[j][i] = rr2[i][j];
261 261
262 rbij[i][j] = gamma2_ij * rbrb; 262 rbij[i][j] = gamma2_ij * rbrb;
263 rbij[j][i] = - rbij[i][j]; 263 rbij[j][i] = - rbij[i][j];
264 264
265 pp2[i][j] = pij2 265 pp2[i][j] = pij2
266 + irelcr * ( - G4Pow::GetInsta 266 + irelcr * ( - G4Pow::GetInstance()->powN ( p4i.e() - p4j.e() , 2 )
267 + gamma2_ij * G4Pow::GetInstan 267 + gamma2_ij * G4Pow::GetInstance()->powN ( ( ( p4i.m2() - p4j.m2() ) / eij ) , 2 ) );
268 268
269 pp2[j][i] = pp2[i][j]; 269 pp2[j][i] = pp2[i][j];
270 270
271 // Gauss term 271 // Gauss term
272 272
273 G4double expa1 = - rr2[i][j] * c0w; 273 G4double expa1 = - rr2[i][j] * c0w;
274 274
275 G4double rh1; 275 G4double rh1;
276 if ( expa1 > epsx ) 276 if ( expa1 > epsx )
277 { 277 {
278 rh1 = G4Exp( expa1 ); 278 rh1 = G4Exp( expa1 );
279 } 279 }
280 else 280 else
281 { 281 {
282 rh1 = 0.0; 282 rh1 = 0.0;
283 } 283 }
284 284
285 G4int ibry = system->GetParticipant(i)-> 285 G4int ibry = system->GetParticipant(i)->GetBaryonNumber();
286 G4int jbry = system->GetParticipant(j)-> 286 G4int jbry = system->GetParticipant(j)->GetBaryonNumber();
287 287
288 rha[i][j] = ibry*jbry*rh1; 288 rha[i][j] = ibry*jbry*rh1;
289 rha[j][i] = rha[i][j]; 289 rha[j][i] = rha[i][j];
290 290
291 // Coulomb terms 291 // Coulomb terms
292 292
293 G4double rrs2 = rr2[i][j] + epscl; 293 G4double rrs2 = rr2[i][j] + epscl;
294 G4double rrs = std::sqrt ( rrs2 ); 294 G4double rrs = std::sqrt ( rrs2 );
295 295
296 G4int icharge = system->GetParticipant(i 296 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
297 G4int jcharge = system->GetParticipant(j 297 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
298 298
299 G4double xerf = 0.0; 299 G4double xerf = 0.0;
300 // T. K. add this protection. 5.8 is goo 300 // T. K. add this protection. 5.8 is good enough for double
301 if ( rrs*c0sw < 5.8 ) 301 if ( rrs*c0sw < 5.8 )
302 { 302 {
303 #if defined WIN32-VC 303 #if defined WIN32-VC
304 xerf = CLHEP::HepStat::erf ( rrs*c0sw 304 xerf = CLHEP::HepStat::erf ( rrs*c0sw );
305 #else 305 #else
306 xerf = std::erf ( rrs*c0sw ); 306 xerf = std::erf ( rrs*c0sw );
307 #endif 307 #endif
308 } 308 }
309 else 309 else
310 { 310 {
311 xerf = 1.0; 311 xerf = 1.0;
312 } 312 }
313 313
314 G4double erfij = xerf/rrs; 314 G4double erfij = xerf/rrs;
315 315
316 rhe[i][j] = icharge*jcharge * erfij; 316 rhe[i][j] = icharge*jcharge * erfij;
317 rhe[j][i] = rhe[i][j]; 317 rhe[j][i] = rhe[i][j];
318 rhc[i][j] = icharge*jcharge * ( - erfij 318 rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2;
319 rhc[j][i] = rhc[i][j]; 319 rhc[j][i] = rhc[i][j];
320 } 320 }
321 } 321 }
322 322
323 void G4LightIonQMDMeanField::CalGraduate() 323 void G4LightIonQMDMeanField::CalGraduate()
324 { 324 {
325 ffr.resize( system->GetTotalNumberOfPartici 325 ffr.resize( system->GetTotalNumberOfParticipant() );
326 ffp.resize( system->GetTotalNumberOfPartici 326 ffp.resize( system->GetTotalNumberOfParticipant() );
327 rh3d.resize( system->GetTotalNumberOfPartic 327 rh3d.resize( system->GetTotalNumberOfParticipant() );
328 rh3d_tau.resize( system->GetTotalNumberOfPa 328 rh3d_tau.resize( system->GetTotalNumberOfParticipant() ); // Skyrme-QMD
329 329
330 for ( G4int i = 0 ; i < system->GetTotalNum 330 for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; ++i )
331 { 331 {
332 G4double rho3 = 0.0; 332 G4double rho3 = 0.0;
333 for ( G4int j = 0 ; j < system->GetTotal 333 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j )
334 { 334 {
335 rho3 += rha[j][i]; 335 rho3 += rha[j][i];
336 } 336 }
337 rh3d[i] = G4Pow::GetInstance()->powA ( r 337 rh3d[i] = G4Pow::GetInstance()->powA ( rho3 , pag );
338 rh3d_tau[i] = G4Pow::GetInstance()->powA 338 rh3d_tau[i] = G4Pow::GetInstance()->powA ( rho3 , pag_tau ); // Skyrme-QMD
339 } 339 }
340 340
341 for ( G4int i = 0 ; i < system->GetTotalNum 341 for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; ++i )
342 { 342 {
343 G4ThreeVector ri = system->GetParticipan 343 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
344 G4LorentzVector p4i = system->GetPartici 344 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum();
345 345
346 G4ThreeVector betai = p4i.v()/p4i.e(); 346 G4ThreeVector betai = p4i.v()/p4i.e();
347 347
348 // R-JQMD 348 // R-JQMD
349 G4double Vi = GetPotential( i ); 349 G4double Vi = GetPotential( i );
350 G4double p_zero = std::sqrt( p4i.e()*p4i 350 G4double p_zero = std::sqrt( p4i.e()*p4i.e() + 2*p4i.m()*Vi);
351 G4ThreeVector betai_R = p4i.v()/p_zero; 351 G4ThreeVector betai_R = p4i.v()/p_zero;
352 G4double mi_R = p4i.m()/p_zero; 352 G4double mi_R = p4i.m()/p_zero;
353 353
354 ffr[i] = betai_R; 354 ffr[i] = betai_R;
355 ffp[i] = G4ThreeVector( 0.0 ); 355 ffp[i] = G4ThreeVector( 0.0 );
356 356
357 for ( G4int j = 0 ; j < system->GetTotal 357 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j )
358 { 358 {
359 G4ThreeVector rj = system->GetPartici 359 G4ThreeVector rj = system->GetParticipant( j )->GetPosition();
360 G4LorentzVector p4j = system->GetPart 360 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
361 361
362 G4double eij = p4i.e() + p4j.e(); 362 G4double eij = p4i.e() + p4j.e();
363 363
364 G4int icharge = system->GetParticipan 364 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
365 G4int jcharge = system->GetParticipan 365 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
366 366
367 G4int inuc = system->GetParticipant(i 367 G4int inuc = system->GetParticipant(i)->GetNuc();
368 G4int jnuc = system->GetParticipant(j 368 G4int jnuc = system->GetParticipant(j)->GetNuc();
369 369
370 G4double fsij = 3.0/(2*wl) - rr2[j][i 370 G4double fsij = 3.0/(2*wl) - rr2[j][i]/(2*wl)/(2*wl); // Add for Skyrme-QMD
371 371
372 G4double ccpp = c0g * rha[j][i] 372 G4double ccpp = c0g * rha[j][i]
373 + c3g * rha[j][i] * ( r 373 + c3g * rha[j][i] * ( rh3d[j] + rh3d[i] )
374 + cg0 * rha[j][i]/wl 374 + cg0 * rha[j][i]/wl
375 + cg0 * rha[j][i] * fsi 375 + cg0 * rha[j][i] * fsij
376 + cgtau0 * rha[j][i] * 376 + cgtau0 * rha[j][i] * ( rh3d_tau[j] + rh3d_tau[i] )
377 + csg * rha[j][i] * jnu 377 + csg * rha[j][i] * jnuc * inuc
378 * ( 1. - 2. * std 378 * ( 1. - 2. * std::abs( jcharge - icharge ) )
379 * (1. - kappas * 379 * (1. - kappas * fsij + kappas / wl)
380 + cl * rhc[j][i]; 380 + cl * rhc[j][i];
381 381
382 ccpp *= mi_R; 382 ccpp *= mi_R;
383 383
384 G4double grbb = - rbij[j][i]; 384 G4double grbb = - rbij[j][i];
385 G4double ccrr = grbb * ccpp / eij; 385 G4double ccrr = grbb * ccpp / eij;
386 386
387 G4ThreeVector rij = ri - rj; 387 G4ThreeVector rij = ri - rj;
388 G4ThreeVector betaij = ( p4i + p4j ) 388 G4ThreeVector betaij = ( p4i + p4j ).v()/eij;
389 G4ThreeVector cij = betaij - betai; 389 G4ThreeVector cij = betaij - betai;
390 390
391 ffr[i] = ffr[i] + 2*ccrr* ( rij + grb 391 ffr[i] = ffr[i] + 2*ccrr* ( rij + grbb*cij );
392 ffp[i] = ffp[i] - 2*ccpp* ( rij + grb 392 ffp[i] = ffp[i] - 2*ccpp* ( rij + grbb*betaij );
393 } 393 }
394 } 394 }
395 } 395 }
396 396
397 G4double G4LightIonQMDMeanField::GetPotential( 397 G4double G4LightIonQMDMeanField::GetPotential( G4int i )
398 { 398 {
399 G4int n = system->GetTotalNumberOfParticipa 399 G4int n = system->GetTotalNumberOfParticipant();
400 400
401 G4double rhoa = 0.0; 401 G4double rhoa = 0.0;
402 G4double rho3 = 0.0; 402 G4double rho3 = 0.0;
403 G4double fsij_rhoa = 0.0; // Skyrme-QMD 403 G4double fsij_rhoa = 0.0; // Skyrme-QMD
404 //G4double fsij_rhos = 0.0; // Skyrme-QMD 404 //G4double fsij_rhos = 0.0; // Skyrme-QMD
405 G4double rho3_tau = 0.0; // Skyrme-QMD 405 G4double rho3_tau = 0.0; // Skyrme-QMD
406 G4double rhos = 0.0; 406 G4double rhos = 0.0;
407 G4double rhoc = 0.0; 407 G4double rhoc = 0.0;
408 408
409 G4int icharge = system->GetParticipant(i)-> 409 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
410 G4int inuc = system->GetParticipant(i)->Get 410 G4int inuc = system->GetParticipant(i)->GetNuc();
411 411
412 for ( G4int j = 0 ; j < n ; ++j ) 412 for ( G4int j = 0 ; j < n ; ++j )
413 { 413 {
414 G4int jcharge = system->GetParticipant(j 414 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
415 G4int jnuc = system->GetParticipant(j)-> 415 G4int jnuc = system->GetParticipant(j)->GetNuc();
416 G4double fsij = 3.0/(2*wl) - rr2[j][i]/( 416 G4double fsij = 3.0/(2*wl) - rr2[j][i]/(2*wl)/(2*wl); // Add for Skyrme-QMD
417 417
418 rhoa += rha[j][i]; 418 rhoa += rha[j][i];
419 fsij_rhoa += fsij * rha[j][i]; // Skyrme 419 fsij_rhoa += fsij * rha[j][i]; // Skyrme-QMD
420 rhoc += rhe[j][i]; 420 rhoc += rhe[j][i];
421 rhos += rha[j][i] * jnuc * inuc 421 rhos += rha[j][i] * jnuc * inuc
422 * ( 1. - 2. * std::abs 422 * ( 1. - 2. * std::abs( jcharge - icharge ) ) // Skyrme-QMD
423 * (1. - kappas * fsij) 423 * (1. - kappas * fsij); // Skyrme-QMD
424 } 424 }
425 425
426 rho3 = G4Pow::GetInstance()->powA ( rhoa , 426 rho3 = G4Pow::GetInstance()->powA ( rhoa , gamm );
427 rho3_tau = G4Pow::GetInstance()->powA ( rho 427 rho3_tau = G4Pow::GetInstance()->powA ( rhoa , eta );
428 428
429 G4double potential = c0 * rhoa 429 G4double potential = c0 * rhoa
430 + c3 * rho3 430 + c3 * rho3
431 + g0 * fsij_rhoa // Sky 431 + g0 * fsij_rhoa // Skyrme-QMD
432 //+ g0iso * fsij_rhos / 432 //+ g0iso * fsij_rhos // Skyrme-QMD
433 + gtau0 * rho3_tau // S 433 + gtau0 * rho3_tau // Skyrme-QMD
434 + cs * rhos 434 + cs * rhos
435 + cl * rhoc; 435 + cl * rhoc;
436 return potential; 436 return potential;
437 } 437 }
438 438
439 G4double G4LightIonQMDMeanField::GetTotalPoten 439 G4double G4LightIonQMDMeanField::GetTotalPotential()
440 { 440 {
441 G4int n = system->GetTotalNumberOfParticipa 441 G4int n = system->GetTotalNumberOfParticipant();
442 442
443 std::vector < G4double > rhoa ( n , 0.0 ); 443 std::vector < G4double > rhoa ( n , 0.0 );
444 std::vector < G4double > rho3 ( n , 0.0 ); 444 std::vector < G4double > rho3 ( n , 0.0 );
445 std::vector < G4double > rho3_tau ( n , 0.0 445 std::vector < G4double > rho3_tau ( n , 0.0 ); // Skyrme-QMD
446 //std::vector < G4double > fsij_rhos ( n , 446 //std::vector < G4double > fsij_rhos ( n , 0.0 ); // Skyrme-QMD
447 std::vector < G4double > fsij_rhoa ( n , 0. 447 std::vector < G4double > fsij_rhoa ( n , 0.0 ); // Skyrme-QMD
448 std::vector < G4double > rhos ( n , 0.0 ); 448 std::vector < G4double > rhos ( n , 0.0 );
449 std::vector < G4double > rhoc ( n , 0.0 ); 449 std::vector < G4double > rhoc ( n , 0.0 );
450 450
451 for ( G4int i = 0 ; i < n ; ++i ) 451 for ( G4int i = 0 ; i < n ; ++i )
452 { 452 {
453 G4int icharge = system->GetParticipant(i 453 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
454 G4int inuc = system->GetParticipant(i)-> 454 G4int inuc = system->GetParticipant(i)->GetNuc();
455 455
456 for ( G4int j = 0 ; j < n ; ++j ) 456 for ( G4int j = 0 ; j < n ; ++j )
457 { 457 {
458 G4int jcharge = system->GetParticipan 458 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
459 G4int jnuc = system->GetParticipant(j 459 G4int jnuc = system->GetParticipant(j)->GetNuc();
460 G4double fsij = 3.0/(2*wl) - rr2[j][i 460 G4double fsij = 3.0/(2*wl) - rr2[j][i]/(2*wl)/(2*wl); // Add for Skyrme-QMD
461 461
462 rhoa[i] += rha[j][i]; 462 rhoa[i] += rha[j][i];
463 fsij_rhoa[i] += fsij * rha[j][i]; // 463 fsij_rhoa[i] += fsij * rha[j][i]; // Skyrme-QMD
464 rhoc[i] += rhe[j][i]; 464 rhoc[i] += rhe[j][i];
465 rhos[i] += rha[j][i] * jnuc * inuc 465 rhos[i] += rha[j][i] * jnuc * inuc
466 //* ( 1 - 2 * st 466 //* ( 1 - 2 * std::abs ( jcharge - icharge ) );
467 * ( 1. - 2. * st 467 * ( 1. - 2. * std::abs( jcharge - icharge ) ) // Skyrme-QMD
468 * (1. - kappas * 468 * (1. - kappas * fsij); // Skyrme-QMD
469 //fsij_rhos[i] += fsij * rha[j][i] * 469 //fsij_rhos[i] += fsij * rha[j][i] * jnuc * inuc
470 //* ( 1. - 2. 470 //* ( 1. - 2. * std::abs( jcharge - icharge ) ) // Skyrme-QMD
471 //* (1. - kapp 471 //* (1. - kappas * fsij); // Skyrme-QMD
472 } 472 }
473 473
474 rho3[i] = G4Pow::GetInstance()->powA ( r 474 rho3[i] = G4Pow::GetInstance()->powA ( rhoa[i] , gamm );
475 rho3_tau[i] = G4Pow::GetInstance()->powA 475 rho3_tau[i] = G4Pow::GetInstance()->powA ( rhoa[i] , eta );
476 } 476 }
477 477
478 G4double potential = c0 * std::accumulate( 478 G4double potential = c0 * std::accumulate( rhoa.cbegin() , rhoa.cend() , 0.0 )
479 + c3 * std::accumulate( 479 + c3 * std::accumulate( rho3.cbegin() , rho3.cend() , 0.0 )
480 + g0 * std::accumulate( 480 + g0 * std::accumulate( fsij_rhoa.cbegin() , fsij_rhoa.cend() , 0.0 )
481 //+ g0iso * std::accumu 481 //+ g0iso * std::accumulate( fsij_rhos.cbegin() , fsij_rhos.cend() , 0.0 )
482 + gtau0 * std::accumula 482 + gtau0 * std::accumulate( rho3_tau.cbegin() , rho3_tau.cend() , 0.0 )
483 + cs * std::accumulate( 483 + cs * std::accumulate( rhos.cbegin() , rhos.cend() , 0.0 )
484 + cl * std::accumulate( 484 + cl * std::accumulate( rhoc.cbegin() , rhoc.cend() , 0.0 );
485 485
486 return potential; 486 return potential;
487 } 487 }
488 488
489 G4double G4LightIonQMDMeanField::GetSingleEner 489 G4double G4LightIonQMDMeanField::GetSingleEnergy( G4int j )
490 { 490 {
491 G4LorentzVector p4j = system->GetParticipa 491 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
492 G4double emass = p4j.m(); 492 G4double emass = p4j.m();
493 G4double ekinal2 = p4j.e()*p4j.e(); 493 G4double ekinal2 = p4j.e()*p4j.e();
494 G4double esingle = std::sqrt(ekinal2 + 2*e 494 G4double esingle = std::sqrt(ekinal2 + 2*emass*GetPotential(j));
495 return esingle; 495 return esingle;
496 } 496 }
497 497
498 G4double G4LightIonQMDMeanField::GetTotalEnerg 498 G4double G4LightIonQMDMeanField::GetTotalEnergy()
499 { 499 {
500 500
501 G4int n = system->GetTotalNumberOfParticip 501 G4int n = system->GetTotalNumberOfParticipant();
502 G4double etotal = 0.0; 502 G4double etotal = 0.0;
503 for ( int j = 0 ; j < n ; j++ ) 503 for ( int j = 0 ; j < n ; j++ )
504 { 504 {
505 G4LorentzVector p4j = system->GetParti 505 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
506 G4double emass = p4j.m(); 506 G4double emass = p4j.m();
507 G4double ekinal2 = p4j.e()*p4j.e(); 507 G4double ekinal2 = p4j.e()*p4j.e();
508 etotal += std::sqrt(ekinal2 + 2*emass* 508 etotal += std::sqrt(ekinal2 + 2*emass*GetPotential(j));
509 } 509 }
510 return etotal; 510 return etotal;
511 511
512 } 512 }
513 513
514 G4double G4LightIonQMDMeanField::calPauliBlock 514 G4double G4LightIonQMDMeanField::calPauliBlockingFactor( G4int i )
515 { 515 {
516 // i is supposed beyond total number of Par 516 // i is supposed beyond total number of Participant()
517 517
518 G4double pf = 0.0; 518 G4double pf = 0.0;
519 G4int icharge = system->GetParticipant(i)-> 519 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
520 520
521 for ( G4int j = 0 ; j < system->GetTotalNum 521 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j )
522 { 522 {
523 G4int jcharge = system->GetParticipant(j 523 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
524 G4int jnuc = system->GetParticipant(j)-> 524 G4int jnuc = system->GetParticipant(j)->GetNuc();
525 525
526 if ( jcharge == icharge && jnuc == 1 ) 526 if ( jcharge == icharge && jnuc == 1 )
527 { 527 {
528 G4double expa = -rr2[i][j]*cpw; 528 G4double expa = -rr2[i][j]*cpw;
529 if ( expa > epsx ) 529 if ( expa > epsx )
530 { 530 {
531 expa = expa - pp2[i][j]*cph; 531 expa = expa - pp2[i][j]*cph;
532 if ( expa > epsx ) 532 if ( expa > epsx )
533 { 533 {
534 pf = pf + G4Exp ( expa ); 534 pf = pf + G4Exp ( expa );
535 } 535 }
536 } 536 }
537 } 537 }
538 } 538 }
539 539
540 return ( pf - 1.0 ) * cpc; 540 return ( pf - 1.0 ) * cpc;
541 } 541 }
542 542
543 G4bool G4LightIonQMDMeanField::IsPauliBlocked( 543 G4bool G4LightIonQMDMeanField::IsPauliBlocked( G4int i )
544 { 544 {
545 G4bool result = false; 545 G4bool result = false;
546 546
547 if ( system->GetParticipant( i )->GetNuc() 547 if ( system->GetParticipant( i )->GetNuc() == 1 )
548 { 548 {
549 G4double pf = calPauliBlockingFactor( i 549 G4double pf = calPauliBlockingFactor( i );
550 G4double rand = G4UniformRand(); 550 G4double rand = G4UniformRand();
551 if ( pf > rand ) { result = true; } 551 if ( pf > rand ) { result = true; }
552 } 552 }
553 553
554 return result; 554 return result;
555 } 555 }
556 556
557 void G4LightIonQMDMeanField::DoPropagation( G4 557 void G4LightIonQMDMeanField::DoPropagation( G4double dt )
558 { 558 {
559 G4double cc2 = 1.0; 559 G4double cc2 = 1.0;
560 G4double cc1 = 1.0 - cc2; 560 G4double cc1 = 1.0 - cc2;
561 G4double cc3 = 1.0 / 2.0 / cc2; 561 G4double cc3 = 1.0 / 2.0 / cc2;
562 562
563 G4double dt3 = dt * cc3; 563 G4double dt3 = dt * cc3;
564 G4double dt1 = dt * ( cc1 - cc3 ); 564 G4double dt1 = dt * ( cc1 - cc3 );
565 G4double dt2 = dt * cc2; 565 G4double dt2 = dt * cc2;
566 566
567 CalGraduate(); 567 CalGraduate();
568 568
569 G4int n = system->GetTotalNumberOfParticipa 569 G4int n = system->GetTotalNumberOfParticipant();
570 570
571 // 1st Step 571 // 1st Step
572 572
573 std::vector< G4ThreeVector > f0r, f0p; 573 std::vector< G4ThreeVector > f0r, f0p;
574 f0r.resize( n ); 574 f0r.resize( n );
575 f0p.resize( n ); 575 f0p.resize( n );
576 576
577 for ( G4int i = 0 ; i < n ; ++i ) 577 for ( G4int i = 0 ; i < n ; ++i )
578 { 578 {
579 G4ThreeVector ri = system->GetParticipan 579 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
580 G4ThreeVector p3i = system->GetParticipa 580 G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum();
581 581
582 ri += dt3* ffr[i]; 582 ri += dt3* ffr[i];
583 p3i += dt3* ffp[i]; 583 p3i += dt3* ffp[i];
584 584
585 f0r[i] = ffr[i]; 585 f0r[i] = ffr[i];
586 f0p[i] = ffp[i]; 586 f0p[i] = ffp[i];
587 587
588 system->GetParticipant( i )->SetPosition 588 system->GetParticipant( i )->SetPosition( ri );
589 system->GetParticipant( i )->SetMomentum 589 system->GetParticipant( i )->SetMomentum( p3i );
590 590
591 // we do not need set total momentum by 591 // we do not need set total momentum by ourselvs
592 } 592 }
593 593
594 // 2nd Step 594 // 2nd Step
595 595
596 Cal2BodyQuantities(); 596 Cal2BodyQuantities();
597 CalGraduate(); 597 CalGraduate();
598 598
599 for ( G4int i = 0 ; i < n ; ++i ) 599 for ( G4int i = 0 ; i < n ; ++i )
600 { 600 {
601 G4ThreeVector ri = system->GetParticipan 601 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
602 G4ThreeVector p3i = system->GetParticipa 602 G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum();
603 603
604 ri += dt1* f0r[i] + dt2* ffr[i]; 604 ri += dt1* f0r[i] + dt2* ffr[i];
605 p3i += dt1* f0p[i] + dt2* ffp[i]; 605 p3i += dt1* f0p[i] + dt2* ffp[i];
606 606
607 system->GetParticipant( i )->SetPosition 607 system->GetParticipant( i )->SetPosition( ri );
608 system->GetParticipant( i )->SetMomentum 608 system->GetParticipant( i )->SetMomentum( p3i );
609 609
610 // we do not need set total momentum by 610 // we do not need set total momentum by ourselvs
611 } 611 }
612 612
613 Cal2BodyQuantities(); 613 Cal2BodyQuantities();
614 } 614 }
615 615
616 std::vector< G4LightIonQMDNucleus* > G4LightIo 616 std::vector< G4LightIonQMDNucleus* > G4LightIonQMDMeanField::DoClusterJudgment()
617 { 617 {
618 Cal2BodyQuantities(); 618 Cal2BodyQuantities();
619 619
620 G4double cpf2 = G4Pow::GetInstance()->A23 ( 620 G4double cpf2 = G4Pow::GetInstance()->A23 ( 1.5 * pi*pi * G4Pow::GetInstance()->powA ( 4.0 * pi * wl , -1.5 ) ) * hbc * hbc;
621 G4double rcc2 = rclds*rclds; 621 G4double rcc2 = rclds*rclds;
622 622
623 G4int n = system->GetTotalNumberOfParticipa 623 G4int n = system->GetTotalNumberOfParticipant();
624 std::vector < G4double > rhoa; 624 std::vector < G4double > rhoa;
625 rhoa.resize ( n ); 625 rhoa.resize ( n );
626 626
627 for ( G4int i = 0 ; i < n ; ++i ) 627 for ( G4int i = 0 ; i < n ; ++i )
628 { 628 {
629 rhoa[i] = 0.0; 629 rhoa[i] = 0.0;
630 630
631 if ( system->GetParticipant( i )->GetBary 631 if ( system->GetParticipant( i )->GetBaryonNumber() == 1 )
632 { 632 {
633 for ( G4int j = 0 ; j < n ; ++j ) 633 for ( G4int j = 0 ; j < n ; ++j )
634 { 634 {
635 if ( system->GetParticipant( j )->Get 635 if ( system->GetParticipant( j )->GetBaryonNumber() == 1 )
636 rhoa[i] += rha[i][j]; 636 rhoa[i] += rha[i][j];
637 } 637 }
638 } 638 }
639 639
640 rhoa[i] = G4Pow::GetInstance()->A13 ( rho 640 rhoa[i] = G4Pow::GetInstance()->A13 ( rhoa[i] + 1 );
641 } 641 }
642 642
643 // identification of the cluster 643 // identification of the cluster
644 std::vector < G4bool > is_already_belong_so 644 std::vector < G4bool > is_already_belong_some_cluster;
645 645
646 // cluster_id participant_id 646 // cluster_id participant_id
647 std::multimap < G4int , G4int > comb_map; 647 std::multimap < G4int , G4int > comb_map;
648 std::multimap < G4int , G4int > assign_map; 648 std::multimap < G4int , G4int > assign_map;
649 assign_map.clear(); 649 assign_map.clear();
650 650
651 std::vector < G4int > mascl; 651 std::vector < G4int > mascl;
652 std::vector < G4int > num; 652 std::vector < G4int > num;
653 mascl.resize ( n ); 653 mascl.resize ( n );
654 num.resize ( n ); 654 num.resize ( n );
655 is_already_belong_some_cluster.resize ( n ) 655 is_already_belong_some_cluster.resize ( n );
656 656
657 std::vector < G4int > is_assigned_to ( n , 657 std::vector < G4int > is_assigned_to ( n , -1 );
658 std::multimap < G4int , G4int > clusters; 658 std::multimap < G4int , G4int > clusters;
659 659
660 for ( G4int i = 0 ; i < n ; ++i ) 660 for ( G4int i = 0 ; i < n ; ++i )
661 { 661 {
662 mascl[i] = 1; 662 mascl[i] = 1;
663 num[i] = 1; 663 num[i] = 1;
664 is_already_belong_some_cluster[i] = false 664 is_already_belong_some_cluster[i] = false;
665 } 665 }
666 666
667 G4int ichek = 1; 667 G4int ichek = 1;
668 G4int id = 0; 668 G4int id = 0;
669 G4int cluster_id = -1; 669 G4int cluster_id = -1;
670 for ( G4int i = 0 ; i < n-1 ; ++i ) 670 for ( G4int i = 0 ; i < n-1 ; ++i )
671 { 671 {
672 G4bool hasThisCompany = false; 672 G4bool hasThisCompany = false;
673 673
674 if ( system->GetParticipant( i )->GetBary 674 if ( system->GetParticipant( i )->GetBaryonNumber() == 1 )
675 { 675 {
676 G4int j1 = i + 1; 676 G4int j1 = i + 1;
677 for ( G4int j = j1 ; j < n ; ++j ) 677 for ( G4int j = j1 ; j < n ; ++j )
678 { 678 {
679 std::vector < G4int > cluster_partici 679 std::vector < G4int > cluster_participants;
680 if ( system->GetParticipant( j )->Get 680 if ( system->GetParticipant( j )->GetBaryonNumber() == 1 )
681 { 681 {
682 G4double rdist2 = rr2[ i ][ j ]; 682 G4double rdist2 = rr2[ i ][ j ];
683 G4double pdist2 = pp2[ i ][ j ]; 683 G4double pdist2 = pp2[ i ][ j ];
684 G4double pcc2 = cpf2 684 G4double pcc2 = cpf2
685 * ( rhoa[ i ] + rhoa[ 685 * ( rhoa[ i ] + rhoa[ j ] )
686 * ( rhoa[ i ] + rhoa[ 686 * ( rhoa[ i ] + rhoa[ j ] );
687 687
688 // Check phase space: close enough? 688 // Check phase space: close enough?
689 if ( rdist2 < rcc2 && pdist2 < pcc2 689 if ( rdist2 < rcc2 && pdist2 < pcc2 )
690 { 690 {
691 if ( is_assigned_to [ j ] == -1 ) 691 if ( is_assigned_to [ j ] == -1 )
692 { 692 {
693 if ( is_assigned_to [ i ] == -1 693 if ( is_assigned_to [ i ] == -1 )
694 { 694 {
695 if ( clusters.size() != 0 ) 695 if ( clusters.size() != 0 )
696 { 696 {
697 id = clusters.rbegin()->fir 697 id = clusters.rbegin()->first + 1;
698 } 698 }
699 else 699 else
700 { 700 {
701 id = 0; 701 id = 0;
702 } 702 }
703 clusters.insert ( std::multim 703 clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , i ) );
704 is_assigned_to [ i ] = id; 704 is_assigned_to [ i ] = id;
705 clusters.insert ( std::multim 705 clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , j ) );
706 is_assigned_to [ j ] = id; 706 is_assigned_to [ j ] = id;
707 } 707 }
708 else 708 else
709 { 709 {
710 clusters.insert ( std::multim 710 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ i ] , j ) );
711 is_assigned_to [ j ] = is_ass 711 is_assigned_to [ j ] = is_assigned_to [ i ];
712 } 712 }
713 } 713 }
714 else 714 else
715 { 715 {
716 // j is already belong to some 716 // j is already belong to some cluster
717 if ( is_assigned_to [ i ] == -1 717 if ( is_assigned_to [ i ] == -1 )
718 { 718 {
719 clusters.insert ( std::multim 719 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , i ) );
720 is_assigned_to [ i ] = is_ass 720 is_assigned_to [ i ] = is_assigned_to [ j ];
721 } 721 }
722 else 722 else
723 { 723 {
724 // i has companion 724 // i has companion
725 if ( is_assigned_to [ i ] != 725 if ( is_assigned_to [ i ] != is_assigned_to [ j ] )
726 { 726 {
727 // move companions to the c 727 // move companions to the cluster
728 std::multimap< G4int , G4in 728 std::multimap< G4int , G4int > clusters_tmp;
729 G4int target_cluster_id; 729 G4int target_cluster_id;
730 if ( is_assigned_to [ i ] > 730 if ( is_assigned_to [ i ] > is_assigned_to [ j ] )
731 { 731 {
732 target_cluster_id = is_as 732 target_cluster_id = is_assigned_to [ i ];
733 } 733 }
734 else 734 else
735 { 735 {
736 target_cluster_id = is_as 736 target_cluster_id = is_assigned_to [ j ];
737 } 737 }
738 for ( auto it = clusters.cb 738 for ( auto it = clusters.cbegin() ; it != clusters.cend() ; ++it )
739 { 739 {
740 if ( it->first == target_ 740 if ( it->first == target_cluster_id )
741 { 741 {
742 is_assigned_to [ it->se 742 is_assigned_to [ it->second ] = is_assigned_to [ j ];
743 clusters_tmp.insert ( s 743 clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , it->second ) );
744 } 744 }
745 else 745 else
746 { 746 {
747 clusters_tmp.insert ( s 747 clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( it->first , it->second ) );
748 } 748 }
749 } 749 }
750 clusters = std::move(cluste << 750 clusters = clusters_tmp;
751 } 751 }
752 } 752 }
753 } 753 }
754 754
755 comb_map.insert( std::multimap<G4 755 comb_map.insert( std::multimap<G4int,G4int>::value_type ( i , j ) );
756 cluster_participants.push_back ( 756 cluster_participants.push_back ( j );
757 757
758 if ( assign_map.find( cluster_id 758 if ( assign_map.find( cluster_id ) == assign_map.end() )
759 { 759 {
760 is_already_belong_some_cluster[ 760 is_already_belong_some_cluster[i] = true;
761 assign_map.insert ( std::multim 761 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , i ) );
762 hasThisCompany = true; 762 hasThisCompany = true;
763 } 763 }
764 assign_map.insert ( std::multimap 764 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , j ) );
765 is_already_belong_some_cluster[j] 765 is_already_belong_some_cluster[j] = true;
766 } 766 }
767 767
768 if ( ichek == i ) 768 if ( ichek == i )
769 { 769 {
770 ++ichek; 770 ++ichek;
771 } 771 }
772 } 772 }
773 } 773 }
774 } 774 }
775 if ( hasThisCompany == true ) { ++cluster 775 if ( hasThisCompany == true ) { ++cluster_id; }
776 } 776 }
777 777
778 // sort 778 // sort
779 // Heavy cluster comes first 779 // Heavy cluster comes first
780 // size cluster_id 780 // size cluster_id
781 std::multimap< G4int , G4int > sorted_clust 781 std::multimap< G4int , G4int > sorted_cluster_map;
782 for ( G4int i = 0 ; i <= id ; ++i ) // << 782 for ( G4int i = 0 ; i <= id ; ++i ) // << "<=" because id is highest cluster nubmer.
783 { 783 {
784 sorted_cluster_map.insert ( std::multimap 784 sorted_cluster_map.insert ( std::multimap<G4int,G4int>::value_type ( (G4int) clusters.count( i ) , i ) );
785 } 785 }
786 786
787 // create nucleus from divided clusters 787 // create nucleus from divided clusters
788 std::vector < G4LightIonQMDNucleus* > resul 788 std::vector < G4LightIonQMDNucleus* > result;
789 for ( auto it = sorted_cluster_map.crbegin( 789 for ( auto it = sorted_cluster_map.crbegin(); it != sorted_cluster_map.crend(); ++it )
790 { 790 {
791 if ( it->first != 0 ) 791 if ( it->first != 0 )
792 { 792 {
793 G4LightIonQMDNucleus* nucleus = new G4L 793 G4LightIonQMDNucleus* nucleus = new G4LightIonQMDNucleus();
794 for ( auto itt = clusters.cbegin(); itt 794 for ( auto itt = clusters.cbegin(); itt != clusters.cend(); ++itt )
795 { 795 {
796 if ( it->second == itt->first ) 796 if ( it->second == itt->first )
797 { 797 {
798 nucleus->SetParticipant( system->Ge 798 nucleus->SetParticipant( system->GetParticipant ( itt->second ) );
799 } 799 }
800 } 800 }
801 result.push_back( nucleus ); 801 result.push_back( nucleus );
802 } 802 }
803 } 803 }
804 804
805 // delete participants from current system 805 // delete participants from current system
806 for ( auto it = result.cbegin(); it != resu 806 for ( auto it = result.cbegin(); it != result.cend(); ++it )
807 { 807 {
808 system->SubtractSystem ( *it ); 808 system->SubtractSystem ( *it );
809 } 809 }
810 810
811 return result; 811 return result;
812 } 812 }
813 813
814 void G4LightIonQMDMeanField::Update() 814 void G4LightIonQMDMeanField::Update()
815 { 815 {
816 SetSystem( system ); 816 SetSystem( system );
817 } 817 }
818 818