Changeset 1347 for trunk/source/processes/hadronic/models/qmd

Timestamp:

Dec 22, 2010, 3:52:27 PM (15 years ago)

Author:

garnier

Message:

geant4 tag 9.4

Location:

trunk/source/processes/hadronic/models/qmd

Files:

• History (modified) (1 diff)
• include/G4QMDMeanField.hh (modified) (1 diff)
• include/G4QMDReaction.hh (modified) (2 diffs)
• src/G4QMDMeanField.cc (modified) (4 diffs)
• src/G4QMDReaction.cc (modified) (2 diffs)

trunk/source/processes/hadronic/models/qmd/History

@@ -13 +13 @@
      * Please list in reverse chronological order (last date on top)
      ---------------------------------------------------------------
+
+18-November-2010 Tatsumi Koi (hadr-qmd-V09-03-00)
+- Enable RQMD
+        include/G4QMDReaction.hh  
+        src/G4QMDReaction.cc  
+        include/G4QMDMeanField.hh
+        src/G4QMDMeanField.cc
 
 20-November-2009 Tatsumi Koi (hadr-qmd-V09-02-03)

trunk/source/processes/hadronic/models/qmd/include/G4QMDMeanField.hh

@@ -60 +60 @@
 
       G4double GetTotalPotential(); 
+      G4double GetPotential( G4int );
 
       void DoPropagation( G4double );

trunk/source/processes/hadronic/models/qmd/include/G4QMDReaction.hh

@@ -74 +74 @@
       void EnableFermiG(){ heg = true; setHighEnergyModel(); };
 
+      void SetTMAX( G4int i ){ maxTime = i; };
+      void SetDT( G4double t ){ deltaT = t; };
+      void SetEF( G4double x ){ envelopF = x; };
+
    private:
 
@@ -89 +93 @@
       G4double deltaT;
       G4int maxTime;
+      G4double envelopF;
 
       G4Evaporation* evaporation;

trunk/source/processes/hadronic/models/qmd/src/G4QMDMeanField.cc

@@ -392 +392 @@
       G4ThreeVector ri = system->GetParticipant( i )->GetPosition();  
       G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum();  
-      
+
       G4ThreeVector betai = p4i.v()/p4i.e();
       
-      ffr[i] = betai;
+//    R-JQMD
+      G4double Vi = GetPotential( i );
+      G4double p_zero = std::sqrt( p4i.e()*p4i.e() + 2*p4i.m()*Vi);
+      G4ThreeVector betai_R = p4i.v()/p_zero;
+      G4double mi_R = p4i.m()/p_zero;
+//
+      ffr[i] = betai_R;
       ffp[i] = G4ThreeVector( 0.0 );
+
+      if ( false ) 
+      {
+         ffr[i] = betai;
+         mi_R = 1.0;
+      }
 
       for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; j ++ )
@@ -417 +429 @@
                            * ( 1. - 2. * std::abs( jcharge - icharge ) )
                        + cl * rhc[j][i];
+         ccpp *= mi_R;
 
 /*
@@ -435 +448 @@
 
 
-           G4ThreeVector rij = ri - rj;   
-           G4ThreeVector betaij =  ( p4i + p4j ).v()/eij;   
-
-           G4ThreeVector cij = betaij - betai;   
-
-           ffr[i] = ffr[i] + 2*ccrr* ( rij + grbb*cij );
+         G4ThreeVector rij = ri - rj;   
+         G4ThreeVector betaij =  ( p4i + p4j ).v()/eij;   
+
+         G4ThreeVector cij = betaij - betai;   
+
+         ffr[i] = ffr[i] + 2*ccrr* ( rij + grbb*cij );
             
-           ffp[i] = ffp[i] - 2*ccpp* ( rij + grbb*betaij );
+         ffp[i] = ffp[i] - 2*ccpp* ( rij + grbb*betaij );
 
       }
@@ -450 +463 @@
    //std::cout << "gradu 1 " << ffr[1] << " " << ffp[1] << std::endl;
 
+}
+
+
+
+G4double G4QMDMeanField::GetPotential( G4int i )
+{
+   G4int n = system->GetTotalNumberOfParticipant();
+
+   G4double rhoa = 0.0;
+   G4double rho3 = 0.0;
+   G4double rhos = 0.0;
+   G4double rhoc = 0.0;
+
+
+   G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus();
+   G4int inuc = system->GetParticipant(i)->GetNuc();
+
+   for ( G4int j = 0 ; j < n ; j ++ )
+   {
+      G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus();
+      G4int jnuc = system->GetParticipant(j)->GetNuc();
+
+      rhoa += rha[j][i];
+      rhoc += rhe[j][i];
+      rhos += rha[j][i] * jnuc * inuc
+                * ( 1 - 2 * std::abs ( jcharge - icharge ) );
+   }
+
+   rho3 = std::pow ( rhoa , gamm );
+
+   G4double potential = c0 * rhoa
+                      + c3 * rho3
+                      + cs * rhos
+                      + cl * rhoc;
+
+   return potential;
 }
 

trunk/source/processes/hadronic/models/qmd/src/G4QMDReaction.cc

@@ -43 +43 @@
 G4QMDReaction::G4QMDReaction()
 : system ( NULL )
-, deltaT ( 1 ) // in fsec
+, deltaT ( 1 ) // in fsec (c=1)
 , maxTime ( 100 ) // will have maxTime-th time step
+, envelopF ( 1.05 ) // 10% for Peripheral reactions
 , gem ( true )
 , frag ( false )
@@ -155 +156 @@
 
 // impact parameter 
-      G4double bmax = 1.05*(bmax_0/fermi);  // 10% for Peripheral reactions
+      //G4double bmax = 1.05*(bmax_0/fermi);  // 10% for Peripheral reactions
+      G4double bmax = envelopF*(bmax_0/fermi);
       G4double b = bmax * std::sqrt ( G4UniformRand() );
 //071112