| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 = G4QMDParamete
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 Se
74 }
75
76 void G4QMDMeanField::SetSystem ( G4QMDSystem*
77 {
78 system = aSystem;
79
80 G4int n = system->GetTotalNumberOfParticipa
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
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()
129
130 for ( G4int j = 1 ; j < system->GetTotalNum
131 {
132 G4ThreeVector rj = system->GetParticipan
133 G4LorentzVector p4j = system->GetPartici
134
135 for ( G4int i = 0 ; i < j ; ++i )
136 {
137 G4ThreeVector ri = system->GetPartici
138 G4LorentzVector p4i = system->GetPart
139
140 G4ThreeVector rij = ri - rj;
141 G4ThreeVector pij = (p4i - p4j).v();
142 G4LorentzVector p4ij = p4i - p4j;
143 G4ThreeVector bij = ( p4i + p4j ).boo
144 G4double gammaij = ( p4i + p4j ).gamm
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*r
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::GetIn
163 + gamma2_ij * G4Pow::GetIns
164
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
182 G4int jbry = system->GetParticipant(j
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->GetParticipan
193 G4int jcharge = system->GetParticipan
194
195 G4double xerf = 0.0;
196 // T. K. add this protection. 5.8 is
197 if ( rrs*c0sw < 5.8 )
198 {
199 #if defined WIN32-VC
200 xerf = CLHEP::HepStat::erf ( rrs*c
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 * ( - erf
215 rhc[j][i] = rhc[i][j];
216 } // i
217 } // j
218 }
219
220 void G4QMDMeanField::Cal2BodyQuantities( G4int
221 {
222 G4ThreeVector ri = system->GetParticipant(
223 G4LorentzVector p4i = system->GetParticipan
224
225 for ( G4int j = 0 ; j < system->GetTotalNum
226 {
227 if ( j == i ) { continue; }
228
229 G4ThreeVector rj = system->GetParticipan
230 G4LorentzVector p4j = system->GetPartici
231
232 G4ThreeVector rij = ri - rj;
233 G4ThreeVector pij = (p4i - p4j).v();
234 G4LorentzVector p4ij = p4i - p4j;
235 G4ThreeVector bij = ( p4i + p4j ).boostV
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::GetInsta
255 + gamma2_ij * G4Pow::GetInstan
256
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)->
274 G4int jbry = system->GetParticipant(j)->
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
285 G4int jcharge = system->GetParticipant(j
286
287 G4double xerf = 0.0;
288 // T. K. add this protection. 5.8 is goo
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
307 rhc[j][i] = rhc[i][j];
308 }
309 }
310
311 void G4QMDMeanField::CalGraduate()
312 {
313 ffr.resize( system->GetTotalNumberOfPartici
314 ffp.resize( system->GetTotalNumberOfPartici
315 rh3d.resize( system->GetTotalNumberOfPartic
316
317 for ( G4int i = 0 ; i < system->GetTotalNum
318 {
319 G4double rho3 = 0.0;
320 for ( G4int j = 0 ; j < system->GetTotal
321 {
322 rho3 += rha[j][i];
323 }
324 rh3d[i] = G4Pow::GetInstance()->powA ( r
325 }
326
327 for ( G4int i = 0 ; i < system->GetTotalNum
328 {
329 G4ThreeVector ri = system->GetParticipan
330 G4LorentzVector p4i = system->GetPartici
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
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->GetTotal
344 {
345 G4ThreeVector rj = system->GetPartici
346 G4LorentzVector p4j = system->GetPart
347
348 G4double eij = p4i.e() + p4j.e();
349
350 G4int icharge = system->GetParticipan
351 G4int jcharge = system->GetParticipan
352
353 G4int inuc = system->GetParticipant(i
354 G4int jnuc = system->GetParticipant(j
355
356 G4double ccpp = c0g * rha[j][i]
357 + c3g * rha[j][i] * ( r
358 + csg * rha[j][i] * jnu
359 * ( 1. - 2. * std
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 )
368 G4ThreeVector cij = betaij - betai;
369
370 ffr[i] = ffr[i] + 2*ccrr* ( rij + grb
371 ffp[i] = ffp[i] - 2*ccpp* ( rij + grb
372 }
373 }
374 }
375
376 G4double G4QMDMeanField::GetPotential( G4int i
377 {
378 G4int n = system->GetTotalNumberOfParticipa
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)->
386 G4int inuc = system->GetParticipant(i)->Get
387
388 for ( G4int j = 0 ; j < n ; ++j )
389 {
390 G4int jcharge = system->GetParticipant(j
391 G4int jnuc = system->GetParticipant(j)->
392
393 rhoa += rha[j][i];
394 rhoc += rhe[j][i];
395 rhos += rha[j][i] * jnuc * inuc
396 * ( 1 - 2 * std::abs (
397 }
398
399 rho3 = G4Pow::GetInstance()->powA ( rhoa ,
400
401 // return potential
402 return c0 * rhoa + c3 * rho3 + cs * rhos +
403 }
404
405 G4double G4QMDMeanField::GetTotalPotential()
406 {
407 G4int n = system->GetTotalNumberOfParticipa
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
417 G4int inuc = system->GetParticipant(i)->
418
419 for ( G4int j = 0 ; j < n ; ++j )
420 {
421 G4int jcharge = system->GetParticipan
422 G4int jnuc = system->GetParticipant(j
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 ( jcharg
428 }
429
430 rho3[i] = G4Pow::GetInstance()->powA ( r
431 }
432
433 // return potential
434 return c0 * std::accumulate( rhoa.cbegin()
435 + c3 * std::accumulate( rho3.cbeg
436 + cs * std::accumulate( rhos.cbeg
437 + cl * std::accumulate( rhoc.cbeg
438 }
439
440 G4double G4QMDMeanField::calPauliBlockingFacto
441 {
442 // i is supposed beyond total number of Par
443
444 G4double pf = 0.0;
445 G4int icharge = system->GetParticipant(i)->
446
447 for ( G4int j = 0 ; j < system->GetTotalNum
448 {
449 G4int jcharge = system->GetParticipant(j
450 G4int jnuc = system->GetParticipant(j)->
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()
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 d
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->GetTotalNumberOfParticipa
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->GetParticipan
506 G4ThreeVector p3i = system->GetParticipa
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
515 system->GetParticipant( i )->SetMomentum
516
517 // we do not need set total momentum by
518 }
519
520 // 2nd Step
521
522 Cal2BodyQuantities();
523 CalGraduate();
524
525 for ( G4int i = 0 ; i < n ; ++i )
526 {
527 G4ThreeVector ri = system->GetParticipan
528 G4ThreeVector p3i = system->GetParticipa
529
530 ri += dt1* f0r[i] + dt2* ffr[i];
531 p3i += dt1* f0p[i] + dt2* ffp[i];
532
533 system->GetParticipant( i )->SetPosition
534 system->GetParticipant( i )->SetMomentum
535
536 // we do not need set total momentum by
537 }
538
539 Cal2BodyQuantities();
540 }
541
542 std::vector< G4QMDNucleus* > G4QMDMeanField::D
543 {
544 Cal2BodyQuantities();
545
546 G4double cpf2 = G4Pow::GetInstance()->A23 (
547 * hbc * hbc;
548 G4double rcc2 = rclds*rclds;
549
550 G4int n = system->GetTotalNumberOfParticipa
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 )->GetBary
559 {
560 for ( G4int j = 0 ; j < n ; ++j )
561 {
562 if ( system->GetParticipant( j )->Get
563 rhoa[i] += rha[i][j];
564 }
565 }
566
567 rhoa[i] = G4Pow::GetInstance()->A13 ( rho
568 }
569
570 // identification of the cluster
571 std::vector < G4bool > is_already_belong_so
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 ,
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 )->GetBary
602 {
603 G4int j1 = i + 1;
604 for ( G4int j = j1 ; j < n ; ++j )
605 {
606 std::vector < G4int > cluster_partici
607 if ( system->GetParticipant( j )->Get
608 {
609 G4double rdist2 = rr2[ i ][ j ];
610 G4double pdist2 = pp2[ i ][ j ];
611 G4double pcc2 = cpf2
612 * ( rhoa[ i ] + rhoa[
613 * ( rhoa[ i ] + rhoa[
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()->fir
625 }
626 else
627 {
628 id = 0;
629 }
630 clusters.insert ( std::multim
631 is_assigned_to [ i ] = id;
632 clusters.insert ( std::multim
633 is_assigned_to [ j ] = id;
634 }
635 else
636 {
637 clusters.insert ( std::multim
638 is_assigned_to [ j ] = is_ass
639 }
640 }
641 else
642 {
643 // j is already belong to some
644 if ( is_assigned_to [ i ] == -1
645 {
646 clusters.insert ( std::multim
647 is_assigned_to [ i ] = is_ass
648 }
649 else
650 {
651 // i has companion
652 if ( is_assigned_to [ i ] !=
653 {
654 // move companions to the c
655 std::multimap< G4int , G4in
656 G4int target_cluster_id;
657 if ( is_assigned_to [ i ] >
658 {
659 target_cluster_id = is_as
660 }
661 else
662 {
663 target_cluster_id = is_as
664 }
665 for ( auto it = clusters.cb
666 {
667 if ( it->first == target_
668 {
669 is_assigned_to [ it->se
670 clusters_tmp.insert(std
671 }
672 else
673 {
674 clusters_tmp.insert(std
675 }
676 }
677 clusters = std::move(cluste
678 }
679 }
680 }
681
682 comb_map.insert( std::multimap<G4
683 cluster_participants.push_back (
684
685 if ( assign_map.find( cluster_id
686 {
687 is_already_belong_some_cluster[
688 assign_map.insert ( std::multim
689 hasThisCompany = true;
690 }
691 assign_map.insert ( std::multimap
692 is_already_belong_some_cluster[j]
693 }
694
695 if ( ichek == i )
696 {
697 ++ichek;
698 }
699 }
700 }
701 }
702 if ( hasThisCompany == true ) { ++cluster
703 }
704
705 // sort
706 // Heavy cluster comes first
707 // size cluster_id
708 std::multimap< G4int , G4int > sorted_clust
709 for ( G4int i = 0 ; i <= id ; ++i ) // <<
710 {
711 sorted_cluster_map.insert ( std::multimap
712 }
713
714 // create nucleus from divided clusters
715 std::vector < G4QMDNucleus* > result;
716 for ( auto it = sorted_cluster_map.crbegin(
717 {
718 if ( it->first != 0 )
719 {
720 G4QMDNucleus* nucleus = new G4QMDNucleu
721 for ( auto itt = clusters.cbegin(); itt
722 {
723 if ( it->second == itt->first )
724 {
725 nucleus->SetParticipant( system->Ge
726 }
727 }
728 result.push_back( nucleus );
729 }
730 }
731
732 // delete participants from current system
733 for ( auto it = result.cbegin(); it != resu
734 {
735 system->SubtractSystem ( *it );
736 }
737
738 return result;
739 }
740
741 void G4QMDMeanField::Update()
742 {
743 SetSystem( system );
744 }
745