Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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 <CLHEP/Random/Stat.h>
34 
35 #include "G4QMDMeanField.hh"
36 #include "G4QMDParameters.hh"
37 
38 #include "G4PhysicalConstants.hh"
39 #include "Randomize.hh"
40 
41 G4QMDMeanField::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 
79 G4QMDMeanField::~G4QMDMeanField()
80 {
81  ;
82 }
83 
84 
85 
86 void 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 
132 void 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 
148 void 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 
248 void 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 
373 void 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 
471 G4double 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 
507 G4double 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 
548 G4double 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 
602 G4bool 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 
618 void 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 
680 std::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 
950 void G4QMDMeanField::Update()
951 {
952  SetSystem( system );
953 }
Generated on Sat May 25 2013 14:34:06 for Geant4 by   doxygen 1.8.4