Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/qmd/src/G4LightIonQMDMeanField.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // 230307 Skyrme-QMD parameters added by Y-H. Sato and A. Haga
 29 // 230307 "CalDensityProfile" and "CalChargeDensityProfile" functions added by Y-H. Sato and A. Haga
 30 // 230307 "GetSingleEnergy" and "GetTotalEnergy" functions added by Y-H. Sato and A. Haga
 31 
 32 #include <map>
 33 #include <algorithm>
 34 #include <numeric>
 35 
 36 #include <cmath>
 37 #include <CLHEP/Random/Stat.h>
 38 
 39 #include "G4LightIonQMDMeanField.hh"
 40 #include "G4LightIonQMDParameters.hh"
 41 #include "G4Exp.hh"
 42 #include "G4Pow.hh"
 43 #include "G4PhysicalConstants.hh"
 44 #include "Randomize.hh"
 45 
 46 G4LightIonQMDMeanField::G4LightIonQMDMeanField()
 47 {
 48    G4LightIonQMDParameters* parameters = G4LightIonQMDParameters::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    eta = parameters->Get_eta(); // Skyrme-QMD
 55    kappas = parameters->Get_kappas(); // Skyrme-QMD
 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    g0 = parameters->Get_g0(); // Skyrme-QMD
 65    g0iso = parameters->Get_g0iso(); // Skyrme-QMD
 66    gtau0 = parameters->Get_gtau0(); // Skyrme-QMD
 67 
 68    // distance
 69    c0w = 1.0/4.0/wl;
 70    c0sw = std::sqrt( c0w );
 71    clw = 2.0 / std::sqrt ( 4.0 * pi * wl );
 72 
 73    // graduate
 74    c0g = - c0 / ( 2.0 * wl );
 75    c3g = - c3 / ( 4.0 * wl ) * gamm;
 76    csg = - cs / ( 2.0 * wl );
 77    pag = gamm - 1;
 78    pag_tau = eta - 1; // Skyrme-QMD
 79    cg0 = - g0 / ( 2.0 * wl );  // Skyrme-QMD
 80    cgtau0 = - gtau0 / ( 4.0 * wl ) * eta;  // Skyrme-QMD
 81 
 82    system = nullptr; // will be set through SetSystem method
 83 }
 84 
 85 void G4LightIonQMDMeanField::SetSystem ( G4QMDSystem* aSystem )
 86 { 
 87    system = aSystem;
 88 
 89    G4int n = system->GetTotalNumberOfParticipant();
 90   
 91    pp2.clear();
 92    rr2.clear();
 93    rbij.clear();
 94    rha.clear();
 95    rhe.clear();
 96    rhc.clear();
 97 
 98    rr2.resize( n );
 99    pp2.resize( n );
100    rbij.resize( n );
101    rha.resize( n );
102    rhe.resize( n );
103    rhc.resize( n );
104 
105    for ( G4int i = 0 ; i < n ; ++i )
106    {
107       rr2[i].resize( n );
108       pp2[i].resize( n );
109       rbij[i].resize( n );
110       rha[i].resize( n );
111       rhe[i].resize( n );
112       rhc[i].resize( n );
113    }
114 
115    ffr.clear();
116    ffp.clear();
117    rh3d.clear();
118    rh3d_tau.clear(); // Skyrme-QMD
119 
120    ffr.resize( n );
121    ffp.resize( n );
122    rh3d.resize( n );
123    rh3d_tau.resize( n ); // Skyrme-QMD
124 
125    Cal2BodyQuantities();
126 }
127 
128 void G4LightIonQMDMeanField::SetNucleus ( G4LightIonQMDNucleus* aNucleus ) 
129 {
130    SetSystem( aNucleus );
131 
132    G4double totalPotential = GetTotalPotential(); 
133    aNucleus->SetTotalPotential( totalPotential );
134    aNucleus->CalEnergyAndAngularMomentumInCM();
135 }
136 
137 void G4LightIonQMDMeanField::Cal2BodyQuantities()
138 {
139    if ( system->GetTotalNumberOfParticipant() < 2 )  { return; }
140 
141    for ( G4int j = 1 ; j < system->GetTotalNumberOfParticipant() ; ++j )
142    {
143       G4ThreeVector rj = system->GetParticipant( j )->GetPosition();
144       G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
145 
146       for ( G4int i = 0 ; i < j ; ++i )
147       {
148          G4ThreeVector ri = system->GetParticipant( i )->GetPosition();
149          G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum();
150 
151          G4ThreeVector rij = ri - rj;
152          G4ThreeVector pij = (p4i - p4j).v();
153          G4LorentzVector p4ij = p4i - p4j;
154          G4ThreeVector bij = ( p4i + p4j ).boostVector();
155          G4double gammaij = ( p4i + p4j ).gamma();
156 
157          G4double eij = ( p4i + p4j ).e();
158 
159          G4double rbrb = rij*bij;
160          G4double rij2 = rij*rij;
161          G4double pij2 = pij*pij;
162 
163          rbrb = irelcr * rbrb;
164          G4double  gamma2_ij = gammaij*gammaij;
165 
166          rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb;
167          rr2[j][i] = rr2[i][j];
168 
169          rbij[i][j] = gamma2_ij * rbrb;
170          rbij[j][i] = - rbij[i][j];
171 
172          pp2[i][j] = pij2
173                    + irelcr * ( - G4Pow::GetInstance()->powN ( p4i.e() - p4j.e() , 2 )
174                    + gamma2_ij * G4Pow::GetInstance()->powN ( ( ( p4i.m2() - p4j.m2() ) / eij ) , 2 ) );
175 
176 
177          pp2[j][i] = pp2[i][j];
178 
179          // Gauss term
180 
181          G4double expa1 = - rr2[i][j] * c0w;
182 
183          G4double rh1;
184          if ( expa1 > epsx )
185          {
186             rh1 = G4Exp( expa1 );
187          }
188          else
189          {
190             rh1 = 0.0;
191          }
192 
193          G4int ibry = system->GetParticipant(i)->GetBaryonNumber();
194          G4int jbry = system->GetParticipant(j)->GetBaryonNumber();
195 
196          rha[i][j] = ibry*jbry*rh1;
197          rha[j][i] = rha[i][j];
198 
199          // Coulomb terms
200 
201          G4double rrs2 = rr2[i][j] + epscl;
202          G4double rrs = std::sqrt ( rrs2 );
203 
204          G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
205          G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
206 
207          G4double xerf = 0.0;
208          // T. K. add this protection. 5.8 is good enough for double
209          if ( rrs*c0sw < 5.8 )
210          {
211 #if defined WIN32-VC
212             xerf = CLHEP::HepStat::erf ( rrs*c0sw );
213 #else
214             xerf = std::erf ( rrs*c0sw );
215 #endif
216          } 
217          else
218          {
219             xerf = 1.0;
220          }
221 
222          G4double erfij = xerf/rrs;
223 
224          rhe[i][j] = icharge*jcharge * erfij;
225          rhe[j][i] = rhe[i][j];
226          rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2;
227          rhc[j][i] = rhc[i][j];
228       }  // i
229    }  // j
230 }
231 
232 void G4LightIonQMDMeanField::Cal2BodyQuantities( G4int i )
233 {
234    G4ThreeVector ri = system->GetParticipant( i )->GetPosition();  
235    G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum();  
236 
237    for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j )
238    {
239       if ( j == i )  { continue; }
240 
241       G4ThreeVector rj = system->GetParticipant( j )->GetPosition();  
242       G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();  
243 
244       G4ThreeVector rij = ri - rj;
245       G4ThreeVector pij = (p4i - p4j).v();
246       G4LorentzVector p4ij = p4i - p4j;
247       G4ThreeVector bij = ( p4i + p4j ).boostVector();
248       G4double gammaij = ( p4i + p4j ).gamma();
249 
250       G4double eij = ( p4i + p4j ).e();
251 
252       G4double rbrb = rij*bij;
253       G4double rij2 = rij*rij;
254       G4double pij2 = pij*pij;
255 
256       rbrb = irelcr * rbrb;
257       G4double  gamma2_ij = gammaij*gammaij;
258 
259       rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb;
260       rr2[j][i] = rr2[i][j];
261 
262       rbij[i][j] = gamma2_ij * rbrb;
263       rbij[j][i] = - rbij[i][j];
264       
265       pp2[i][j] = pij2
266                 + irelcr * ( - G4Pow::GetInstance()->powN ( p4i.e() - p4j.e() , 2 )
267                 + gamma2_ij * G4Pow::GetInstance()->powN ( ( ( p4i.m2() - p4j.m2() ) / eij ) , 2 ) );
268 
269       pp2[j][i] = pp2[i][j];
270 
271       // Gauss term
272 
273       G4double expa1 = - rr2[i][j] * c0w;  
274 
275       G4double rh1;
276       if ( expa1 > epsx ) 
277       {
278          rh1 = G4Exp( expa1 );
279       }
280       else 
281       {
282          rh1 = 0.0;  
283       } 
284 
285       G4int ibry = system->GetParticipant(i)->GetBaryonNumber();  
286       G4int jbry = system->GetParticipant(j)->GetBaryonNumber();  
287 
288       rha[i][j] = ibry*jbry*rh1;  
289       rha[j][i] = rha[i][j]; 
290 
291       // Coulomb terms
292 
293       G4double rrs2 = rr2[i][j] + epscl; 
294       G4double rrs = std::sqrt ( rrs2 ); 
295 
296       G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
297       G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
298 
299       G4double xerf = 0.0;
300       // T. K. add this protection. 5.8 is good enough for double
301       if ( rrs*c0sw < 5.8 ) 
302       {
303 #if defined WIN32-VC
304          xerf = CLHEP::HepStat::erf ( rrs*c0sw );
305 #else
306          xerf = std::erf ( rrs*c0sw );
307 #endif
308       } 
309       else 
310       {
311          xerf = 1.0;
312       }
313 
314       G4double erfij = xerf/rrs; 
315       
316       rhe[i][j] = icharge*jcharge * erfij;
317       rhe[j][i] = rhe[i][j];
318       rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2;
319       rhc[j][i] = rhc[i][j];
320    }
321 }
322 
323 void G4LightIonQMDMeanField::CalGraduate()
324 {
325    ffr.resize( system->GetTotalNumberOfParticipant() );
326    ffp.resize( system->GetTotalNumberOfParticipant() );
327    rh3d.resize( system->GetTotalNumberOfParticipant() );
328    rh3d_tau.resize( system->GetTotalNumberOfParticipant() ); // Skyrme-QMD
329 
330    for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; ++i )
331    {
332       G4double rho3 = 0.0;
333       for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j )
334       {
335          rho3 += rha[j][i]; 
336       }
337       rh3d[i] = G4Pow::GetInstance()->powA ( rho3 , pag );
338       rh3d_tau[i] = G4Pow::GetInstance()->powA ( rho3 , pag_tau ); // Skyrme-QMD
339    }
340 
341    for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; ++i )
342    {
343       G4ThreeVector ri = system->GetParticipant( i )->GetPosition();  
344       G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum();  
345 
346       G4ThreeVector betai = p4i.v()/p4i.e();
347       
348       // R-JQMD
349       G4double Vi = GetPotential( i );
350       G4double p_zero = std::sqrt( p4i.e()*p4i.e() + 2*p4i.m()*Vi);
351       G4ThreeVector betai_R = p4i.v()/p_zero;
352       G4double mi_R = p4i.m()/p_zero;
353 
354       ffr[i] = betai_R;
355       ffp[i] = G4ThreeVector( 0.0 );
356 
357       for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j )
358       {
359          G4ThreeVector rj = system->GetParticipant( j )->GetPosition();  
360          G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();  
361 
362          G4double eij = p4i.e() + p4j.e(); 
363 
364          G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
365          G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
366 
367          G4int inuc = system->GetParticipant(i)->GetNuc();
368          G4int jnuc = system->GetParticipant(j)->GetNuc();
369 
370          G4double fsij = 3.0/(2*wl) - rr2[j][i]/(2*wl)/(2*wl); // Add for Skyrme-QMD
371 
372          G4double ccpp = c0g * rha[j][i]
373                        + c3g * rha[j][i] * ( rh3d[j] + rh3d[i] )
374                        + cg0 * rha[j][i]/wl
375                        + cg0 * rha[j][i] * fsij
376                        + cgtau0 * rha[j][i] * ( rh3d_tau[j] + rh3d_tau[i] )
377                        + csg * rha[j][i] * jnuc * inuc
378                              * ( 1. - 2. * std::abs( jcharge - icharge ) )
379                              * (1. - kappas * fsij + kappas / wl)
380                        + cl * rhc[j][i];
381 
382          ccpp *= mi_R;
383 
384          G4double grbb = - rbij[j][i];
385          G4double ccrr = grbb * ccpp / eij;
386 
387          G4ThreeVector rij = ri - rj;   
388          G4ThreeVector betaij =  ( p4i + p4j ).v()/eij;
389          G4ThreeVector cij = betaij - betai;   
390 
391          ffr[i] = ffr[i] + 2*ccrr* ( rij + grbb*cij );
392          ffp[i] = ffp[i] - 2*ccpp* ( rij + grbb*betaij );
393       }
394    }
395 }
396 
397 G4double G4LightIonQMDMeanField::GetPotential( G4int i )
398 {
399    G4int n = system->GetTotalNumberOfParticipant();
400 
401    G4double rhoa = 0.0;
402    G4double rho3 = 0.0;
403    G4double fsij_rhoa = 0.0; // Skyrme-QMD
404    //G4double fsij_rhos = 0.0; // Skyrme-QMD
405    G4double rho3_tau = 0.0; // Skyrme-QMD
406    G4double rhos = 0.0;
407    G4double rhoc = 0.0;
408 
409    G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
410    G4int inuc = system->GetParticipant(i)->GetNuc();
411 
412    for ( G4int j = 0 ; j < n ; ++j )
413    {
414       G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
415       G4int jnuc = system->GetParticipant(j)->GetNuc();
416       G4double fsij = 3.0/(2*wl) - rr2[j][i]/(2*wl)/(2*wl); // Add for Skyrme-QMD
417 
418       rhoa += rha[j][i];
419       fsij_rhoa += fsij * rha[j][i]; // Skyrme-QMD
420       rhoc += rhe[j][i];
421       rhos += rha[j][i] * jnuc * inuc
422                         * ( 1. - 2. * std::abs( jcharge - icharge ) )  // Skyrme-QMD
423                         * (1. - kappas * fsij);  // Skyrme-QMD
424    }
425 
426    rho3 = G4Pow::GetInstance()->powA ( rhoa , gamm );
427    rho3_tau = G4Pow::GetInstance()->powA ( rhoa , eta );
428 
429    G4double potential = c0 * rhoa
430                       + c3 * rho3
431                       + g0 * fsij_rhoa  // Skyrme-QMD
432                       //+ g0iso * fsij_rhos  // Skyrme-QMD
433                       + gtau0 * rho3_tau  // Skyrme-QMD
434                       + cs * rhos
435                       + cl * rhoc;
436    return potential;
437 }
438 
439 G4double G4LightIonQMDMeanField::GetTotalPotential()
440 {
441    G4int n = system->GetTotalNumberOfParticipant();
442 
443    std::vector < G4double > rhoa ( n , 0.0 ); 
444    std::vector < G4double > rho3 ( n , 0.0 );
445    std::vector < G4double > rho3_tau ( n , 0.0 ); // Skyrme-QMD
446    //std::vector < G4double > fsij_rhos ( n , 0.0 ); // Skyrme-QMD
447    std::vector < G4double > fsij_rhoa ( n , 0.0 ); // Skyrme-QMD
448    std::vector < G4double > rhos ( n , 0.0 );
449    std::vector < G4double > rhoc ( n , 0.0 ); 
450 
451    for ( G4int i = 0 ; i < n ; ++i )
452    {
453       G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
454       G4int inuc = system->GetParticipant(i)->GetNuc();
455 
456       for ( G4int j = 0 ; j < n ; ++j )
457       {
458          G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
459          G4int jnuc = system->GetParticipant(j)->GetNuc();
460          G4double fsij = 3.0/(2*wl) - rr2[j][i]/(2*wl)/(2*wl); // Add for Skyrme-QMD
461 
462          rhoa[i] += rha[j][i];
463          fsij_rhoa[i] += fsij * rha[j][i]; // Skyrme-QMD
464          rhoc[i] += rhe[j][i];
465          rhos[i] += rha[j][i] * jnuc * inuc
466                               //* ( 1 - 2 * std::abs ( jcharge - icharge ) );
467                               * ( 1. - 2. * std::abs( jcharge - icharge ) )  // Skyrme-QMD
468                               * (1. - kappas * fsij);  // Skyrme-QMD
469          //fsij_rhos[i] += fsij * rha[j][i] * jnuc * inuc
470                                 //* ( 1. - 2. * std::abs( jcharge - icharge ) )  // Skyrme-QMD
471                                 //* (1. - kappas * fsij);  // Skyrme-QMD
472       }
473 
474       rho3[i] = G4Pow::GetInstance()->powA ( rhoa[i] , gamm );
475       rho3_tau[i] = G4Pow::GetInstance()->powA ( rhoa[i] , eta );
476    }
477 
478     G4double potential = c0 * std::accumulate( rhoa.cbegin() , rhoa.cend() , 0.0 )
479                        + c3 * std::accumulate( rho3.cbegin() , rho3.cend() , 0.0 )
480                        + g0 * std::accumulate( fsij_rhoa.cbegin() , fsij_rhoa.cend() , 0.0 )
481                        //+ g0iso * std::accumulate( fsij_rhos.cbegin() , fsij_rhos.cend() , 0.0 )
482                        + gtau0 * std::accumulate( rho3_tau.cbegin() , rho3_tau.cend() , 0.0 )
483                        + cs * std::accumulate( rhos.cbegin() , rhos.cend() , 0.0 )
484                        + cl * std::accumulate( rhoc.cbegin() , rhoc.cend() , 0.0 );
485 
486    return potential;
487 }
488 
489 G4double G4LightIonQMDMeanField::GetSingleEnergy( G4int j )
490 {
491     G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
492     G4double emass = p4j.m();
493     G4double ekinal2 = p4j.e()*p4j.e();
494     G4double esingle = std::sqrt(ekinal2 + 2*emass*GetPotential(j));
495     return esingle;
496 }
497 
498 G4double G4LightIonQMDMeanField::GetTotalEnergy()
499 {
500     
501     G4int n = system->GetTotalNumberOfParticipant();
502     G4double etotal = 0.0;
503     for ( int j = 0 ; j < n ; j++ )
504     {
505         G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum();
506         G4double emass = p4j.m();
507         G4double ekinal2 = p4j.e()*p4j.e();
508         etotal += std::sqrt(ekinal2 + 2*emass*GetPotential(j));
509     }
510     return etotal;
511 
512 }
513 
514 G4double G4LightIonQMDMeanField::calPauliBlockingFactor( G4int i )
515 {
516    // i is supposed beyond total number of Participant()
517 
518    G4double pf = 0.0;
519    G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
520 
521    for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; ++j )
522    {
523       G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
524       G4int jnuc = system->GetParticipant(j)->GetNuc();
525 
526       if ( jcharge == icharge && jnuc == 1 )
527       {
528          G4double expa = -rr2[i][j]*cpw;
529          if ( expa > epsx ) 
530          {
531             expa = expa - pp2[i][j]*cph;
532             if ( expa > epsx ) 
533             {
534                pf = pf + G4Exp ( expa );
535             }
536          }
537       }
538    }
539 
540    return ( pf - 1.0 ) * cpc;
541 }
542 
543 G4bool G4LightIonQMDMeanField::IsPauliBlocked( G4int i )
544 {
545     G4bool result = false; 
546     
547     if ( system->GetParticipant( i )->GetNuc() == 1 )
548     {
549        G4double pf = calPauliBlockingFactor( i );
550        G4double rand = G4UniformRand(); 
551        if ( pf > rand ) { result = true; }
552     }
553 
554     return result; 
555 }
556 
557 void G4LightIonQMDMeanField::DoPropagation( G4double dt )
558 {
559    G4double cc2 = 1.0; 
560    G4double cc1 = 1.0 - cc2; 
561    G4double cc3 = 1.0 / 2.0 / cc2; 
562 
563    G4double dt3 = dt * cc3;
564    G4double dt1 = dt * ( cc1 - cc3 );
565    G4double dt2 = dt * cc2;
566 
567    CalGraduate(); 
568 
569    G4int n = system->GetTotalNumberOfParticipant();
570 
571    // 1st Step
572 
573    std::vector< G4ThreeVector > f0r, f0p;
574    f0r.resize( n );
575    f0p.resize( n );
576 
577    for ( G4int i = 0 ; i < n ; ++i )
578    {
579       G4ThreeVector ri = system->GetParticipant( i )->GetPosition();  
580       G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum();  
581 
582       ri += dt3* ffr[i];
583       p3i += dt3* ffp[i];
584 
585       f0r[i] = ffr[i];
586       f0p[i] = ffp[i];
587       
588       system->GetParticipant( i )->SetPosition( ri );  
589       system->GetParticipant( i )->SetMomentum( p3i );  
590 
591       // we do not need set total momentum by ourselvs
592    }
593 
594    // 2nd Step
595 
596    Cal2BodyQuantities();
597    CalGraduate(); 
598 
599    for ( G4int i = 0 ; i < n ; ++i )
600    {
601       G4ThreeVector ri = system->GetParticipant( i )->GetPosition();  
602       G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum();  
603 
604       ri += dt1* f0r[i] + dt2* ffr[i];
605       p3i += dt1* f0p[i] + dt2* ffp[i];
606 
607       system->GetParticipant( i )->SetPosition( ri );  
608       system->GetParticipant( i )->SetMomentum( p3i );  
609 
610       // we do not need set total momentum by ourselvs
611    }
612 
613    Cal2BodyQuantities();
614 }
615 
616 std::vector< G4LightIonQMDNucleus* > G4LightIonQMDMeanField::DoClusterJudgment()
617 {
618    Cal2BodyQuantities(); 
619 
620    G4double cpf2 = G4Pow::GetInstance()->A23 ( 1.5 * pi*pi * G4Pow::GetInstance()->powA ( 4.0 * pi * wl , -1.5 ) ) * hbc * hbc;
621    G4double rcc2 = rclds*rclds;
622 
623    G4int n = system->GetTotalNumberOfParticipant();
624    std::vector < G4double > rhoa;
625    rhoa.resize ( n );
626 
627    for ( G4int i = 0 ; i < n ; ++i )
628    {
629      rhoa[i] = 0.0;
630 
631      if ( system->GetParticipant( i )->GetBaryonNumber() == 1 )
632      {
633        for ( G4int j = 0 ; j < n ; ++j )
634        {
635          if ( system->GetParticipant( j )->GetBaryonNumber() == 1 )
636          rhoa[i] += rha[i][j];
637        }
638      }
639 
640      rhoa[i] = G4Pow::GetInstance()->A13 ( rhoa[i] + 1 );
641    }
642 
643    // identification of the cluster
644    std::vector < G4bool > is_already_belong_some_cluster;
645 
646    //         cluster_id   participant_id
647    std::multimap < G4int , G4int > comb_map; 
648    std::multimap < G4int , G4int > assign_map; 
649    assign_map.clear(); 
650 
651    std::vector < G4int > mascl;
652    std::vector < G4int > num;
653    mascl.resize ( n );
654    num.resize ( n );
655    is_already_belong_some_cluster.resize ( n );
656 
657    std::vector < G4int > is_assigned_to ( n , -1 );
658    std::multimap < G4int , G4int > clusters;
659 
660    for ( G4int i = 0 ; i < n ; ++i )
661    {
662      mascl[i] = 1;
663      num[i] = 1;
664      is_already_belong_some_cluster[i] = false;
665    }
666 
667    G4int ichek = 1;
668    G4int id = 0;
669    G4int cluster_id = -1;  
670    for ( G4int i = 0 ; i < n-1 ; ++i )
671    {
672      G4bool hasThisCompany = false;
673 
674      if ( system->GetParticipant( i )->GetBaryonNumber() == 1 )
675      {
676        G4int j1 = i + 1;
677        for ( G4int j = j1 ; j < n ; ++j )
678        {
679          std::vector < G4int > cluster_participants;
680          if ( system->GetParticipant( j )->GetBaryonNumber() == 1 )  
681          {
682            G4double rdist2 = rr2[ i ][ j ];
683            G4double pdist2 = pp2[ i ][ j ];
684            G4double pcc2 = cpf2
685                          * ( rhoa[ i ] + rhoa[ j ] )
686                          * ( rhoa[ i ] + rhoa[ j ] );
687 
688            // Check phase space: close enough?
689            if ( rdist2 < rcc2 && pdist2 < pcc2 )
690            {
691              if ( is_assigned_to [ j ] == -1 )
692              {
693                if ( is_assigned_to [ i ] == -1 )
694                {
695                  if ( clusters.size() != 0 )
696                  {
697                    id = clusters.rbegin()->first + 1;
698                  }
699                  else
700                  {
701                    id = 0;
702                  }
703                  clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , i ) );
704                  is_assigned_to [ i ] = id;
705                  clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , j ) );
706                  is_assigned_to [ j ] = id;
707                }
708                else
709                {
710                  clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ i ] , j ) );
711                  is_assigned_to [ j ] = is_assigned_to [ i ];
712                }
713              }
714              else
715              {
716                // j is already belong to some cluster
717                if ( is_assigned_to [ i ] == -1 )
718                {
719                  clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , i ) );
720                  is_assigned_to [ i ] = is_assigned_to [ j ];
721                }
722                else
723                {
724                  // i has companion
725                  if ( is_assigned_to [ i ] != is_assigned_to [ j ] )
726                  {
727                    // move companions to the cluster
728                    std::multimap< G4int , G4int > clusters_tmp;
729                    G4int target_cluster_id;
730                    if ( is_assigned_to [ i ] > is_assigned_to [ j ] )
731                    {
732                      target_cluster_id = is_assigned_to [ i ];
733                    }
734                    else
735                    {
736                      target_cluster_id = is_assigned_to [ j ];
737                    }
738                    for ( auto it = clusters.cbegin() ; it != clusters.cend() ; ++it )
739                    {
740                      if ( it->first == target_cluster_id )
741                      {
742                        is_assigned_to [ it->second ] = is_assigned_to [ j ];
743                        clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type (  is_assigned_to [ j ] , it->second ) );
744                      }
745                      else
746                      {
747                        clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( it->first , it->second ) );
748                      }
749                    }
750                    clusters = std::move(clusters_tmp);
751                  }
752                }
753              }
754 
755              comb_map.insert( std::multimap<G4int,G4int>::value_type ( i , j ) );
756              cluster_participants.push_back ( j );
757 
758              if ( assign_map.find( cluster_id ) == assign_map.end() )
759              {
760                is_already_belong_some_cluster[i] = true;
761                assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , i ) );
762                hasThisCompany = true;
763              }
764              assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , j ) );
765              is_already_belong_some_cluster[j] = true;
766            }
767 
768            if ( ichek == i )
769            {
770              ++ichek;
771            }
772          }
773        }
774      }
775      if ( hasThisCompany == true ) { ++cluster_id; }
776    }
777 
778    // sort
779    // Heavy cluster comes first
780    //             size    cluster_id
781    std::multimap< G4int , G4int > sorted_cluster_map;
782    for ( G4int i = 0 ; i <= id ; ++i )  // << "<=" because id is highest cluster nubmer.
783    {
784      sorted_cluster_map.insert ( std::multimap<G4int,G4int>::value_type ( (G4int) clusters.count( i ) , i ) );
785    }
786 
787    // create nucleus from divided clusters
788    std::vector < G4LightIonQMDNucleus* > result;
789    for ( auto it = sorted_cluster_map.crbegin(); it != sorted_cluster_map.crend(); ++it )
790    {
791      if ( it->first != 0 )
792      {
793        G4LightIonQMDNucleus* nucleus = new G4LightIonQMDNucleus();
794        for ( auto itt = clusters.cbegin(); itt != clusters.cend(); ++itt )
795        {
796          if ( it->second == itt->first )
797          {
798            nucleus->SetParticipant( system->GetParticipant ( itt->second ) );
799          }
800        }
801        result.push_back( nucleus );
802      }
803    }
804 
805    // delete participants from current system
806    for ( auto it = result.cbegin(); it != result.cend(); ++it )
807    {
808      system->SubtractSystem ( *it );
809    }
810    
811    return result;
812 }
813 
814 void G4LightIonQMDMeanField::Update() 
815 { 
816    SetSystem( system ); 
817 }
818