Index: trunk/source/processes/hadronic/models/qmd/src/G4QMDCollision.cc
===================================================================
--- trunk/source/processes/hadronic/models/qmd/src/G4QMDCollision.cc	(revision 962)
+++ trunk/source/processes/hadronic/models/qmd/src/G4QMDCollision.cc	(revision 1055)
@@ -28,7 +28,8 @@
 //        Add several required updating of Mean Filed 
 //        Modified handling of absorption case by T. Koi
+// 090126 Fix in absorption case by T. Koi  
+// 090331 Fix for gamma participant by T. Koi 
 //
 #include "G4QMDCollision.hh"
-#include "G4ParticleDefinition.hh"
 #include "G4Scatterer.hh"
 #include "Randomize.hh"
@@ -260,4 +261,6 @@
       G4double rmi =  theSystem->GetParticipant( i )->GetMass();
       G4ParticleDefinition* pdi =  theSystem->GetParticipant( i )->GetDefinition();
+//090331 gamma 
+      if ( pdi->GetPDGMass() == 0.0 ) continue;
 
       //std::cout << " p4i00 " << p4i << std::endl;
@@ -295,4 +298,6 @@
          G4double rmj =  theSystem->GetParticipant( j )->GetMass();
          G4ParticleDefinition* pdj =  theSystem->GetParticipant( j )->GetDefinition();
+//090331 gamma 
+         if ( pdj->GetPDGMass() == 0.0 ) continue;
 
          G4double rr2 = theMeanField->GetRR2( i , j );
@@ -332,5 +337,5 @@
          G4double pidr = p4i.vect()*dr;
          G4double pjdr = p4j.vect()*dr;
-  
+
          G4double aij = 1.0 - ( rmi*rmj /pij ) * ( rmi*rmj /pij ); 
          G4double bij = pidr / rmi - pjdr*rmi/pij;
@@ -660,5 +665,7 @@
       else 
       {
-         if ( theMeanField->IsPauliBlocked ( i ) == false ) 
+         //if ( theMeanField->IsPauliBlocked ( i ) == false ) 
+         //090126                            i-1 cause jth is erased
+         if ( theMeanField->IsPauliBlocked ( i-1 ) == false ) 
          { 
             //G4cout << "Absorption Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl;
Index: trunk/source/processes/hadronic/models/qmd/src/G4QMDReaction.cc
===================================================================
--- trunk/source/processes/hadronic/models/qmd/src/G4QMDReaction.cc	(revision 962)
+++ trunk/source/processes/hadronic/models/qmd/src/G4QMDReaction.cc	(revision 1055)
@@ -31,4 +31,6 @@
 //            UseFrag (chage criterion of a inelastic reaction)
 //        Fix bug in nucleon projectiles  by T. Koi    
+// 090122 Be8 -> Alpha + Alpha 
+// 090331 Change member shenXS and genspaXS object to pointer 
 //
 #include "G4QMDReaction.hh"
@@ -45,4 +47,8 @@
 , frag ( false )
 {
+
+   //090331
+   shenXS = new G4IonsShenCrossSection();
+   //genspaXS = new G4GeneralSpaceNNCrossSection();
    meanField = new G4QMDMeanField();
    collision = new G4QMDCollision();
@@ -99,6 +105,7 @@
      G4double aTemp = projectile.GetMaterial()->GetTemperature();
 
-     //G4double xs_0 = shenXS.GetCrossSection ( proj_dp , targ_ele , aTemp );
-     G4double xs_0 = genspaXS.GetCrossSection ( proj_dp , targ_ele , aTemp );
+     //090331
+     G4double xs_0 = shenXS->GetCrossSection ( proj_dp , targ_ele , aTemp );
+     //G4double xs_0 = genspaXS->GetCrossSection ( proj_dp , targ_ele , aTemp );
      G4double bmax_0 = std::sqrt( xs_0 / pi );
      //std::cout << "bmax_0 in fm (fermi) " <<  bmax_0/fermi << std::endl;
@@ -470,10 +477,42 @@
           // Secondary from this nucleus (*it) 
           G4ParticleDefinition* pd = (*itt)->GetDefinition();
+
           G4LorentzVector p4 ( (*itt)->GetMomentum()/GeV , (*itt)->GetTotalEnergy()/GeV );  //in nucleus(*it) rest system
           G4LorentzVector p4_CM = CLHEP::boostOf( p4 , -nucleus_p4CM.findBoostToCM() );  // Back to CM
           G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB ); // Back to LAB  
 
-          G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );  
-          theParticleChange.AddSecondary( dp ); 
+
+//090122
+          //theParticleChange.AddSecondary( dp ); 
+          if ( !( pd->GetAtomicNumber() == 4 && pd->GetAtomicMass() == 8 ) )
+          {
+             G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );  
+             theParticleChange.AddSecondary( dp ); 
+          }
+          else
+          {
+             //Be8 -> Alpha + Alpha + Q
+             G4ThreeVector randomized_direction( G4UniformRand() , G4UniformRand() , G4UniformRand() );
+             randomized_direction = randomized_direction.unit();
+             G4double q_decay = (*itt)->GetMass() - 2*G4Alpha::Alpha()->GetPDGMass();
+             G4double p_decay = std::sqrt ( std::pow(G4Alpha::Alpha()->GetPDGMass()+q_decay/2,2) - std::pow(G4Alpha::Alpha()->GetPDGMass() , 2 ) ); 
+             G4LorentzVector p4_a1 ( p_decay*randomized_direction , G4Alpha::Alpha()->GetPDGMass()+q_decay/2 );  //in Be8 rest system
+             
+             G4LorentzVector p4_a1_Be8 = CLHEP::boostOf ( p4_a1/GeV , -p4.findBoostToCM() );
+             G4LorentzVector p4_a1_CM = CLHEP::boostOf ( p4_a1_Be8 , -nucleus_p4CM.findBoostToCM() );
+             G4LorentzVector p4_a1_LAB = CLHEP::boostOf ( p4_a1_CM , boostBackToLAB );
+
+             G4LorentzVector p4_a2 ( -p_decay*randomized_direction , G4Alpha::Alpha()->GetPDGMass()+q_decay/2 );  //in Be8 rest system
+             
+             G4LorentzVector p4_a2_Be8 = CLHEP::boostOf ( p4_a2/GeV , -p4.findBoostToCM() );
+             G4LorentzVector p4_a2_CM = CLHEP::boostOf ( p4_a2_Be8 , -nucleus_p4CM.findBoostToCM() );
+             G4LorentzVector p4_a2_LAB = CLHEP::boostOf ( p4_a2_CM , boostBackToLAB );
+             
+             G4DynamicParticle* dp1 = new G4DynamicParticle( G4Alpha::Alpha() , p4_a1_LAB*GeV );  
+             G4DynamicParticle* dp2 = new G4DynamicParticle( G4Alpha::Alpha() , p4_a2_LAB*GeV );  
+             theParticleChange.AddSecondary( dp1 ); 
+             theParticleChange.AddSecondary( dp2 ); 
+          }
+//090122
 
 /*
@@ -698,3 +737,2 @@
 
 }
-
