Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4QMDMeanField.cc
Go to the documentation of this file.
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// 081120 Add Update by T. Koi
27//
28
29#include <map>
30#include <algorithm>
31#include <numeric>
32
33#include <cmath>
34#include <CLHEP/Random/Stat.h>
35
36#include "G4QMDMeanField.hh"
37#include "G4QMDParameters.hh"
38#include "G4Exp.hh"
39#include "G4Pow.hh"
40#include "G4PhysicalConstants.hh"
41#include "Randomize.hh"
42
43G4QMDMeanField::G4QMDMeanField()
44: rclds ( 4.0 ) // distance for cluster judgement
45, epsx ( -20.0 ) // gauss term
46, epscl ( 0.0001 ) // coulomb term
47, irelcr ( 1 )
48{
49
50 G4QMDParameters* parameters = G4QMDParameters::GetInstance();
51 wl = parameters->Get_wl();
52 cl = parameters->Get_cl();
53 rho0 = parameters->Get_rho0();
54 hbc = parameters->Get_hbc();
55 gamm = parameters->Get_gamm();
56
57 cpw = parameters->Get_cpw();
58 cph = parameters->Get_cph();
59 cpc = parameters->Get_cpc();
60
61 c0 = parameters->Get_c0();
62 c3 = parameters->Get_c3();
63 cs = parameters->Get_cs();
64
65// distance
66 c0w = 1.0/4.0/wl;
67 //c3w = 1.0/4.0/wl; //no need
68 c0sw = std::sqrt( c0w );
69 clw = 2.0 / std::sqrt ( 4.0 * pi * wl );
70
71// graduate
72 c0g = - c0 / ( 2.0 * wl );
73 c3g = - c3 / ( 4.0 * wl ) * gamm;
74 csg = - cs / ( 2.0 * wl );
75 pag = gamm - 1;
76
77 system = NULL; // will be set through SetSystem method
78}
79
80
81
82G4QMDMeanField::~G4QMDMeanField()
83{
84 ;
85}
86
87
88
89void G4QMDMeanField::SetSystem ( G4QMDSystem* aSystem )
90{
91
92 //std::cout << "QMDMeanField SetSystem" << std::endl;
93
94 system = aSystem;
95
96 G4int n = system->GetTotalNumberOfParticipant();
97
98 pp2.clear();
99 rr2.clear();
100 rbij.clear();
101 rha.clear();
102 rhe.clear();
103 rhc.clear();
104
105 rr2.resize( n );
106 pp2.resize( n );
107 rbij.resize( n );
108 rha.resize( n );
109 rhe.resize( n );
110 rhc.resize( n );
111
112 for ( int i = 0 ; i < n ; i++ )
113 {
114 rr2[i].resize( n );
115 pp2[i].resize( n );
116 rbij[i].resize( n );
117 rha[i].resize( n );
118 rhe[i].resize( n );
119 rhc[i].resize( n );
120 }
121
122
123 ffr.clear();
124 ffp.clear();
125 rh3d.clear();
126
127 ffr.resize( n );
128 ffp.resize( n );
129 rh3d.resize( n );
130
131 Cal2BodyQuantities();
132
133}
134
135void G4QMDMeanField::SetNucleus ( G4QMDNucleus* aNucleus )
136{
137
138 //std::cout << "QMDMeanField SetNucleus" << std::endl;
139
140 SetSystem( aNucleus );
141
142 G4double totalPotential = GetTotalPotential();
143 aNucleus->SetTotalPotential( totalPotential );
144
145 aNucleus->CalEnergyAndAngularMomentumInCM();
146
147}
148
149
150
151void G4QMDMeanField::Cal2BodyQuantities()
152{
153
154 if ( system->GetTotalNumberOfParticipant() < 2 ) return;
155
156 for ( G4int j = 1 ; j < system->GetTotalNumberOfParticipant() ; j++ )
157 {
158
159 G4ThreeVector rj = system->GetParticipant( j )->GetPosition();
160 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
161
162 for ( G4int i = 0 ; i < j ; i++ )
163 {
164
165 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
166 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum();
167
168 G4ThreeVector rij = ri - rj;
169 G4ThreeVector pij = (p4i - p4j).v();
170 G4LorentzVector p4ij = p4i - p4j;
171 G4ThreeVector bij = ( p4i + p4j ).boostVector();
172 G4double gammaij = ( p4i + p4j ).gamma();
173
174 G4double eij = ( p4i + p4j ).e();
175
176 G4double rbrb = rij*bij;
177// G4double bij2 = bij*bij;
178 G4double rij2 = rij*rij;
179 G4double pij2 = pij*pij;
180
181 rbrb = irelcr * rbrb;
182 G4double gamma2_ij = gammaij*gammaij;
183
184
185 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb;
186 rr2[j][i] = rr2[i][j];
187
188 rbij[i][j] = gamma2_ij * rbrb;
189 rbij[j][i] = - rbij[i][j];
190
191 pp2[i][j] = pij2
192 + irelcr * ( - G4Pow::GetInstance()->powN ( p4i.e() - p4j.e() , 2 )
193 + gamma2_ij * G4Pow::GetInstance()->powN ( ( ( p4i.m2() - p4j.m2() ) / eij ) , 2 ) );
194
195
196 pp2[j][i] = pp2[i][j];
197
198// Gauss term
199
200 G4double expa1 = - rr2[i][j] * c0w;
201
202 G4double rh1;
203 if ( expa1 > epsx )
204 {
205 rh1 = G4Exp( expa1 );
206 }
207 else
208 {
209 rh1 = 0.0;
210 }
211
212 G4int ibry = system->GetParticipant(i)->GetBaryonNumber();
213 G4int jbry = system->GetParticipant(j)->GetBaryonNumber();
214
215
216 rha[i][j] = ibry*jbry*rh1;
217 rha[j][i] = rha[i][j];
218
219// Coulomb terms
220
221 G4double rrs2 = rr2[i][j] + epscl;
222 G4double rrs = std::sqrt ( rrs2 );
223
224 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
225 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
226
227 G4double xerf = 0.0;
228 // T. K. add this protection. 5.8 is good enough for double
229 if ( rrs*c0sw < 5.8 ) {
230 //erf = G4RandStat::erf ( rrs*c0sw );
231 //Restore to CLHEP for avoiding compilation error in MT
232 //erf = CLHEP::HepStat::erf ( rrs*c0sw );
233 //Use cmath
234#if defined WIN32-VC
235 xerf = CLHEP::HepStat::erf ( rrs*c0sw );
236#else
237 xerf = erf ( rrs*c0sw );
238#endif
239 } else {
240 xerf = 1.0;
241 }
242
243 G4double erfij = xerf/rrs;
244
245
246 rhe[i][j] = icharge*jcharge * erfij;
247
248 rhe[j][i] = rhe[i][j];
249
250 rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2;
251
252 rhc[j][i] = rhc[i][j];
253
254 } // i
255 } // j
256}
257
258
259
260void G4QMDMeanField::Cal2BodyQuantities( G4int i )
261{
262
263 //std::cout << "Cal2BodyQuantities " << i << std::endl;
264
265 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
266 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum();
267
268
269 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; j ++ )
270 {
271 if ( j == i ) continue;
272
273 G4ThreeVector rj = system->GetParticipant( j )->GetPosition();
274 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
275
276 G4ThreeVector rij = ri - rj;
277 G4ThreeVector pij = (p4i - p4j).v();
278 G4LorentzVector p4ij = p4i - p4j;
279 G4ThreeVector bij = ( p4i + p4j ).boostVector();
280 G4double gammaij = ( p4i + p4j ).gamma();
281
282 G4double eij = ( p4i + p4j ).e();
283
284 G4double rbrb = rij*bij;
285// G4double bij2 = bij*bij;
286 G4double rij2 = rij*rij;
287 G4double pij2 = pij*pij;
288
289 rbrb = irelcr * rbrb;
290 G4double gamma2_ij = gammaij*gammaij;
291
292/*
293 G4double rbrb = 0.0;
294 G4double beta2_ij = 0.0;
295 G4double rij2 = 0.0;
296 G4double pij2 = 0.0;
297
298//
299 G4LorentzVector p4ip4j = p4i + p4j;
300 G4double eij = p4ip4j.e();
301
302 G4ThreeVector r = ri - rj;
303 G4LorentzVector p4 = p4i - p4j;
304
305 rbrb = r.x()*p4ip4j.x()/eij
306 + r.y()*p4ip4j.y()/eij
307 + r.z()*p4ip4j.z()/eij;
308
309 beta2_ij = ( p4ip4j.x()*p4ip4j.x() + p4ip4j.y()*p4ip4j.y() + p4ip4j.z()*p4ip4j.z() ) / ( eij*eij );
310 rij2 = r*r;
311 pij2 = p4.v()*p4.v();
312
313 rbrb = irelcr * rbrb;
314
315 G4double gamma2_ij = 1 / ( 1 - beta2_ij );
316*/
317
318 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb;
319 rr2[j][i] = rr2[i][j];
320
321 rbij[i][j] = gamma2_ij * rbrb;
322 rbij[j][i] = - rbij[i][j];
323
324 pp2[i][j] = pij2
325 + irelcr * ( - G4Pow::GetInstance()->powN ( p4i.e() - p4j.e() , 2 )
326 + gamma2_ij * G4Pow::GetInstance()->powN ( ( ( p4i.m2() - p4j.m2() ) / eij ) , 2 ) );
327
328 pp2[j][i] = pp2[i][j];
329
330// Gauss term
331
332 G4double expa1 = - rr2[i][j] * c0w;
333
334 G4double rh1;
335 if ( expa1 > epsx )
336 {
337 rh1 = G4Exp( expa1 );
338 }
339 else
340 {
341 rh1 = 0.0;
342 }
343
344 G4int ibry = system->GetParticipant(i)->GetBaryonNumber();
345 G4int jbry = system->GetParticipant(j)->GetBaryonNumber();
346
347
348 rha[i][j] = ibry*jbry*rh1;
349 rha[j][i] = rha[i][j];
350
351// Coulomb terms
352
353 G4double rrs2 = rr2[i][j] + epscl;
354 G4double rrs = std::sqrt ( rrs2 );
355
356 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
357 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
358
359 G4double xerf = 0.0;
360 // T. K. add this protection. 5.8 is good enough for double
361 if ( rrs*c0sw < 5.8 ) {
362 //xerf = G4RandStat::erf ( rrs*c0sw );
363 //Use cmath
364#if defined WIN32-VC
365 xerf = CLHEP::HepStat::erf ( rrs*c0sw );
366#else
367 xerf = erf ( rrs*c0sw );
368#endif
369 } else {
370 xerf = 1.0;
371 }
372
373 G4double erfij = xerf/rrs;
374
375
376 rhe[i][j] = icharge*jcharge * erfij;
377
378 rhe[j][i] = rhe[i][j];
379
380// G4double clw;
381
382 rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2;
383
384 rhc[j][i] = rhc[i][j];
385
386 }
387
388}
389
390
391
392void G4QMDMeanField::CalGraduate()
393{
394
395 ffr.resize( system->GetTotalNumberOfParticipant() );
396 ffp.resize( system->GetTotalNumberOfParticipant() );
397 rh3d.resize( system->GetTotalNumberOfParticipant() );
398
399 for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; i ++ )
400 {
401 G4double rho3 = 0.0;
402 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; j ++ )
403 {
404 rho3 += rha[j][i];
405 }
406 rh3d[i] = G4Pow::GetInstance()->powA ( rho3 , pag );
407 }
408
409
410 for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; i ++ )
411 {
412
413 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
414 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum();
415
416 G4ThreeVector betai = p4i.v()/p4i.e();
417
418// R-JQMD
419 G4double Vi = GetPotential( i );
420 G4double p_zero = std::sqrt( p4i.e()*p4i.e() + 2*p4i.m()*Vi);
421 G4ThreeVector betai_R = p4i.v()/p_zero;
422 G4double mi_R = p4i.m()/p_zero;
423//
424 ffr[i] = betai_R;
425 ffp[i] = G4ThreeVector( 0.0 );
426
427 if ( false )
428 {
429 ffr[i] = betai;
430 mi_R = 1.0;
431 }
432
433 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; j ++ )
434 {
435
436 G4ThreeVector rj = system->GetParticipant( j )->GetPosition();
437 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
438
439 G4double eij = p4i.e() + p4j.e();
440
441 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
442 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
443
444 G4int inuc = system->GetParticipant(i)->GetNuc();
445 G4int jnuc = system->GetParticipant(j)->GetNuc();
446
447 G4double ccpp = c0g * rha[j][i]
448 + c3g * rha[j][i] * ( rh3d[j] + rh3d[i] )
449 + csg * rha[j][i] * jnuc * inuc
450 * ( 1. - 2. * std::abs( jcharge - icharge ) )
451 + cl * rhc[j][i];
452 ccpp *= mi_R;
453
454/*
455 G4cout << c0g << " " << c3g << " " << csg << " " << cl << G4endl;
456 G4cout << "ccpp " << i << " " << j << " " << ccpp << G4endl;
457 G4cout << "rha[j][i] " << rha[j][i] << G4endl;
458 G4cout << "rh3d " << rh3d[j] << " " << rh3d[i] << G4endl;
459 G4cout << "rhc[j][i] " << rhc[j][i] << G4endl;
460*/
461
462 G4double grbb = - rbij[j][i];
463 G4double ccrr = grbb * ccpp / eij;
464
465/*
466 G4cout << "ccrr " << ccrr << G4endl;
467 G4cout << "grbb " << grbb << G4endl;
468*/
469
470
471 G4ThreeVector rij = ri - rj;
472 G4ThreeVector betaij = ( p4i + p4j ).v()/eij;
473
474 G4ThreeVector cij = betaij - betai;
475
476 ffr[i] = ffr[i] + 2*ccrr* ( rij + grbb*cij );
477
478 ffp[i] = ffp[i] - 2*ccpp* ( rij + grbb*betaij );
479
480 }
481 }
482
483 //std::cout << "gradu 0 " << ffr[0] << " " << ffp[0] << std::endl;
484 //std::cout << "gradu 1 " << ffr[1] << " " << ffp[1] << std::endl;
485
486}
487
488
489
490G4double G4QMDMeanField::GetPotential( G4int i )
491{
492 G4int n = system->GetTotalNumberOfParticipant();
493
494 G4double rhoa = 0.0;
495 G4double rho3 = 0.0;
496 G4double rhos = 0.0;
497 G4double rhoc = 0.0;
498
499
500 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
501 G4int inuc = system->GetParticipant(i)->GetNuc();
502
503 for ( G4int j = 0 ; j < n ; j ++ )
504 {
505 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
506 G4int jnuc = system->GetParticipant(j)->GetNuc();
507
508 rhoa += rha[j][i];
509 rhoc += rhe[j][i];
510 rhos += rha[j][i] * jnuc * inuc
511 * ( 1 - 2 * std::abs ( jcharge - icharge ) );
512 }
513
514 rho3 = G4Pow::GetInstance()->powA ( rhoa , gamm );
515
516 G4double potential = c0 * rhoa
517 + c3 * rho3
518 + cs * rhos
519 + cl * rhoc;
520
521 return potential;
522}
523
524
525
526G4double G4QMDMeanField::GetTotalPotential()
527{
528
529 G4int n = system->GetTotalNumberOfParticipant();
530
531 std::vector < G4double > rhoa ( n , 0.0 );
532 std::vector < G4double > rho3 ( n , 0.0 );
533 std::vector < G4double > rhos ( n , 0.0 );
534 std::vector < G4double > rhoc ( n , 0.0 );
535
536
537 for ( G4int i = 0 ; i < n ; i ++ )
538 {
539 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
540 G4int inuc = system->GetParticipant(i)->GetNuc();
541
542 for ( G4int j = 0 ; j < n ; j ++ )
543 {
544 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
545 G4int jnuc = system->GetParticipant(j)->GetNuc();
546
547 rhoa[i] += rha[j][i];
548 rhoc[i] += rhe[j][i];
549 rhos[i] += rha[j][i] * jnuc * inuc
550 * ( 1 - 2 * std::abs ( jcharge - icharge ) );
551 }
552
553 rho3[i] = G4Pow::GetInstance()->powA ( rhoa[i] , gamm );
554 }
555
556 G4double potential = c0 * std::accumulate( rhoa.begin() , rhoa.end() , 0.0 )
557 + c3 * std::accumulate( rho3.begin() , rho3.end() , 0.0 )
558 + cs * std::accumulate( rhos.begin() , rhos.end() , 0.0 )
559 + cl * std::accumulate( rhoc.begin() , rhoc.end() , 0.0 );
560
561 return potential;
562
563}
564
565
566
567G4double G4QMDMeanField::calPauliBlockingFactor( G4int i )
568{
569
570 G4double pf = 0.0;
571// i is supposed beyond total number of Participant()
572 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
573
574 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; j++ )
575 {
576
577 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
578 G4int jnuc = system->GetParticipant(j)->GetNuc();
579
580 if ( jcharge == icharge && jnuc == 1 )
581 {
582
583/*
584 G4cout << "Pauli i j " << i << " " << j << G4endl;
585 G4cout << "Pauli icharge " << icharge << G4endl;
586 G4cout << "Pauli jcharge " << jcharge << G4endl;
587*/
588 G4double expa = -rr2[i][j]*cpw;
589
590
591 if ( expa > epsx )
592 {
593 expa = expa - pp2[i][j]*cph;
594/*
595 G4cout << "Pauli cph " << cph << G4endl;
596 G4cout << "Pauli pp2 " << pp2[i][j] << G4endl;
597 G4cout << "Pauli expa " << expa << G4endl;
598 G4cout << "Pauli epsx " << epsx << G4endl;
599*/
600 if ( expa > epsx )
601 {
602// std::cout << "Pauli phase " << pf << std::endl;
603 pf = pf + G4Exp ( expa );
604 }
605 }
606 }
607
608 }
609
610
611 pf = ( pf - 1.0 ) * cpc;
612
613 //std::cout << "Pauli pf " << pf << std::endl;
614
615 return pf;
616
617}
618
619
620
621G4bool G4QMDMeanField::IsPauliBlocked( G4int i )
622{
623 G4bool result = false;
624
625 if ( system->GetParticipant( i )->GetNuc() == 1 )
626 {
627 G4double pf = calPauliBlockingFactor( i );
628 G4double rand = G4UniformRand();
629 if ( pf > rand ) result = true;
630 }
631
632 return result;
633}
634
635
636
637void G4QMDMeanField::DoPropagation( G4double dt )
638{
639
640 G4double cc2 = 1.0;
641 G4double cc1 = 1.0 - cc2;
642 G4double cc3 = 1.0 / 2.0 / cc2;
643
644 G4double dt3 = dt * cc3;
645 G4double dt1 = dt * ( cc1 - cc3 );
646 G4double dt2 = dt * cc2;
647
648 CalGraduate();
649
650 G4int n = system->GetTotalNumberOfParticipant();
651
652// 1st Step
653
654 std::vector< G4ThreeVector > f0r, f0p;
655 f0r.resize( n );
656 f0p.resize( n );
657
658 for ( G4int i = 0 ; i < n ; i++ )
659 {
660 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
661 G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum();
662
663 ri += dt3* ffr[i];
664 p3i += dt3* ffp[i];
665
666 f0r[i] = ffr[i];
667 f0p[i] = ffp[i];
668
669 system->GetParticipant( i )->SetPosition( ri );
670 system->GetParticipant( i )->SetMomentum( p3i );
671
672// we do not need set total momentum by ourselvs
673 }
674
675// 2nd Step
676 Cal2BodyQuantities();
677 CalGraduate();
678
679 for ( G4int i = 0 ; i < n ; i++ )
680 {
681 G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
682 G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum();
683
684 ri += dt1* f0r[i] + dt2* ffr[i];
685 p3i += dt1* f0p[i] + dt2* ffp[i];
686
687 system->GetParticipant( i )->SetPosition( ri );
688 system->GetParticipant( i )->SetMomentum( p3i );
689
690// we do not need set total momentum by ourselvs
691 }
692
693 Cal2BodyQuantities();
694
695}
696
697
698
699std::vector< G4QMDNucleus* > G4QMDMeanField::DoClusterJudgment()
700{
701
702 //std::cout << "MeanField DoClusterJudgemnt" << std::endl;
703
704 Cal2BodyQuantities();
705
706 G4double cpf2 = G4Pow::GetInstance()->A23 ( 1.5 * pi*pi * G4Pow::GetInstance()->powA ( 4.0 * pi * wl , -1.5 ) ) * hbc * hbc;
707 G4double rcc2 = rclds*rclds;
708
709 G4int n = system->GetTotalNumberOfParticipant();
710 std::vector < G4double > rhoa;
711 rhoa.resize ( n );
712
713 for ( G4int i = 0 ; i < n ; i++ )
714 {
715 rhoa[i] = 0.0;
716
717 if ( system->GetParticipant( i )->GetBaryonNumber() == 1 )
718 {
719 for ( G4int j = 0 ; j < n ; j++ )
720 {
721 if ( system->GetParticipant( j )->GetBaryonNumber() == 1 )
722 rhoa[i] += rha[i][j];
723 }
724 }
725
726 rhoa[i] = G4Pow::GetInstance()->A13 ( rhoa[i] + 1 );
727
728 }
729
730// identification of the cluster
731
732 std::map < G4int , std::vector < G4int > > cluster_map;
733 std::vector < G4bool > is_already_belong_some_cluster;
734
735 // cluster_id participant_id
736 std::multimap < G4int , G4int > comb_map;
737 std::multimap < G4int , G4int > assign_map;
738 assign_map.clear();
739
740 std::vector < G4int > mascl;
741 std::vector < G4int > num;
742 mascl.resize ( n );
743 num.resize ( n );
744 is_already_belong_some_cluster.resize ( n );
745
746 std::vector < G4int > is_assigned_to ( n , -1 );
747 std::multimap < G4int , G4int > clusters;
748
749 for ( G4int i = 0 ; i < n ; i++ )
750 {
751 mascl[i] = 1;
752 num[i] = 1;
753
754 is_already_belong_some_cluster[i] = false;
755 }
756
757
758 G4int ichek = 1;
759
760
761 G4int id = 0;
762 G4int cluster_id = -1;
763 for ( G4int i = 0 ; i < n-1 ; i++ )
764 {
765
766 G4bool hasThisCompany = false;
767// Check only for bryons?
768// std::cout << "Check Baryon " << i << std::endl;
769
770 if ( system->GetParticipant( i )->GetBaryonNumber() == 1 )
771 {
772
773// if ( is_already_belong_some_cluster[i] != true )
774// {
775 //G4int j1 = ichek + 1;
776 G4int j1 = i + 1;
777 for ( G4int j = j1 ; j < n ; j++ )
778 {
779
780 std::vector < G4int > cluster_participants;
781 if ( system->GetParticipant( j )->GetBaryonNumber() == 1 )
782 {
783 G4double rdist2 = rr2[ i ][ j ];
784 G4double pdist2 = pp2[ i ][ j ];
785 //G4double rdist2 = rr2[ num[i] ][ num[j] ];
786 //G4double pdist2 = pp2[ num[i] ][ num[j] ];
787 G4double pcc2 = cpf2
788 * ( rhoa[ i ] + rhoa[ j ] )
789 * ( rhoa[ i ] + rhoa[ j ] );
790
791// Check phase space: close enough?
792 if ( rdist2 < rcc2 && pdist2 < pcc2 )
793 {
794
795/*
796 G4cout << "G4QMDRESULT "
797 << i << " " << j << " " << id << " "
798 << is_assigned_to [ i ] << " " << is_assigned_to [ j ]
799 << G4endl;
800*/
801
802 if ( is_assigned_to [ j ] == -1 )
803 {
804 if ( is_assigned_to [ i ] == -1 )
805 {
806 if ( clusters.size() != 0 )
807 {
808 id = clusters.rbegin()->first + 1;
809 //std::cout << "id is increare " << id << std::endl;
810 }
811 else
812 {
813 id = 0;
814 }
815 clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , i ) );
816 is_assigned_to [ i ] = id;
817 clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , j ) );
818 is_assigned_to [ j ] = id;
819 }
820 else
821 {
822 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ i ] , j ) );
823 is_assigned_to [ j ] = is_assigned_to [ i ];
824 }
825 }
826 else
827 {
828// j is already belong to some cluester
829 if ( is_assigned_to [ i ] == -1 )
830 {
831 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , i ) );
832 is_assigned_to [ i ] = is_assigned_to [ j ];
833 }
834 else
835 {
836// i has companion
837 if ( is_assigned_to [ i ] != is_assigned_to [ j ] )
838 {
839// move companions to the cluster
840//
841 //std::cout << "combine " << is_assigned_to [ i ] << " to " << is_assigned_to [ j ] << std::endl;
842 std::multimap< G4int , G4int > clusters_tmp;
843 G4int target_cluster_id;
844 if ( is_assigned_to [ i ] > is_assigned_to [ j ] )
845 target_cluster_id = is_assigned_to [ i ];
846 else
847 target_cluster_id = is_assigned_to [ j ];
848
849 for ( std::multimap< G4int , G4int >::iterator it
850 = clusters.begin() ; it != clusters.end() ; it++ )
851 {
852
853 //std::cout << it->first << " " << it->second << " " << target_cluster_id << std::endl;
854 if ( it->first == target_cluster_id )
855 {
856 //std::cout << "move " << it->first << " " << it->second << std::endl;
857 is_assigned_to [ it->second ] = is_assigned_to [ j ];
858 clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , it->second ) );
859 }
860 else
861 {
862 clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( it->first , it->second ) );
863 }
864 }
865
866 clusters = clusters_tmp;
867 //id = clusters.rbegin()->first;
868 //id = target_cluster_id;
869 //std::cout << "id " << id << std::endl;
870 }
871 }
872 }
873
874 //std::cout << "combination " << i << " " << j << std::endl;
875 comb_map.insert( std::multimap<G4int,G4int>::value_type ( i , j ) );
876 cluster_participants.push_back ( j );
877
878
879
880 if ( assign_map.find( cluster_id ) == assign_map.end() )
881 {
882 is_already_belong_some_cluster[i] = true;
883 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , i ) );
884 hasThisCompany = true;
885 }
886 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , j ) );
887 is_already_belong_some_cluster[j] = true;
888
889 }
890
891 if ( ichek == i )
892 {
893 ichek++;
894 }
895 }
896
897 if ( cluster_participants.size() > 0 )
898 {
899// cluster , participant
900 cluster_map.insert ( std::pair < G4int , std::vector < G4int > > ( i , cluster_participants ) );
901 }
902 }
903// }
904 }
905 if ( hasThisCompany == true ) cluster_id++;
906 }
907
908 //std::cout << " id " << id << std::endl;
909
910// sort
911// Heavy cluster comes first
912// size cluster_id
913 std::multimap< G4int , G4int > sorted_cluster_map;
914 for ( G4int i = 0 ; i <= id ; i++ ) // << "<=" because id is highest cluster nubmer.
915 {
916
917 //std::cout << i << " cluster has " << clusters.count( i ) << " nucleons." << std::endl;
918 sorted_cluster_map.insert ( std::multimap<G4int,G4int>::value_type ( (G4int) clusters.count( i ) , i ) );
919
920 }
921
922
923// create nucleus from devided clusters
924 std::vector < G4QMDNucleus* > result;
925 for ( std::multimap < G4int , G4int >::reverse_iterator it
926 = sorted_cluster_map.rbegin() ; it != sorted_cluster_map.rend() ; it ++)
927 {
928
929 //G4cout << "Add Participants to cluseter " << it->second << G4endl;
930
931 if ( it->first != 0 )
932 {
933 G4QMDNucleus* nucleus = new G4QMDNucleus();
934 for ( std::multimap < G4int , G4int >::iterator itt
935 = clusters.begin() ; itt != clusters.end() ; itt ++)
936 {
937
938 if ( it->second == itt->first )
939 {
940 nucleus->SetParticipant( system->GetParticipant ( itt->second ) );
941 //G4cout << "Add Participants " << itt->second << " " << system->GetParticipant ( itt->second )->GetPosition() << G4endl;
942 }
943
944 }
945 result.push_back( nucleus );
946 }
947
948 }
949
950// delete participants from current system
951
952 for ( std::vector < G4QMDNucleus* > ::iterator it
953 = result.begin() ; it != result.end() ; it++ )
954 {
955 system->SubtractSystem ( *it );
956 }
957
958 return result;
959
960}
961
962
963
964void G4QMDMeanField::Update()
965{
966 SetSystem( system );
967}
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:36
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
CLHEP::HepLorentzVector::v
Hep3Vector v() const
CLHEP::HepStat::erf
static double erf(double x)
Definition: flatToGaussian.cc:268
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4Pow::A13
G4double A13(G4double A) const
Definition: G4Pow.cc:116
G4Pow::powN
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:162
G4Pow::powA
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:230
G4Pow::A23
G4double A23(G4double A) const
Definition: G4Pow.hh:131
G4QMDMeanField::GetTotalPotential
G4double GetTotalPotential()
Definition: G4QMDMeanField.cc:526
G4QMDMeanField::SetNucleus
void SetNucleus(G4QMDNucleus *aSystem)
Definition: G4QMDMeanField.cc:135
G4QMDMeanField::DoPropagation
void DoPropagation(G4double)
Definition: G4QMDMeanField.cc:637
G4QMDMeanField::Cal2BodyQuantities
void Cal2BodyQuantities()
Definition: G4QMDMeanField.cc:151
G4QMDMeanField::GetPotential
G4double GetPotential(G4int)
Definition: G4QMDMeanField.cc:490
G4QMDMeanField::DoClusterJudgment
std::vector< G4QMDNucleus * > DoClusterJudgment()
Definition: G4QMDMeanField.cc:699
G4QMDMeanField::SetSystem
void SetSystem(G4QMDSystem *aSystem)
Definition: G4QMDMeanField.cc:89
G4QMDMeanField::IsPauliBlocked
G4bool IsPauliBlocked(G4int)
Definition: G4QMDMeanField.cc:621
G4QMDNucleus::SetTotalPotential
void SetTotalPotential(G4double x)
Definition: G4QMDNucleus.hh:62
G4QMDNucleus::CalEnergyAndAngularMomentumInCM
void CalEnergyAndAngularMomentumInCM()
Definition: G4QMDNucleus.cc:138
G4QMDParameters::Get_cpc
G4double Get_cpc()
Definition: G4QMDParameters.hh:65
G4QMDParameters::GetInstance
static G4QMDParameters * GetInstance()
Definition: G4QMDParameters.hh:49
G4QMDParameters::Get_hbc
G4double Get_hbc()
Definition: G4QMDParameters.hh:59
G4QMDParameters::Get_rho0
G4double Get_rho0()
Definition: G4QMDParameters.hh:60
G4QMDParameters::Get_gamm
G4double Get_gamm()
Definition: G4QMDParameters.hh:61
G4QMDParameters::Get_cpw
G4double Get_cpw()
Definition: G4QMDParameters.hh:62
G4QMDParameters::Get_cph
G4double Get_cph()
Definition: G4QMDParameters.hh:63
G4QMDParticipant::GetPosition
G4ThreeVector GetPosition()
Definition: G4QMDParticipant.hh:59
G4QMDParticipant::SetPosition
void SetPosition(G4ThreeVector r)
Definition: G4QMDParticipant.hh:58
G4QMDParticipant::Get4Momentum
G4LorentzVector Get4Momentum()
Definition: G4QMDParticipant.cc:50
G4QMDParticipant::GetChargeInUnitOfEplus
G4int GetChargeInUnitOfEplus()
Definition: G4QMDParticipant.hh:73
G4QMDParticipant::GetMomentum
G4ThreeVector GetMomentum()
Definition: G4QMDParticipant.hh:62
G4QMDParticipant::SetMomentum
void SetMomentum(G4ThreeVector p)
Definition: G4QMDParticipant.hh:61
G4QMDSystem::GetParticipant
G4QMDParticipant * GetParticipant(G4int i)
Definition: G4QMDSystem.hh:83
G4QMDSystem::SubtractSystem
void SubtractSystem(G4QMDSystem *)
Definition: G4QMDSystem.cc:59
G4QMDSystem::GetTotalNumberOfParticipant
G4int GetTotalNumberOfParticipant()
Definition: G4QMDSystem.hh:78
G4QMDSystem::SetParticipant
void SetParticipant(G4QMDParticipant *particle)
Definition: G4QMDSystem.hh:52