ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QMDMeanField.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4QMDMeanField.cc
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 
43 G4QMDMeanField::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 
82 G4QMDMeanField::~G4QMDMeanField()
83 {
84  ;
85 }
86 
87 
88 
89 void 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 
135 void 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 
151 void 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 
260 void 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 
392 void 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 
490 G4double 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 
526 G4double 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 
567 G4double 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 
621 G4bool 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 
637 void 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 
699 std::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 nclst = 1;
759  G4int ichek = 1;
760 
761 
762  G4int id = 0;
763  G4int cluster_id = -1;
764  for ( G4int i = 0 ; i < n-1 ; i++ )
765  {
766 
767  G4bool hasThisCompany = false;
768 // Check only for bryons?
769 // std::cout << "Check Baryon " << i << std::endl;
770 
771  if ( system->GetParticipant( i )->GetBaryonNumber() == 1 )
772  {
773 
774 // if ( is_already_belong_some_cluster[i] != true )
775 // {
776  //G4int j1 = ichek + 1;
777  G4int j1 = i + 1;
778  for ( G4int j = j1 ; j < n ; j++ )
779  {
780 
781  std::vector < G4int > cluster_participants;
782  if ( system->GetParticipant( j )->GetBaryonNumber() == 1 )
783  {
784  G4double rdist2 = rr2[ i ][ j ];
785  G4double pdist2 = pp2[ i ][ j ];
786  //G4double rdist2 = rr2[ num[i] ][ num[j] ];
787  //G4double pdist2 = pp2[ num[i] ][ num[j] ];
788  G4double pcc2 = cpf2
789  * ( rhoa[ i ] + rhoa[ j ] )
790  * ( rhoa[ i ] + rhoa[ j ] );
791 
792 // Check phase space: close enough?
793  if ( rdist2 < rcc2 && pdist2 < pcc2 )
794  {
795 
796 /*
797  G4cout << "G4QMDRESULT "
798  << i << " " << j << " " << id << " "
799  << is_assigned_to [ i ] << " " << is_assigned_to [ j ]
800  << G4endl;
801 */
802 
803  if ( is_assigned_to [ j ] == -1 )
804  {
805  if ( is_assigned_to [ i ] == -1 )
806  {
807  if ( clusters.size() != 0 )
808  {
809  id = clusters.rbegin()->first + 1;
810  //std::cout << "id is increare " << id << std::endl;
811  }
812  else
813  {
814  id = 0;
815  }
816  clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , i ) );
817  is_assigned_to [ i ] = id;
818  clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , j ) );
819  is_assigned_to [ j ] = id;
820  }
821  else
822  {
823  clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ i ] , j ) );
824  is_assigned_to [ j ] = is_assigned_to [ i ];
825  }
826  }
827  else
828  {
829 // j is already belong to some cluester
830  if ( is_assigned_to [ i ] == -1 )
831  {
832  clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , i ) );
833  is_assigned_to [ i ] = is_assigned_to [ j ];
834  }
835  else
836  {
837 // i has companion
838  if ( is_assigned_to [ i ] != is_assigned_to [ j ] )
839  {
840 // move companions to the cluster
841 //
842  //std::cout << "combine " << is_assigned_to [ i ] << " to " << is_assigned_to [ j ] << std::endl;
843  std::multimap< G4int , G4int > clusters_tmp;
844  G4int target_cluster_id;
845  if ( is_assigned_to [ i ] > is_assigned_to [ j ] )
846  target_cluster_id = is_assigned_to [ i ];
847  else
848  target_cluster_id = is_assigned_to [ j ];
849 
850  for ( std::multimap< G4int , G4int >::iterator it
851  = clusters.begin() ; it != clusters.end() ; it++ )
852  {
853 
854  //std::cout << it->first << " " << it->second << " " << target_cluster_id << std::endl;
855  if ( it->first == target_cluster_id )
856  {
857  //std::cout << "move " << it->first << " " << it->second << std::endl;
858  is_assigned_to [ it->second ] = is_assigned_to [ j ];
859  clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , it->second ) );
860  }
861  else
862  {
863  clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( it->first , it->second ) );
864  }
865  }
866 
867  clusters = clusters_tmp;
868  //id = clusters.rbegin()->first;
869  //id = target_cluster_id;
870  //std::cout << "id " << id << std::endl;
871  }
872  }
873  }
874 
875  //std::cout << "combination " << i << " " << j << std::endl;
876  comb_map.insert( std::multimap<G4int,G4int>::value_type ( i , j ) );
877  cluster_participants.push_back ( j );
878 
879 
880 
881  if ( assign_map.find( cluster_id ) == assign_map.end() )
882  {
883  is_already_belong_some_cluster[i] = true;
884  assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , i ) );
885  hasThisCompany = true;
886  }
887  assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , j ) );
888  is_already_belong_some_cluster[j] = true;
889 
890  }
891 
892  if ( ichek == i )
893  {
894  nclst++;
895  ichek++;
896  }
897  }
898 
899  if ( cluster_participants.size() > 0 )
900  {
901 // cluster , participant
902  cluster_map.insert ( std::pair < G4int , std::vector < G4int > > ( i , cluster_participants ) );
903  }
904  }
905 // }
906  }
907  if ( hasThisCompany == true ) cluster_id++;
908  }
909 
910  //std::cout << " id " << id << std::endl;
911 
912 // sort
913 // Heavy cluster comes first
914 // size cluster_id
915  std::multimap< G4int , G4int > sorted_cluster_map;
916  for ( G4int i = 0 ; i <= id ; i++ ) // << "<=" because id is highest cluster nubmer.
917  {
918 
919  //std::cout << i << " cluster has " << clusters.count( i ) << " nucleons." << std::endl;
920  sorted_cluster_map.insert ( std::multimap<G4int,G4int>::value_type ( (G4int) clusters.count( i ) , i ) );
921 
922  }
923 
924 
925 // create nucleus from devided clusters
926  std::vector < G4QMDNucleus* > result;
927  for ( std::multimap < G4int , G4int >::reverse_iterator it
928  = sorted_cluster_map.rbegin() ; it != sorted_cluster_map.rend() ; it ++)
929  {
930 
931  //G4cout << "Add Participants to cluseter " << it->second << G4endl;
932 
933  if ( it->first != 0 )
934  {
935  G4QMDNucleus* nucleus = new G4QMDNucleus();
936  for ( std::multimap < G4int , G4int >::iterator itt
937  = clusters.begin() ; itt != clusters.end() ; itt ++)
938  {
939 
940  if ( it->second == itt->first )
941  {
942  nucleus->SetParticipant( system->GetParticipant ( itt->second ) );
943  //G4cout << "Add Participants " << itt->second << " " << system->GetParticipant ( itt->second )->GetPosition() << G4endl;
944  }
945 
946  }
947  result.push_back( nucleus );
948  }
949 
950  }
951 
952 // delete participants from current system
953 
954  for ( std::vector < G4QMDNucleus* > ::iterator it
955  = result.begin() ; it != result.end() ; it++ )
956  {
957  system->SubtractSystem ( *it );
958  }
959 
960  return result;
961 
962 }
963 
964 
965 
966 void G4QMDMeanField::Update()
967 {
968  SetSystem( system );
969 }