| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4QMDMeanField implementation
27 //
28 // Author: Tatsumi Koi (SLAC), 29 March 2007
29 // --------------------------------------------------------------------
30
31 #include <map>
32 #include <algorithm>
33 #include <numeric>
34
35 #include <cmath>
36 #include <CLHEP/Random/Stat.h>
37
38 #include "G4QMDMeanField.hh"
39 #include "G4QMDParameters.hh"
40 #include "G4Exp.hh"
41 #include "G4Pow.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "Randomize.hh"
44
45 G4QMDMeanField::G4QMDMeanField()
46 {
47 G4QMDParameters* parameters = G4QMDParameters::GetInstance();
48 wl = parameters->Get_wl();
49 cl = parameters->Get_cl();
50 rho0 = parameters->Get_rho0();
51 hbc = parameters->Get_hbc();
52 gamm = parameters->Get_gamm();
53
54 cpw = parameters->Get_cpw();
55 cph = parameters->Get_cph();
56 cpc = parameters->Get_cpc();
57
58 c0 = parameters->Get_c0();
59 c3 = parameters->Get_c3();
60 cs = parameters->Get_cs();
61
62 // distance
63 c0w = 1.0/4.0/wl;
64 c0sw = std::sqrt( c0w );
65 clw = 2.0 / std::sqrt ( 4.0 * pi * wl );
66
67 // graduate
68 c0g = - c0 / ( 2.0 * wl );
69 c3g = - c3 / ( 4.0 * wl ) * gamm;
70 csg = - cs / ( 2.0 * wl );
71 pag = gamm - 1;
72
73 system = nullptr; // will be set through SetSystem() method
74 }
75
76 void G4QMDMeanField::SetSystem ( G4QMDSystem* aSystem )
77 {
78 system = aSystem;
79
80 G4int n = system->GetTotalNumberOfParticipant();
81
82 pp2.clear();
83 rr2.clear();
84 rbij.clear();
85 rha.clear();
86 rhe.clear();
87 rhc.clear();
88
89 rr2.resize( n );
90 pp2.resize( n );
91 rbij.resize( n );
92 rha.resize( n );
93 rhe.resize( n );
94 rhc.resize( n );
95
96 for ( G4int i = 0 ; i < n ; ++i )
97 {
98 rr2[i].resize( n );
99 pp2[i].resize( n );
100 rbij[i].resize( n );
101 rha[i].resize( n );
102 rhe[i].resize( n );
103 rhc[i].resize( n );
104 }
105
106 ffr.clear();
107 ffp.clear();
108 rh3d.clear();
109
110 ffr.resize( n );
111 ffp.resize( n );
112 rh3d.resize( n );
113
114 Cal2BodyQuantities();
115 }
116
117 void G4QMDMeanField::SetNucleus ( G4QMDNucleus* aNucleus )
118 {
119 SetSystem( aNucleus );
120
121 G4double totalPotential = GetTotalPotential();
122 aNucleus->SetTotalPotential( totalPotential );
123 aNucleus->CalEnergyAndAngularMomentumInCM();
124 }
125
126 void G4QMDMeanField::Cal2BodyQuantities()
127 {
128 if ( system->GetTotalNumberOfParticipant() < 2 ) { return; }
129
130 for ( G4int j = 1 ; j < system->GetTotalNumberOfParticipant() ; ++j )
131 {
132 G4ThreeVector rj = system->GetParticipant( j )->GetPosition();
133 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
134
135 for ( G4int i = 0 ; i < j ; ++i )
136 {
137 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
138 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum();
139
140 G4ThreeVector rij = ri - rj;
141 G4ThreeVector pij = (p4i - p4j).v();
142 G4LorentzVector p4ij = p4i - p4j;
143 G4ThreeVector bij = ( p4i + p4j ).boostVector();
144 G4double gammaij = ( p4i + p4j ).gamma();
145
146 G4double eij = ( p4i + p4j ).e();
147
148 G4double rbrb = rij*bij;
149 G4double rij2 = rij*rij;
150 G4double pij2 = pij*pij;
151
152 rbrb = irelcr * rbrb;
153 G4double gamma2_ij = gammaij*gammaij;
154
155 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb;
156 rr2[j][i] = rr2[i][j];
157
158 rbij[i][j] = gamma2_ij * rbrb;
159 rbij[j][i] = - rbij[i][j];
160
161 pp2[i][j] = pij2
162 + irelcr * ( - G4Pow::GetInstance()->powN ( p4i.e() - p4j.e() , 2 )
163 + gamma2_ij * G4Pow::GetInstance()->powN ( ( ( p4i.m2() - p4j.m2() )
164 / eij ) , 2 ) );
165 pp2[j][i] = pp2[i][j];
166
167 // Gauss term
168
169 G4double expa1 = - rr2[i][j] * c0w;
170
171 G4double rh1;
172 if ( expa1 > epsx )
173 {
174 rh1 = G4Exp( expa1 );
175 }
176 else
177 {
178 rh1 = 0.0;
179 }
180
181 G4int ibry = system->GetParticipant(i)->GetBaryonNumber();
182 G4int jbry = system->GetParticipant(j)->GetBaryonNumber();
183
184 rha[i][j] = ibry*jbry*rh1;
185 rha[j][i] = rha[i][j];
186
187 // Coulomb terms
188
189 G4double rrs2 = rr2[i][j] + epscl;
190 G4double rrs = std::sqrt ( rrs2 );
191
192 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
193 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
194
195 G4double xerf = 0.0;
196 // T. K. add this protection. 5.8 is good enough for double
197 if ( rrs*c0sw < 5.8 )
198 {
199 #if defined WIN32-VC
200 xerf = CLHEP::HepStat::erf ( rrs*c0sw );
201 #else
202 xerf = std::erf ( rrs*c0sw );
203 #endif
204 }
205 else
206 {
207 xerf = 1.0;
208 }
209
210 G4double erfij = xerf/rrs;
211
212 rhe[i][j] = icharge*jcharge * erfij;
213 rhe[j][i] = rhe[i][j];
214 rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2;
215 rhc[j][i] = rhc[i][j];
216 } // i
217 } // j
218 }
219
220 void G4QMDMeanField::Cal2BodyQuantities( G4int i )
221 {
222 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
223 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum();
224
225 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j )
226 {
227 if ( j == i ) { continue; }
228
229 G4ThreeVector rj = system->GetParticipant( j )->GetPosition();
230 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
231
232 G4ThreeVector rij = ri - rj;
233 G4ThreeVector pij = (p4i - p4j).v();
234 G4LorentzVector p4ij = p4i - p4j;
235 G4ThreeVector bij = ( p4i + p4j ).boostVector();
236 G4double gammaij = ( p4i + p4j ).gamma();
237
238 G4double eij = ( p4i + p4j ).e();
239
240 G4double rbrb = rij*bij;
241 G4double rij2 = rij*rij;
242 G4double pij2 = pij*pij;
243
244 rbrb = irelcr * rbrb;
245 G4double gamma2_ij = gammaij*gammaij;
246
247 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb;
248 rr2[j][i] = rr2[i][j];
249
250 rbij[i][j] = gamma2_ij * rbrb;
251 rbij[j][i] = - rbij[i][j];
252
253 pp2[i][j] = pij2
254 + irelcr * ( - G4Pow::GetInstance()->powN ( p4i.e() - p4j.e() , 2 )
255 + gamma2_ij * G4Pow::GetInstance()->powN ( ( ( p4i.m2() - p4j.m2() )
256 / eij ) , 2 ) );
257 pp2[j][i] = pp2[i][j];
258
259 // Gauss term
260
261 G4double expa1 = - rr2[i][j] * c0w;
262
263 G4double rh1;
264 if ( expa1 > epsx )
265 {
266 rh1 = G4Exp( expa1 );
267 }
268 else
269 {
270 rh1 = 0.0;
271 }
272
273 G4int ibry = system->GetParticipant(i)->GetBaryonNumber();
274 G4int jbry = system->GetParticipant(j)->GetBaryonNumber();
275
276 rha[i][j] = ibry*jbry*rh1;
277 rha[j][i] = rha[i][j];
278
279 // Coulomb terms
280
281 G4double rrs2 = rr2[i][j] + epscl;
282 G4double rrs = std::sqrt ( rrs2 );
283
284 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
285 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
286
287 G4double xerf = 0.0;
288 // T. K. add this protection. 5.8 is good enough for double
289 if ( rrs*c0sw < 5.8 )
290 {
291 #if defined WIN32-VC
292 xerf = CLHEP::HepStat::erf ( rrs*c0sw );
293 #else
294 xerf = std::erf ( rrs*c0sw );
295 #endif
296 }
297 else
298 {
299 xerf = 1.0;
300 }
301
302 G4double erfij = xerf/rrs;
303
304 rhe[i][j] = icharge*jcharge * erfij;
305 rhe[j][i] = rhe[i][j];
306 rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2;
307 rhc[j][i] = rhc[i][j];
308 }
309 }
310
311 void G4QMDMeanField::CalGraduate()
312 {
313 ffr.resize( system->GetTotalNumberOfParticipant() );
314 ffp.resize( system->GetTotalNumberOfParticipant() );
315 rh3d.resize( system->GetTotalNumberOfParticipant() );
316
317 for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; ++i )
318 {
319 G4double rho3 = 0.0;
320 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j )
321 {
322 rho3 += rha[j][i];
323 }
324 rh3d[i] = G4Pow::GetInstance()->powA ( rho3 , pag );
325 }
326
327 for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; ++i )
328 {
329 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
330 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum();
331
332 G4ThreeVector betai = p4i.v()/p4i.e();
333
334 // R-JQMD
335 G4double Vi = GetPotential( i );
336 G4double p_zero = std::sqrt( p4i.e()*p4i.e() + 2*p4i.m()*Vi);
337 G4ThreeVector betai_R = p4i.v()/p_zero;
338 G4double mi_R = p4i.m()/p_zero;
339
340 ffr[i] = betai_R;
341 ffp[i] = G4ThreeVector( 0.0 );
342
343 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j )
344 {
345 G4ThreeVector rj = system->GetParticipant( j )->GetPosition();
346 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
347
348 G4double eij = p4i.e() + p4j.e();
349
350 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
351 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
352
353 G4int inuc = system->GetParticipant(i)->GetNuc();
354 G4int jnuc = system->GetParticipant(j)->GetNuc();
355
356 G4double ccpp = c0g * rha[j][i]
357 + c3g * rha[j][i] * ( rh3d[j] + rh3d[i] )
358 + csg * rha[j][i] * jnuc * inuc
359 * ( 1. - 2. * std::abs( jcharge - icharge ) )
360 + cl * rhc[j][i];
361 ccpp *= mi_R;
362
363 G4double grbb = - rbij[j][i];
364 G4double ccrr = grbb * ccpp / eij;
365
366 G4ThreeVector rij = ri - rj;
367 G4ThreeVector betaij = ( p4i + p4j ).v()/eij;
368 G4ThreeVector cij = betaij - betai;
369
370 ffr[i] = ffr[i] + 2*ccrr* ( rij + grbb*cij );
371 ffp[i] = ffp[i] - 2*ccpp* ( rij + grbb*betaij );
372 }
373 }
374 }
375
376 G4double G4QMDMeanField::GetPotential( G4int i )
377 {
378 G4int n = system->GetTotalNumberOfParticipant();
379
380 G4double rhoa = 0.0;
381 G4double rho3 = 0.0;
382 G4double rhos = 0.0;
383 G4double rhoc = 0.0;
384
385 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
386 G4int inuc = system->GetParticipant(i)->GetNuc();
387
388 for ( G4int j = 0 ; j < n ; ++j )
389 {
390 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
391 G4int jnuc = system->GetParticipant(j)->GetNuc();
392
393 rhoa += rha[j][i];
394 rhoc += rhe[j][i];
395 rhos += rha[j][i] * jnuc * inuc
396 * ( 1 - 2 * std::abs ( jcharge - icharge ) );
397 }
398
399 rho3 = G4Pow::GetInstance()->powA ( rhoa , gamm );
400
401 // return potential
402 return c0 * rhoa + c3 * rho3 + cs * rhos + cl * rhoc;
403 }
404
405 G4double G4QMDMeanField::GetTotalPotential()
406 {
407 G4int n = system->GetTotalNumberOfParticipant();
408
409 std::vector < G4double > rhoa ( n , 0.0 );
410 std::vector < G4double > rho3 ( n , 0.0 );
411 std::vector < G4double > rhos ( n , 0.0 );
412 std::vector < G4double > rhoc ( n , 0.0 );
413
414 for ( G4int i = 0 ; i < n ; ++i )
415 {
416 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
417 G4int inuc = system->GetParticipant(i)->GetNuc();
418
419 for ( G4int j = 0 ; j < n ; ++j )
420 {
421 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
422 G4int jnuc = system->GetParticipant(j)->GetNuc();
423
424 rhoa[i] += rha[j][i];
425 rhoc[i] += rhe[j][i];
426 rhos[i] += rha[j][i] * jnuc * inuc
427 * ( 1 - 2 * std::abs ( jcharge - icharge ) );
428 }
429
430 rho3[i] = G4Pow::GetInstance()->powA ( rhoa[i] , gamm );
431 }
432
433 // return potential
434 return c0 * std::accumulate( rhoa.cbegin() , rhoa.cend() , 0.0 )
435 + c3 * std::accumulate( rho3.cbegin() , rho3.cend() , 0.0 )
436 + cs * std::accumulate( rhos.cbegin() , rhos.cend() , 0.0 )
437 + cl * std::accumulate( rhoc.cbegin() , rhoc.cend() , 0.0 );
438 }
439
440 G4double G4QMDMeanField::calPauliBlockingFactor( G4int i )
441 {
442 // i is supposed beyond total number of Participant()
443
444 G4double pf = 0.0;
445 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
446
447 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j )
448 {
449 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
450 G4int jnuc = system->GetParticipant(j)->GetNuc();
451
452 if ( jcharge == icharge && jnuc == 1 )
453 {
454 G4double expa = -rr2[i][j]*cpw;
455 if ( expa > epsx )
456 {
457 expa = expa - pp2[i][j]*cph;
458 if ( expa > epsx )
459 {
460 pf = pf + G4Exp ( expa );
461 }
462 }
463 }
464 }
465
466 return ( pf - 1.0 ) * cpc;
467 }
468
469 G4bool G4QMDMeanField::IsPauliBlocked( G4int i )
470 {
471 G4bool result = false;
472
473 if ( system->GetParticipant( i )->GetNuc() == 1 )
474 {
475 G4double pf = calPauliBlockingFactor( i );
476 G4double rand = G4UniformRand();
477 if ( pf > rand ) { result = true; }
478 }
479
480 return result;
481 }
482
483 void G4QMDMeanField::DoPropagation( G4double dt )
484 {
485 G4double cc2 = 1.0;
486 G4double cc1 = 1.0 - cc2;
487 G4double cc3 = 1.0 / 2.0 / cc2;
488
489 G4double dt3 = dt * cc3;
490 G4double dt1 = dt * ( cc1 - cc3 );
491 G4double dt2 = dt * cc2;
492
493 CalGraduate();
494
495 G4int n = system->GetTotalNumberOfParticipant();
496
497 // 1st Step
498
499 std::vector< G4ThreeVector > f0r, f0p;
500 f0r.resize( n );
501 f0p.resize( n );
502
503 for ( G4int i = 0 ; i < n ; ++i )
504 {
505 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
506 G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum();
507
508 ri += dt3* ffr[i];
509 p3i += dt3* ffp[i];
510
511 f0r[i] = ffr[i];
512 f0p[i] = ffp[i];
513
514 system->GetParticipant( i )->SetPosition( ri );
515 system->GetParticipant( i )->SetMomentum( p3i );
516
517 // we do not need set total momentum by ourselves
518 }
519
520 // 2nd Step
521
522 Cal2BodyQuantities();
523 CalGraduate();
524
525 for ( G4int i = 0 ; i < n ; ++i )
526 {
527 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
528 G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum();
529
530 ri += dt1* f0r[i] + dt2* ffr[i];
531 p3i += dt1* f0p[i] + dt2* ffp[i];
532
533 system->GetParticipant( i )->SetPosition( ri );
534 system->GetParticipant( i )->SetMomentum( p3i );
535
536 // we do not need set total momentum by ourselves
537 }
538
539 Cal2BodyQuantities();
540 }
541
542 std::vector< G4QMDNucleus* > G4QMDMeanField::DoClusterJudgment()
543 {
544 Cal2BodyQuantities();
545
546 G4double cpf2 = G4Pow::GetInstance()->A23 ( 1.5 * pi*pi * G4Pow::GetInstance()->powA ( 4.0 * pi * wl , -1.5 ) )
547 * hbc * hbc;
548 G4double rcc2 = rclds*rclds;
549
550 G4int n = system->GetTotalNumberOfParticipant();
551 std::vector < G4double > rhoa;
552 rhoa.resize ( n );
553
554 for ( G4int i = 0 ; i < n ; ++i )
555 {
556 rhoa[i] = 0.0;
557
558 if ( system->GetParticipant( i )->GetBaryonNumber() == 1 )
559 {
560 for ( G4int j = 0 ; j < n ; ++j )
561 {
562 if ( system->GetParticipant( j )->GetBaryonNumber() == 1 )
563 rhoa[i] += rha[i][j];
564 }
565 }
566
567 rhoa[i] = G4Pow::GetInstance()->A13 ( rhoa[i] + 1 );
568 }
569
570 // identification of the cluster
571 std::vector < G4bool > is_already_belong_some_cluster;
572
573 // cluster_id participant_id
574 std::multimap < G4int , G4int > comb_map;
575 std::multimap < G4int , G4int > assign_map;
576 assign_map.clear();
577
578 std::vector < G4int > mascl;
579 std::vector < G4int > num;
580 mascl.resize ( n );
581 num.resize ( n );
582 is_already_belong_some_cluster.resize ( n );
583
584 std::vector < G4int > is_assigned_to ( n , -1 );
585 std::multimap < G4int , G4int > clusters;
586
587 for ( G4int i = 0 ; i < n ; ++i )
588 {
589 mascl[i] = 1;
590 num[i] = 1;
591 is_already_belong_some_cluster[i] = false;
592 }
593
594 G4int ichek = 1;
595 G4int id = 0;
596 G4int cluster_id = -1;
597 for ( G4int i = 0 ; i < n-1 ; ++i )
598 {
599 G4bool hasThisCompany = false;
600
601 if ( system->GetParticipant( i )->GetBaryonNumber() == 1 )
602 {
603 G4int j1 = i + 1;
604 for ( G4int j = j1 ; j < n ; ++j )
605 {
606 std::vector < G4int > cluster_participants;
607 if ( system->GetParticipant( j )->GetBaryonNumber() == 1 )
608 {
609 G4double rdist2 = rr2[ i ][ j ];
610 G4double pdist2 = pp2[ i ][ j ];
611 G4double pcc2 = cpf2
612 * ( rhoa[ i ] + rhoa[ j ] )
613 * ( rhoa[ i ] + rhoa[ j ] );
614
615 // Check phase space: close enough?
616 if ( rdist2 < rcc2 && pdist2 < pcc2 )
617 {
618 if ( is_assigned_to [ j ] == -1 )
619 {
620 if ( is_assigned_to [ i ] == -1 )
621 {
622 if ( clusters.size() != 0 )
623 {
624 id = clusters.rbegin()->first + 1;
625 }
626 else
627 {
628 id = 0;
629 }
630 clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , i ) );
631 is_assigned_to [ i ] = id;
632 clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , j ) );
633 is_assigned_to [ j ] = id;
634 }
635 else
636 {
637 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ i ] , j ) );
638 is_assigned_to [ j ] = is_assigned_to [ i ];
639 }
640 }
641 else
642 {
643 // j is already belong to some cluster
644 if ( is_assigned_to [ i ] == -1 )
645 {
646 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , i ) );
647 is_assigned_to [ i ] = is_assigned_to [ j ];
648 }
649 else
650 {
651 // i has companion
652 if ( is_assigned_to [ i ] != is_assigned_to [ j ] )
653 {
654 // move companions to the cluster
655 std::multimap< G4int , G4int > clusters_tmp;
656 G4int target_cluster_id;
657 if ( is_assigned_to [ i ] > is_assigned_to [ j ] )
658 {
659 target_cluster_id = is_assigned_to [ i ];
660 }
661 else
662 {
663 target_cluster_id = is_assigned_to [ j ];
664 }
665 for ( auto it = clusters.cbegin() ; it != clusters.cend() ; ++it )
666 {
667 if ( it->first == target_cluster_id )
668 {
669 is_assigned_to [ it->second ] = is_assigned_to [ j ];
670 clusters_tmp.insert(std::multimap<G4int,G4int>::value_type(is_assigned_to[j], it->second));
671 }
672 else
673 {
674 clusters_tmp.insert(std::multimap<G4int,G4int>::value_type(it->first, it->second));
675 }
676 }
677 clusters = std::move(clusters_tmp);
678 }
679 }
680 }
681
682 comb_map.insert( std::multimap<G4int,G4int>::value_type ( i , j ) );
683 cluster_participants.push_back ( j );
684
685 if ( assign_map.find( cluster_id ) == assign_map.end() )
686 {
687 is_already_belong_some_cluster[i] = true;
688 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , i ) );
689 hasThisCompany = true;
690 }
691 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , j ) );
692 is_already_belong_some_cluster[j] = true;
693 }
694
695 if ( ichek == i )
696 {
697 ++ichek;
698 }
699 }
700 }
701 }
702 if ( hasThisCompany == true ) { ++cluster_id; }
703 }
704
705 // sort
706 // Heavy cluster comes first
707 // size cluster_id
708 std::multimap< G4int , G4int > sorted_cluster_map;
709 for ( G4int i = 0 ; i <= id ; ++i ) // << "<=" because id is highest cluster nubmer.
710 {
711 sorted_cluster_map.insert ( std::multimap<G4int,G4int>::value_type ( (G4int) clusters.count( i ) , i ) );
712 }
713
714 // create nucleus from divided clusters
715 std::vector < G4QMDNucleus* > result;
716 for ( auto it = sorted_cluster_map.crbegin(); it != sorted_cluster_map.crend(); ++it )
717 {
718 if ( it->first != 0 )
719 {
720 G4QMDNucleus* nucleus = new G4QMDNucleus();
721 for ( auto itt = clusters.cbegin(); itt != clusters.cend(); ++itt )
722 {
723 if ( it->second == itt->first )
724 {
725 nucleus->SetParticipant( system->GetParticipant ( itt->second ) );
726 }
727 }
728 result.push_back( nucleus );
729 }
730 }
731
732 // delete participants from current system
733 for ( auto it = result.cbegin(); it != result.cend(); ++it )
734 {
735 system->SubtractSystem ( *it );
736 }
737
738 return result;
739 }
740
741 void G4QMDMeanField::Update()
742 {
743 SetSystem( system );
744 }
745