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