69 clw = 2.0 / std::sqrt ( 4.0 *
pi *
wl );
112 for (
int i = 0 ; i <
n ; i++ )
162 for (
G4int i = 0 ; i < j ; i++ )
172 G4double gammaij = ( p4i + p4j ).gamma();
182 G4double gamma2_ij = gammaij*gammaij;
185 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb;
188 rbij[i][j] = gamma2_ij * rbrb;
205 rh1 =
G4Exp( expa1 );
216 rha[i][j] = ibry*jbry*rh1;
229 if ( rrs*
c0sw < 5.8 ) {
237 xerf = erf ( rrs*c0sw );
246 rhe[i][j] = icharge*jcharge * erfij;
250 rhc[i][j] = icharge*jcharge * ( - erfij +
clw * rh1 ) / rrs2;
252 rhc[j][i] = rhc[i][j];
271 if ( j == i )
continue;
280 G4double gammaij = ( p4i + p4j ).gamma();
290 G4double gamma2_ij = gammaij*gammaij;
318 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb;
321 rbij[i][j] = gamma2_ij * rbrb;
337 rh1 =
G4Exp( expa1 );
348 rha[i][j] = ibry*jbry*rh1;
361 if ( rrs*
c0sw < 5.8 ) {
367 xerf = erf ( rrs*c0sw );
376 rhe[i][j] = icharge*jcharge * erfij;
382 rhc[i][j] = icharge*jcharge * ( - erfij +
clw * rh1 ) / rrs2;
384 rhc[j][i] = rhc[i][j];
420 G4double p_zero = std::sqrt( p4i.
e()*p4i.
e() + 2*p4i.
m()*Vi);
449 +
csg *
rha[j][i] * jnuc * inuc
450 * ( 1. - 2. * std::abs( jcharge - icharge ) )
476 ffr[i] =
ffr[i] + 2*ccrr* ( rij + grbb*cij );
478 ffp[i] = ffp[i] - 2*ccpp* ( rij + grbb*betaij );
503 for (
G4int j = 0 ; j <
n ; j ++ )
510 rhos +=
rha[j][i] * jnuc * inuc
511 * ( 1 - 2 * std::abs ( jcharge - icharge ) );
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 );
537 for (
G4int i = 0 ; i <
n ; i ++ )
542 for (
G4int j = 0 ; j <
n ; j ++ )
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 ) );
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 );
580 if ( jcharge == icharge && jnuc == 1 )
593 expa = expa -
pp2[i][j]*
cph;
603 pf = pf +
G4Exp ( expa );
611 pf = ( pf - 1.0 ) *
cpc;
629 if ( pf > rand ) result =
true;
654 std::vector< G4ThreeVector > f0r, f0p;
658 for (
G4int i = 0 ; i <
n ; i++ )
679 for (
G4int i = 0 ; i <
n ; i++ )
684 ri += dt1* f0r[i] + dt2*
ffr[i];
685 p3i += dt1* f0p[i] + dt2*
ffp[i];
710 std::vector < G4double > rhoa;
713 for (
G4int i = 0 ; i <
n ; i++ )
719 for (
G4int j = 0 ; j <
n ; j++ )
722 rhoa[i] +=
rha[i][j];
732 std::map < G4int , std::vector < G4int > > cluster_map;
733 std::vector < G4bool > is_already_belong_some_cluster;
736 std::multimap < G4int , G4int > comb_map;
737 std::multimap < G4int , G4int > assign_map;
740 std::vector < G4int > mascl;
741 std::vector < G4int > num;
744 is_already_belong_some_cluster.resize ( n );
746 std::vector < G4int > is_assigned_to ( n , -1 );
747 std::multimap < G4int , G4int > clusters;
749 for (
G4int i = 0 ; i <
n ; i++ )
754 is_already_belong_some_cluster[i] =
false;
763 G4int cluster_id = -1;
764 for (
G4int i = 0 ; i < n-1 ; i++ )
767 G4bool hasThisCompany =
false;
778 for (
G4int j = j1 ; j <
n ; j++ )
781 std::vector < G4int > cluster_participants;
789 * ( rhoa[ i ] + rhoa[ j ] )
790 * ( rhoa[ i ] + rhoa[ j ] );
793 if ( rdist2 < rcc2 && pdist2 < pcc2 )
803 if ( is_assigned_to [ j ] == -1 )
805 if ( is_assigned_to [ i ] == -1 )
807 if ( clusters.size() != 0 )
809 id = clusters.rbegin()->first + 1;
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;
823 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ i ] , j ) );
824 is_assigned_to [ j ] = is_assigned_to [ i ];
830 if ( is_assigned_to [ i ] == -1 )
832 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , i ) );
833 is_assigned_to [ i ] = is_assigned_to [ j ];
838 if ( is_assigned_to [ i ] != is_assigned_to [ j ] )
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 ];
848 target_cluster_id = is_assigned_to [ j ];
850 for ( std::multimap< G4int , G4int >::iterator it
851 = clusters.begin() ; it != clusters.end() ; it++ )
855 if ( it->first == target_cluster_id )
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 ) );
863 clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( it->first , it->second ) );
867 clusters = clusters_tmp;
876 comb_map.insert( std::multimap<G4int,G4int>::value_type ( i , j ) );
877 cluster_participants.push_back ( j );
881 if ( assign_map.find( cluster_id ) == assign_map.end() )
883 is_already_belong_some_cluster[i] =
true;
884 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , i ) );
885 hasThisCompany =
true;
887 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , j ) );
888 is_already_belong_some_cluster[j] =
true;
899 if ( cluster_participants.size() > 0 )
902 cluster_map.insert ( std::pair <
G4int , std::vector < G4int > > ( i , cluster_participants ) );
907 if ( hasThisCompany ==
true ) cluster_id++;
915 std::multimap< G4int , G4int > sorted_cluster_map;
916 for (
G4int i = 0 ; i <= id ; i++ )
920 sorted_cluster_map.insert ( std::multimap<G4int,G4int>::value_type ( (
G4int) clusters.count( i ) , i ) );
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 ++)
933 if ( it->first != 0 )
936 for ( std::multimap < G4int , G4int >::iterator itt
937 = clusters.begin() ; itt != clusters.end() ; itt ++)
940 if ( it->second == itt->first )
947 result.push_back( nucleus );
954 for ( std::vector < G4QMDNucleus* > ::iterator it
955 = result.begin() ; it != result.end() ; it++ )
G4int GetChargeInUnitOfEplus()
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
CLHEP::Hep3Vector G4ThreeVector
G4double A13(G4double A) const
std::vector< std::vector< G4double > > rha
G4double GetTotalPotential()
void CalEnergyAndAngularMomentumInCM()
G4int GetTotalNumberOfParticipant()
void SubtractSystem(G4QMDSystem *)
G4LorentzVector Get4Momentum()
G4double A23(G4double A) const
void SetMomentum(G4ThreeVector p)
void SetTotalPotential(G4double x)
std::vector< G4ThreeVector > ffp
G4double powN(G4double x, G4int n) const
std::vector< std::vector< G4double > > rr2
static G4Pow * GetInstance()
std::vector< G4ThreeVector > ffr
G4double powA(G4double A, G4double y) const
std::vector< std::vector< G4double > > rbij
void SetPosition(G4ThreeVector r)
void DoPropagation(G4double)
std::vector< std::vector< G4double > > rhc
std::vector< std::vector< G4double > > rhe
G4double G4ParticleHPJENDLHEData::G4double result
G4double GetPotential(G4int)
void SetParticipant(G4QMDParticipant *particle)
G4ThreeVector GetMomentum()
std::vector< G4double > rh3d
void SetSystem(G4QMDSystem *aSystem)
std::vector< G4QMDNucleus * > DoClusterJudgment()
G4double calPauliBlockingFactor(G4int)
static double erf(double x)
static constexpr double pi
G4QMDParticipant * GetParticipant(G4int i)
void Cal2BodyQuantities()
static G4QMDParameters * GetInstance()
G4bool IsPauliBlocked(G4int)
G4ThreeVector GetPosition()
void SetNucleus(G4QMDNucleus *aSystem)
std::vector< std::vector< G4double > > pp2