Changeset 1228 for trunk/source/processes/hadronic/models/qmd/src

Timestamp:

Jan 8, 2010, 11:56:51 AM (15 years ago)

Author:

garnier

Message:

update geant4.9.3 tag

File:

• trunk/source/processes/hadronic/models/qmd/src/G4QMDReaction.cc (modified) (8 diffs)

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

@@ -52 +52 @@
    shenXS = new G4IonsShenCrossSection();
    //genspaXS = new G4GeneralSpaceNNCrossSection();
+   piNucXS = new G4PiNuclearCrossSection();
    meanField = new G4QMDMeanField();
    collision = new G4QMDCollision();
@@ -107 +108 @@
 
      //090331
-     G4double xs_0 = shenXS->GetCrossSection ( proj_dp , targ_ele , aTemp );
+  
+   G4VCrossSectionDataSet* theXS = shenXS;
+
+   if ( proj_pd->GetParticleType() == "meson" ) theXS = piNucXS;
+
+   G4double xs_0 = theXS->GetCrossSection ( proj_dp , targ_ele , aTemp );
+
      //G4double xs_0 = genspaXS->GetCrossSection ( proj_dp , targ_ele , aTemp );
      G4double bmax_0 = std::sqrt( xs_0 / pi );
@@ -158 +165 @@
 
       G4double plab = projectile.GetTotalMomentum()/GeV;
-      G4double elab = (projectile.GetKineticEnergy() + proj_pd->GetPDGMass() + targ_pd->GetPDGMass() )/GeV;
+      G4double elab = ( projectile.GetKineticEnergy() + proj_pd->GetPDGMass() + targ_pd->GetPDGMass() )/GeV;
 
       calcOffSetOfCollision( b , proj_pd , targ_pd , plab , elab , bmax , boostCMtoNN );
@@ -177 +184 @@
          //proj->ShowParticipants();
 
-      }
-
-      meanField->SetSystem ( proj );
-      proj->SetTotalPotential( meanField->GetTotalPotential() );
-      proj->CalEnergyAndAngularMomentumInCM();
+
+         meanField->SetSystem ( proj );
+         proj->SetTotalPotential( meanField->GetTotalPotential() );
+         proj->CalEnergyAndAngularMomentumInCM();
+
+      }
 
 // Target
@@ -220 +228 @@
       G4LorentzVector targ4pCM = CLHEP::boostOf ( targ4p , boostToReac );
 
-
 //    Projectile
       if ( proj != NULL )
@@ -251 +258 @@
 //       projectile is particle
 
-         G4int i = targ->GetTotalNumberOfParticipant(); 
+         // avoid multiple set in "elastic" loop
+         if ( system->GetTotalNumberOfParticipant() == targ->GetTotalNumberOfParticipant() )
+         {
+
+            G4int i = targ->GetTotalNumberOfParticipant(); 
       
-         G4ThreeVector p0( 0 ); 
-         G4ThreeVector r0( 0 );
-
-
-         G4ThreeVector p ( p0.x() + coulomb_collision_px_proj 
-                         , p0.y() 
-                         , p0.z() * coulomb_collision_gamma_proj + coulomb_collision_pz_proj ); 
-
-         G4ThreeVector r ( r0.x() + coulomb_collision_rx_proj 
-                         , r0.y() 
-                         , r0.z() / coulomb_collision_gamma_proj + coulomb_collision_rz_proj ); 
-
-         system->SetParticipant( new G4QMDParticipant( (G4ParticleDefinition*)projectile.GetDefinition() , p , r ) );
-         system->GetParticipant ( i )->SetProjectile();
-      }
+            G4ThreeVector p0( 0 ); 
+            G4ThreeVector r0( 0 );
+
+            G4ThreeVector p ( p0.x() + coulomb_collision_px_proj 
+                            , p0.y() 
+                            , p0.z() * coulomb_collision_gamma_proj + coulomb_collision_pz_proj ); 
+
+            G4ThreeVector r ( r0.x() + coulomb_collision_rx_proj 
+                            , r0.y() 
+                            , r0.z() / coulomb_collision_gamma_proj + coulomb_collision_rz_proj ); 
+
+            system->SetParticipant( new G4QMDParticipant( (G4ParticleDefinition*)projectile.GetDefinition() , p , r ) );
+            // This is not important becase only 1 projectile particle.
+            system->GetParticipant ( i )->SetProjectile();
+         }
+
+      }
+      //system->ShowParticipants();
 
       delete targ;
       delete proj;
-
 
    meanField->SetSystem ( system );
@@ -611 +624 @@
 G4double ptot , G4double etot , G4double bmax , G4ThreeVector boostToCM )
 {
+
    G4double mass_proj = pd_proj->GetPDGMass()/GeV;
    G4double mass_targ = pd_targ->GetPDGMass()/GeV;
@@ -623 +637 @@
    G4double eccm = stot - ( mass_proj + mass_targ );
    
-   G4int zp = pd_proj->GetAtomicNumber();
-   G4int ap = pd_proj->GetAtomicMass();
+   G4int zp = 1;
+   G4int ap = 1;
+   if ( pd_proj->GetParticleType() == "nucleus" )
+   {
+      zp = pd_proj->GetAtomicNumber();
+      ap = pd_proj->GetAtomicMass();
+   }
+   else 
+   {
+      // proton, neutron, mesons
+      zp = int ( pd_proj->GetPDGCharge()/eplus + 0.5 );  
+      // ap = 1;
+   }
+   
+
    G4int zt = pd_targ->GetAtomicNumber();
    G4int at = pd_targ->GetAtomicMass();
+
 
    //G4double rmax0 = 8.0;  // T.K dicide parameter value  // for low energy