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

Timestamp:

Apr 6, 2009, 12:30:29 PM (15 years ago)

Author:

garnier

Message:

update processes

Location:

trunk/source/processes/hadronic/models/qmd

Files:

• History (modified) (1 diff)
• include/G4QMDCollision.hh (modified) (3 diffs)
• include/G4QMDMeanField.hh (modified) (2 diffs)
• include/G4QMDParticipant.hh (modified) (3 diffs)
• include/G4QMDReaction.hh (modified) (5 diffs)
• include/G4QMDSystem.hh (modified) (2 diffs)
• src/G4QMDCollision.cc (modified) (16 diffs)
• src/G4QMDGroundStateNucleus.cc (modified) (2 diffs)
• src/G4QMDMeanField.cc (modified) (4 diffs)
• src/G4QMDNucleus.cc (modified) (3 diffs)
• src/G4QMDReaction.cc (modified) (13 diffs)
• src/G4QMDSystem.cc (modified) (2 diffs)

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

@@ -13 +13 @@
      * Please list in reverse chronological order (last date on top)
      ---------------------------------------------------------------
+
+20-November-2008 Tatsumi Koi (hadr-qmd-V09-01-08)
+- Add Update 
+        include/G4QMDMeanField.hh
+        src/G4QMDMeanField.cc
+- Add "hit" flag and related methods
+        include/G4QMDParticipant.hh  
+- Add Erase and Insert Participant methods
+        include/G4QMDSystem.hh  
+        src/G4QMDSystem.cc
+- Add deltaT in signature of CalKinematicsOfBinaryCollisions
+  Add several required updating of Mean Filed 
+  Modify handling of absorption case 
+        include/G4QMDCollision.hh  
+        src/G4QMDCollision.cc  
+- Add deltaT in signature of CalKinematicsOfBinaryCollisions
+        src/G4QMDReaction.cc  
+
+07-November-2008 Tatsumi Koi (hadr-qmd-V09-01-07)
+- Add UnUseGEM and UseFRAG 
+  Fix bug in neucleon projectile
+        include/G4QMDReaction.hh
+        src/G4QMDReaction.cc  
+- Migrate to particles-V09-01-09 
+        src/G4QMDNucleus.cc  
+
+27-Oct-2008 Tatsumi Koi (hadr-qmd-V09-01-06)
+- Migrate to particles-V09-01-09 by T. Koi
+  Z,A to A,Z for functions of NucleiProperties 
+        src/G4QMDGroundStateNucleus.cc  
+        src/G4QMDNucleus.cc  
+
+24-Oct-2008 Tatsumi Koi (hadr-qmd-V09-01-05)
+- Migrate to particles-V09-01-09 by T. Koi
+        src/G4QMDGroundStateNucleus.cc  
+        src/G4QMDNucleus.cc  
+
+12-Jun-2008 Tatsumi Koi (hadr-qmd-V09-01-04)
+- Delete unnecessary dependency and unused functions
+  Change criterion of reaction by T. Koi
+        include/G4QMDReaction.hh
+        src/G4QMDReaction.cc  
+
+06-Jun-2008 Tatsumi Koi (hadr-qmd-V09-01-03)
+- Fix memory leaks by T. Koi 
+        include/G4QMDSystem.hh  
+        src/G4QMDReaction.cc  
+
+03-Jun-2008 Tatsumi Koi (hadr-qmd-V09-01-02)
+- Fix memory leaks by T. Koi 
+        src/G4QMDReaction.cc  
+        src/G4QMDCollision.cc  
+
+05-May-2008 Tatsumi Koi (hadr-qmd-V09-01-01)
+- Fixed and changed sampling method of impact parameter
+        src/G4QMDReaction.cc  
+
 26-Nov-2007 Tatsumi Koi (hadr-qmd-V09-00-05)
 - Fix of typo which introduced at (hadr-qmd-V09-00-04)

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

@@ -34 +34 @@
 //      Creation date: 9 April 2007
 // -----------------------------------------------------------------------------
+//
+// 081120 Add deltaT in signature of CalKinematicsOfBinaryCollisions
 
 #ifndef G4QMDCollision_hh
@@ -49 +51 @@
       ~G4QMDCollision();
 
-      void CalKinematicsOfBinaryCollisions();
+      void CalKinematicsOfBinaryCollisions( G4double );
       G4bool CalFinalStateOfTheBinaryCollision( G4int , G4int );
       G4bool CalFinalStateOfTheBinaryCollisionJQMD( G4double , G4double , G4ThreeVector , G4double , G4double , G4ThreeVector , G4double , G4int , G4int );
@@ -55 +57 @@
 
       void SetMeanField ( G4QMDMeanField* meanfield ){ theMeanField = meanfield; theSystem = meanfield->GetSystem(); }
-
 
    private:

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

@@ -34 +34 @@
 //      Creation date: 29 March 2007
 // -----------------------------------------------------------------------------
+// 081120 Add Update
 
 #ifndef G4QMDMeanField_hh
@@ -74 +75 @@
       std::vector< G4double > GetDepthOfPotential();  
 
+      void Update();
+
    private:
       G4double calPauliBlockingFactor( G4int ); 

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

@@ -34 +34 @@
 //      Creation date: 29 March 2007
 // -----------------------------------------------------------------------------
+//
+// 081120 Add hit flag and related methods
 
 #ifndef G4QMDParticipant_hh
@@ -70 +72 @@
 
       void UnsetInitialMark() { projectile = false; target = false; }
+      void UnsetHitMark() { hit = false; }
+      G4bool IsThisHit() { return hit; }
+      void SetHitMark() { hit = true; }
+
       void SetProjectile() { projectile = true; }
       void SetTarget() { target = true; }
@@ -82 +88 @@
       G4bool projectile; 
       G4bool target; 
+      G4bool hit; 
 };
 

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

@@ -34 +34 @@
 //      Creation date: 02 April 2007
 // -----------------------------------------------------------------------------
+//
+// 081107 Add UnUseGEM (then use the default channel of G4Evaporation)
+//            UseFrag (chage criterion of a inelastic reaction)
+//
 
 #ifndef G4QMDReaction_hh
@@ -63 +67 @@
       G4HadFinalState *ApplyYourself( const G4HadProjectile &aTrack, G4Nucleus & targetNucleus );
 
+      void UnUseGEM(){ gem = false; setEvaporationCh(); };
+      void UseFRAG(){ frag = true; };
+
    private:
-      void setInitialCondition( G4QMDSystem* , G4QMDSystem* );
+
+      void setEvaporationCh();
 
       G4QMDMeanField* meanField;
@@ -70 +78 @@
       G4QMDCollision* collision;
 
-      void doPropagation();
       void doCollision();
       std::vector< G4QMDSystem* > doClusterJudgment();
@@ -80 +87 @@
       G4Evaporation* evaporation;
       G4ExcitationHandler* excitationHandler;
+
 //      G4VPreCompoundModel* preco; 
 
@@ -102 +110 @@
       G4IonsShenCrossSection genspaXS;
 
+      G4bool gem;
+      G4bool frag;
+
 };
 

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

@@ -34 +34 @@
 //      Creation date: 29 March 2007
 // -----------------------------------------------------------------------------
+//
+// 080602 Fix memory leaks by T. Koi 
+// 081120 Add EraseParticipant and InsertParticipant Methods by T. Koi
 
 #ifndef G4QMDSystem_hh
@@ -51 +54 @@
       void SubtractSystem ( G4QMDSystem* );
 
-      void DeleteParticipant( G4int i ) { participants.erase( std::find ( participants.begin() , participants.end() , participants[ i ] ) ); };
+      G4QMDParticipant* EraseParticipant( G4int i ) { G4QMDParticipant* particle =  participants[ i ]; participants.erase( std::find ( participants.begin() , participants.end() , participants[ i ] ) ) ; return particle; };
+      void DeleteParticipant( G4int i ) { delete participants[ i ] ; participants.erase( std::find ( participants.begin() , participants.end() , participants[ i ] ) ); };
+      void InsertParticipant( G4QMDParticipant* particle , G4int j );
 
       G4int GetTotalNumberOfParticipant() { return participants.size(); };

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// 080602 Fix memory leaks by T. Koi 
+// 081120 Add deltaT in signature of CalKinematicsOfBinaryCollisions
+//        Add several required updating of Mean Filed 
+//        Modified handling of absorption case by T. Koi
+//
 #include "G4QMDCollision.hh"
 #include "G4ParticleDefinition.hh"
 #include "G4Scatterer.hh"
 #include "Randomize.hh"
-#include "G4PionZero.hh"
 
 G4QMDCollision::G4QMDCollision()
@@ -49 +53 @@
 
 
-void G4QMDCollision::CalKinematicsOfBinaryCollisions()
+void G4QMDCollision::CalKinematicsOfBinaryCollisions( G4double dt )
 {
-
+   G4double deltaT = dt; 
 
    G4int n = theSystem->GetTotalNumberOfParticipant(); 
+//081118
+   //G4int nb = 0;
+   for ( G4int i = 0 ; i < n ; i++ )
+   {
+      theSystem->GetParticipant( i )->UnsetHitMark();
+      theSystem->GetParticipant( i )->UnsetHitMark();
+      //nb += theSystem->GetParticipant( i )->GetBaryonNumber();
+   }
+   //G4cout << "nb = " << nb << " n = " << n << G4endl;
+
 
 //071101
    for ( G4int i = 0 ; i < n ; i++ )
    {
+
       //std::cout << i << " " << theSystem->GetParticipant( i )->GetDefinition()->GetParticleName() << " " << theSystem->GetParticipant( i )->GetPosition() << std::endl;
+
       if ( theSystem->GetParticipant( i )->GetDefinition()->IsShortLived() )
       {
+
+         G4bool decayed = false; 
+
          G4ParticleDefinition* pd0 = theSystem->GetParticipant( i )->GetDefinition();
          G4ThreeVector p0 = theSystem->GetParticipant( i )->GetMomentum();
@@ -67 +86 @@
          G4LorentzVector p40 = theSystem->GetParticipant( i )->Get4Momentum();
 
-         G4double epot = theMeanField->GetTotalPotential();
-         G4double eini = epot + p40.e();
+         G4double eini = theMeanField->GetTotalPotential() + p40.e();
 
          G4int n0 = theSystem->GetTotalNumberOfParticipant(); 
          G4int i0 = 0; 
-G4bool isThisEnergyOK = false;
+
+         G4bool isThisEnergyOK = false;
+
          for ( G4int ii = 0 ; ii < 4 ; ii++ )
-{ 
-
-         //G4LorentzVector p4 = theSystem->GetParticipant( i )->Get4Momentum();
-         G4LorentzVector p400 = p40;
-
-         p400 *= GeV;
-         //G4KineticTrack kt( theSystem->GetParticipant( i )->GetDefinition() , 0.0 , (theSystem->GetParticipant( i )->GetPosition())*fermi , p4 );
-         G4KineticTrack kt( pd0 , 0.0 , r0*fermi , p400 );
-//         std::cout << "G4KineticTrack " << i << " " <<  kt.GetDefinition()->GetParticleName() <<  kt.GetPosition() << std::endl;
-         G4KineticTrackVector* secs = NULL;
-         secs = kt.Decay();
-         G4int id = 0;
-         G4double et = 0;
-         if ( secs )
-         {
-            for ( G4KineticTrackVector::iterator it 
-                  = secs->begin() ; it != secs->end() ; it++ )
+         { 
+
+            //G4LorentzVector p4 = theSystem->GetParticipant( i )->Get4Momentum();
+            G4LorentzVector p400 = p40;
+
+            p400 *= GeV;
+            //G4KineticTrack kt( theSystem->GetParticipant( i )->GetDefinition() , 0.0 , (theSystem->GetParticipant( i )->GetPosition())*fermi , p4 );
+            G4KineticTrack kt( pd0 , 0.0 , r0*fermi , p400 );
+            //std::cout << "G4KineticTrack " << i << " " <<  kt.GetDefinition()->GetParticleName() <<  kt.GetPosition() << std::endl;
+            G4KineticTrackVector* secs = NULL;
+            secs = kt.Decay();
+            G4int id = 0;
+            G4double et = 0;
+            if ( secs )
             {
-//              std::cout << "G4KineticTrack" 
-//                 << " " << (*it)->GetDefinition()->GetParticleName()
-//                 << " " << (*it)->Get4Momentum()
-//                 << " " << (*it)->GetPosition()/fermi
-//                         << std::endl;
-                if ( id == 0 ) 
-                {
-                   theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() );
-                   theSystem->GetParticipant( i )->SetMomentum( (*it)->Get4Momentum().v()/GeV );
-                   theSystem->GetParticipant( i )->SetPosition( (*it)->GetPosition()/fermi );
-                   //theMeanField->Cal2BodyQuantities( i ); 
-                   et += (*it)->Get4Momentum().e()/GeV;
-                }
-                if ( id > 0 )
-                {
-                   // Append end;
-                   theSystem->SetParticipant ( new G4QMDParticipant ( (*it)->GetDefinition() , (*it)->Get4Momentum().v()/GeV , (*it)->GetPosition()/fermi ) );
-                   et += (*it)->Get4Momentum().e()/GeV;
-                   if ( id > 1 )
+               for ( G4KineticTrackVector::iterator it 
+                     = secs->begin() ; it != secs->end() ; it++ )
+               {
+/*
+                  G4cout << "G4KineticTrack" 
+                  << " " << (*it)->GetDefinition()->GetParticleName()
+                  << " " << (*it)->Get4Momentum()
+                  << " " << (*it)->GetPosition()/fermi
+                  << G4endl;
+*/
+                   if ( id == 0 ) 
                    {
- //                     std::cout << "NAGISA id >2; id= " << id  << std::endl; 
+                      theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() );
+                      theSystem->GetParticipant( i )->SetMomentum( (*it)->Get4Momentum().v()/GeV );
+                      theSystem->GetParticipant( i )->SetPosition( (*it)->GetPosition()/fermi );
+                      //theMeanField->Cal2BodyQuantities( i ); 
+                      et += (*it)->Get4Momentum().e()/GeV;
                    }
-                }
-                id++; 
-
-                delete *it;
+                   if ( id > 0 )
+                   {
+                      // Append end;
+                      theSystem->SetParticipant ( new G4QMDParticipant ( (*it)->GetDefinition() , (*it)->Get4Momentum().v()/GeV , (*it)->GetPosition()/fermi ) );
+                      et += (*it)->Get4Momentum().e()/GeV;
+                      if ( id > 1 )
+                      {
+                         //081118
+                         //G4cout << "G4QMDCollision id >2; id= " << id  << G4endl; 
+                      }
+                   }
+                   id++; // number of daughter particles
+
+                   delete *it;
+               }
+
+               theMeanField->Update();
+               i0 = id-1; // 0 enter to i
+
+               delete secs;
             }
-            theMeanField->SetSystem ( theSystem );
-            i0 = id-1; // 0 enter to i
-         }
-
-//       EnergyCheck  
-
-         G4double epot = theMeanField->GetTotalPotential();
-         G4double efin = epot + et; 
-         //std::cout <<  std::abs ( eini - efin ) - epse << std::endl; 
-//         std::cout <<  std::abs ( eini - efin ) - epse*10 << std::endl; 
-//071031
-//                                        *10 TK  
-         if ( std::abs ( eini - efin ) < epse*10 ) 
-         {
-            // Energy OK 
-//            std::cout << "Decay Succeeded Energy OK" << std::endl;
-            isThisEnergyOK = true;
-            break; 
-         }
-         else
-         {
-            for ( G4int i0i = 0 ; i0i < id-1 ; i0i++ )
+
+//          EnergyCheck  
+
+            G4double efin = theMeanField->GetTotalPotential() + et; 
+            //std::cout <<  std::abs ( eini - efin ) - epse << std::endl; 
+//            std::cout <<  std::abs ( eini - efin ) - epse*10 << std::endl; 
+//                                           *10 TK  
+            if ( std::abs ( eini - efin ) < epse*10 ) 
             {
-//               std::cout << "Decay Energitically Blocked deleteing " << i0i+n0  << std::endl;
-               theSystem->DeleteParticipant( i0i+n0 );
+               // Energy OK 
+               isThisEnergyOK = true;
+               break; 
             }
-         }
-}
+            else
+            {
+
+               theSystem->GetParticipant( i )->SetDefinition( pd0 );
+               theSystem->GetParticipant( i )->SetPosition( r0 );
+               theSystem->GetParticipant( i )->SetMomentum( p0 );
+
+               for ( G4int i0i = 0 ; i0i < id-1 ; i0i++ )
+               {
+                  //081118
+                  //std::cout << "Decay Energitically Blocked deleteing " << i0i+n0  << std::endl;
+                  theSystem->DeleteParticipant( i0i+n0 );
+               }
+               //081103
+               theMeanField->Update();
+            }
+
+         }
+
 
 //       Pauli Check 
          if ( isThisEnergyOK == true )
          {
-//          if ( theMeanField->IsPauliBlocked ( i ) != true ) 
+            if ( theMeanField->IsPauliBlocked ( i ) != true ) 
             {
-               bool allOK = true; 
+
+               G4bool allOK = true; 
                for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
                {
@@ -163 +197 @@
                }
 
-//               if ( allOK ) std::cout << "Decay Succeeded" << std::endl;
-               if ( allOK ) continue; //Do not Pauli Blocked
+               if ( allOK ) 
+               {
+                  decayed = true; //Decay Succeeded
+               }
             }
+
          }
 //       
 
-//         std::cout << "Decay Blocked" << std::endl;
-         theSystem->GetParticipant( i )->SetDefinition( pd0 );
-         theSystem->GetParticipant( i )->SetPosition( r0 );
-         theSystem->GetParticipant( i )->SetMomentum( p0 );
-
-         if ( isThisEnergyOK == true )
-         {
-         for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
-         {
-//            std::cout << "Decay Blocked deleteing " << i0i+n0  << std::endl;
-            theSystem->DeleteParticipant( i0i+n0 );
-         }
-         }
-         
-      }
-   }
+         if ( decayed )
+         {
+            //081119
+            //G4cout << "Decay Suceeded! " << std::endl;
+            theSystem->GetParticipant( i )->SetHitMark();
+            for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
+            {
+                theSystem->GetParticipant( i0i+n0 )->SetHitMark();
+            }
+
+         }
+         else
+         {
+
+//          Decay Blocked and re-enter orginal participant;
+
+            if ( isThisEnergyOK == true )  // for false case already done
+            {
+
+               theSystem->GetParticipant( i )->SetDefinition( pd0 );
+               theSystem->GetParticipant( i )->SetPosition( r0 );
+               theSystem->GetParticipant( i )->SetMomentum( p0 );
+
+               for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
+               {
+                  //081118
+                  //std::cout << "Decay Blocked deleteing " << i0i+n0  << std::endl;
+                  theSystem->DeleteParticipant( i0i+n0 );
+               }
+               //081103
+               theMeanField->Update();
+            }
+
+         }
+
+      }  //shortlive
+   }  // go next participant 
 //071101
 
 
    n = theSystem->GetTotalNumberOfParticipant(); 
-   //std::cout << "Collision n " << n << std::endl;
-
-   std::vector< G4bool > isCollided ( n , false );
-
-   for ( G4int i = 1 ; i < n ; i++ )
+
+//081118
+   //for ( G4int i = 1 ; i < n ; i++ )
+   for ( G4int i = 1 ; i < theSystem->GetTotalNumberOfParticipant() ; i++ )
    {
 
       //std::cout << "Collision i " << i << std::endl;
+      if ( theSystem->GetParticipant( i )->IsThisHit() ) continue; 
 
       G4ThreeVector ri =  theSystem->GetParticipant( i )->GetPosition();
@@ -206 +264 @@
       for ( G4int j = 0 ; j < i ; j++ )
       {
-//         std::cout << "Collision " << i << " " << j << std::endl;
+
 
 /*
@@ -216 +274 @@
 
          // Only 1 Collision allowed for each particle in a time step. 
-         if ( isCollided[ i ] == true ) continue;
-         if ( isCollided[ j ] == true ) continue;
+         //081119
+         if ( theSystem->GetParticipant( i )->IsThisHit() ) continue; 
+         if ( theSystem->GetParticipant( j )->IsThisHit() ) continue; 
+
+         //std::cout << "Collision " << i << " " << j << std::endl;
 
          // Do not allow collision between nucleons in target/projectile til its first collision.
@@ -285 +346 @@
          G4double tj = (-pjdr/rmj - bji / aij ) * p4j.e() / rmj;
 
-         G4double deltaT = 0.0;
-         deltaT = 1.0; // TK  
 
 /*
@@ -319 +378 @@
          //G4bool energetically_forbidden = !( CalFinalStateOfTheBinaryCollisionJQMD ( sig , cutoff , pcm , prcm , srt, beta , gamma , i , j ) ); // JQMD Elastic 
 
+/*
          G4bool pauli_blocked = false;
-         if ( energetically_forbidden != true ) 
+         if ( energetically_forbidden == false ) // result true 
          { 
             if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) 
@@ -330 +390 @@
          else
          {
+            if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) 
+               pauli_blocked = false;
             //std::cout << "G4QMDRESULT Collsion Blocked " << std::endl;
          } 
+*/
 
 /*
@@ -338 +401 @@
                       << std::endl;
 */
-
- 
-         if ( energetically_forbidden == true ||  pauli_blocked == true )
-         {
-//       Collsion not allowed then re enter orginal participants 
-//       Now only momentum, becasuse we only consider elastic scattering of nucleons
+//       081118
+         //if ( energetically_forbidden == true || pauli_blocked == true )
+         if ( energetically_forbidden == true )
+         {
+
+            //G4cout << " energetically_forbidden  " << G4endl;
+//          Collsion not allowed then re enter orginal participants 
+//          Now only momentum, becasuse we only consider elastic scattering of nucleons
 
             theSystem->GetParticipant( i )->SetMomentum( p4i.vect() );
             theSystem->GetParticipant( i )->SetDefinition( pdi );
             theSystem->GetParticipant( i )->SetPosition( ri );
+
             theSystem->GetParticipant( j )->SetMomentum( p4j.vect() );
             theSystem->GetParticipant( j )->SetDefinition( pdj );
             theSystem->GetParticipant( j )->SetPosition( rj );
 
+            theMeanField->Cal2BodyQuantities( i ); 
+            theMeanField->Cal2BodyQuantities( j ); 
+
          }
          else 
          {
-//       Collsion allowed (really happened) 
+
+            
+           G4bool absorption = false; 
+           if ( n == theSystem->GetTotalNumberOfParticipant()+1 ) absorption = true;
+           if ( absorption ) 
+           {
+              //G4cout << "Absorption happend " << G4endl; 
+              i = i-1; 
+              n = n-1;
+           } 
+              
+//          Collsion allowed (really happened) 
 
             // Unset Projectile/Target flag
             theSystem->GetParticipant( i )->UnsetInitialMark();
-            theSystem->GetParticipant( j )->UnsetInitialMark();
-
-            isCollided[ i ] = true; 
-            isCollided[ j ] = true; 
+            if ( !absorption ) theSystem->GetParticipant( j )->UnsetInitialMark();
+
+            theSystem->GetParticipant( i )->SetHitMark();
+            if ( !absorption ) theSystem->GetParticipant( j )->SetHitMark();
 
             theSystem->IncrementCollisionCounter();
@@ -390 +470 @@
                       << std::endl;
 */
-
-         }
-
-//         theMeanField
-
-      }
+             
+
+         }
+
+      }
+
    }
 
-//071106
-   n = theSystem->GetTotalNumberOfParticipant(); 
-   G4bool isThisModefied = false;
-   for ( G4int i = 0 ; i < n ; i++ )
-   {
-      if ( theSystem->GetParticipant( i )->GetDefinition() == G4PionZero::PionZero() )
-      {
-         if ( theSystem->GetParticipant( i )->GetPosition().mag() > 1.0e9 )
-         {
-//            std::cout << "Deleting " << i << " " << theSystem->GetParticipant( i )->GetPosition().mag() << std::endl;
-            theSystem->DeleteParticipant( i );
-            isThisModefied = true;
-         }
-      } 
-   }
-   if ( isThisModefied == true ) theMeanField->SetSystem ( theSystem );
-//071106
 
 }
@@ -423 +486 @@
 {
 
-   //G4cout << "CalFinalStateOfTheBinaryCollision " << G4endl;
-
-   G4bool result = true;
-
-   G4LorentzVector p4i =  theSystem->GetParticipant( i )->Get4Momentum();
-   G4LorentzVector p4j =  theSystem->GetParticipant( j )->Get4Momentum();
+//081103
+   //G4cout << "CalFinalStateOfTheBinaryCollision " << i << " " << j << " " << theSystem->GetTotalNumberOfParticipant() << G4endl;
+
+   G4bool result = false;
+   G4bool energyOK = false; 
+   G4bool pauliOK = false; 
+   G4bool abs = false;
+   G4QMDParticipant* absorbed = NULL;  
+
+   G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum();
+   G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum();
+
+//071031
+
+   G4double epot = theMeanField->GetTotalPotential();
+
+   G4double eini = epot + p4i.e() + p4j.e();
 
 //071031
    // will use KineticTrack
-   G4LorentzVector p4ix = p4i*GeV;
-   G4LorentzVector p4jx = p4j*GeV;
-   G4ThreeVector rix = (theSystem->GetParticipant( i )->GetPosition())*fermi; 
-   G4ThreeVector rjx = (theSystem->GetParticipant( j )->GetPosition())*fermi; 
-//071031
-
-   G4double epot = theMeanField->GetTotalPotential();
-
-   G4double eini = epot + p4i.e() + p4j.e();
-
-
-//071031
    G4ParticleDefinition* pdi0 =theSystem->GetParticipant( i )->GetDefinition();
    G4ParticleDefinition* pdj0 =theSystem->GetParticipant( j )->GetDefinition();
-   G4ThreeVector ri0 =(theSystem->GetParticipant( i )->GetPosition())*fermi;
-   G4ThreeVector rj0 =(theSystem->GetParticipant( j )->GetPosition())*fermi;
+   G4LorentzVector p4i0 = p4i*GeV;
+   G4LorentzVector p4j0 = p4j*GeV;
+   G4ThreeVector ri0 = ( theSystem->GetParticipant( i )->GetPosition() )*fermi;
+   G4ThreeVector rj0 = ( theSystem->GetParticipant( j )->GetPosition() )*fermi;
 
    for ( G4int iitry = 0 ; iitry < 4 ; iitry++ )
    {
 
-      G4KineticTrack kt1( pdi0 , 0.0 , ri0 , p4ix );
-      G4KineticTrack kt2( pdj0 , 0.0 , rj0 , p4jx );
+      abs = false;
+
+      G4KineticTrack kt1( pdi0 , 0.0 , ri0 , p4i0 );
+      G4KineticTrack kt2( pdj0 , 0.0 , rj0 , p4j0 );
+
       G4LorentzVector p4ix_new; 
       G4LorentzVector p4jx_new; 
@@ -462 +529 @@
       //std::cout << "G4QMDSCATTERER BEFORE " << kt2.GetDefinition()->GetParticleName() << " " << kt2.Get4Momentum()/GeV << " " << kt2.GetPosition()/fermi << std::endl;
       //std::cout << "THESCATTERER " << theScatterer->GetCrossSection ( kt1 , kt2 )/millibarn << " " << elastic << " " << sig << std::endl;
+
+
       if ( secs )
       {
@@ -488 +557 @@
             }
          }
-         else
-         {
-            //std::cout << "NAGISA pion absrorption " << secs->front()->GetDefinition()->GetParticleName() << std::endl;
+         else if ( secs->size() == 1 )
+         {
+//081118
+            abs = true;
+            //G4cout << "G4QMDCollision pion absrorption " << secs->front()->GetDefinition()->GetParticleName() << G4endl;
             //secs->front()->Decay();
             theSystem->GetParticipant( i )->SetDefinition( secs->front()->GetDefinition() );
             p4ix_new = secs->front()->Get4Momentum()/GeV;
             theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() );
-            
-              //std::cout << "THESCATTERER " << (*it)->GetDefinition()->GetParticleName() << std::endl;
-              p4jx_new( 0 ); 
-              //theSystem->GetParticipant( j )->SetDefinition( G4Gamma::Gamma() );
-              //theSystem->GetParticipant( j )->SetDefinition( G4Neutron::Neutron() );
-              theSystem->GetParticipant( j )->SetDefinition( G4PionZero::PionZero() );
-              theSystem->GetParticipant( j )->SetMomentum( G4ThreeVector( G4UniformRand() )*eV );
-              theSystem->GetParticipant( j )->SetPosition( G4ThreeVector( 1000, 1000, 1000 )*km );
 
          } 
 
-         if ( secs->size() > 2 ) std::cout << "NAGISA secs size > 2;  " << secs->size() << std::endl;
+//081118
+         if ( secs->size() > 2 ) 
+         {
+
+            G4cout << "G4QMDCollision secs size > 2;  " << secs->size() << G4endl;
+
+            for ( G4KineticTrackVector::iterator it 
+                = secs->begin() ; it != secs->end() ; it++ )
+            {
+               G4cout << "G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() << " " << (*it)->Get4Momentum()/GeV << G4endl;
+            }
+
+         }
 
          // deleteing KineticTrack
@@ -514 +589 @@
             delete *it;
          }
+
+         delete secs;
       }
 //071031
 
-      theMeanField->Cal2BodyQuantities( i ); 
-      theMeanField->Cal2BodyQuantities( j ); 
+      if ( !abs )
+      { 
+         theMeanField->Cal2BodyQuantities( i ); 
+         theMeanField->Cal2BodyQuantities( j ); 
+      } 
+      else
+      {
+         absorbed = theSystem->EraseParticipant( j ); 
+         theMeanField->Update();
+      }
 
       epot = theMeanField->GetTotalPotential();
@@ -540 +625 @@
          //std::cout << "collisions after " << theSystem->GetParticipant( i )->Get4Momentum() << " " << theSystem->GetParticipant( j )->Get4Momentum() << std::endl;
          //std::cout << "collisions dif " << ( p4i + p4j ) - ( theSystem->GetParticipant( i )->Get4Momentum() + theSystem->GetParticipant( j )->Get4Momentum() ) << std::endl;
-         //std::cout << "collisions before " << rix/fermi << " " << rjx/fermi << std::endl;
+         //std::cout << "collisions before " << ri0/fermi << " " << rj0/fermi << std::endl;
          //std::cout << "collisions after " << theSystem->GetParticipant( i )->GetPosition() << " " << theSystem->GetParticipant( j )->GetPosition() << std::endl;
+         energyOK = true;
+         break;
+      }
+      else
+      {
+         //G4cout << "Energy Not OK " << G4endl;
+         if ( abs )
+         {
+            //G4cout << "TKDB reinsert j " << G4endl;
+            theSystem->InsertParticipant( absorbed , j );   
+            theMeanField->Update();
+         }
+         // do not need reinsert in no absroption case 
       }
 //071031
-
-      if ( std::abs ( eini - efin ) < epse ) return result;  // Collison OK 
-
    }
 
-//  Energetically forbidden collision
-    result = false;
-
-    return result;
+// Energetically forbidden collision
+
+   if ( energyOK )
+   {
+      // Pauli Check 
+      //G4cout << "Pauli Checking " << theSystem->GetTotalNumberOfParticipant() << G4endl;
+      if ( !abs ) 
+      {
+         if ( !( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) ) 
+         {
+            //G4cout << "Binary Collision Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl;
+            pauliOK = true;
+         }
+      }
+      else 
+      {
+         if ( theMeanField->IsPauliBlocked ( i ) == false ) 
+         { 
+            //G4cout << "Absorption Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl;
+            delete absorbed;
+            pauliOK = true;
+         }
+      }
+      
+
+      if ( pauliOK ) 
+      {
+         result = true;
+      }
+      else
+      {
+         //G4cout << "Pauli Blocked" << G4endl;
+         if ( abs )
+         {
+            //G4cout << "TKDB reinsert j pauli block" << G4endl;
+            theSystem->InsertParticipant( absorbed , j );   
+            theMeanField->Update(); 
+         }
+      }
+   }
+
+   return result;
 
 } 

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
+//
 #include "G4QMDGroundStateNucleus.hh"
 
@@ -119 +121 @@
    //std::cout << "G4QMDGroundStateNucleus::packNucleons" << std::endl;
 
-   ebini = - ( G4NucleiPropertiesTable::GetBindingEnergy( GetAtomicNumber() , GetMassNumber() ) ) / GetMassNumber();
+   ebini = - G4NucleiProperties::GetBindingEnergy( GetMassNumber() , GetAtomicNumber() ) / GetMassNumber();
 
    G4double ebin00 = ebini * 0.001;

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

@@ -23 +23 @@
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
+//
+// 081120 Add Update by T. Koi
 //
 #include "G4QMDMeanField.hh"
@@ -861 +863 @@
    {
 
-      //std::cout << "Add Participants to cluseter " << it->second << std::endl;
+      //G4cout << "Add Participants to cluseter " << it->second << G4endl;
 
       if ( it->first != 0 ) 
@@ -873 +875 @@
             {
                nucleus->SetParticipant( system->GetParticipant ( itt->second ) );
-               //std::cout << "Add Participants " << itt->second << " "  << system->GetParticipant ( itt->second )->GetPosition() << std::endl;
+               //G4cout << "Add Participants " << itt->second << " "  << system->GetParticipant ( itt->second )->GetPosition() << G4endl;
              }
 
@@ -893 +895 @@
    
 }
+
+
+
+void G4QMDMeanField::Update() 
+{ 
+   SetSystem( system ); 
+}

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
+//
 #include "G4QMDNucleus.hh"
 #include "G4Proton.hh"
@@ -91 +93 @@
 {
 
-   G4double mass = G4NucleiPropertiesTable::GetNuclearMass( GetAtomicNumber() , GetMassNumber() );
+   G4double mass = G4NucleiProperties::GetNuclearMass( GetMassNumber() , GetAtomicNumber() );
    
    if ( mass == 0.0 )
@@ -214 +216 @@
    //G4cout << "PotentialEnergy in GeV " << potentialEnergy << G4endl;
    //G4cout << "BindingEnergy in GeV " << bindingEnergy << G4endl;
-   //G4cout << "G4BindingEnergy in GeV " << G4NucleiPropertiesTable::GetBindingEnergy( GetAtomicNumber() , GetMassNumber() )/GeV << G4endl;
-
-   excitationEnergy = bindingEnergy + G4NucleiPropertiesTable::GetBindingEnergy( GetAtomicNumber() , GetMassNumber() )/GeV;
+   //G4cout << "G4BindingEnergy in GeV " << G4NucleiProperties::GetBindingEnergy( GetAtomicNumber() , GetMassNumber() )/GeV << G4endl;
+
+   excitationEnergy = bindingEnergy + G4NucleiProperties::GetBindingEnergy( GetMassNumber() , GetAtomicNumber() )/GeV;
    //G4cout << "excitationEnergy in GeV " << excitationEnergy << G4endl;
    if ( excitationEnergy < 0 ) excitationEnergy = 0.0; 

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// 080505 Fixed and changed sampling method of impact parameter by T. Koi 
+// 080602 Fix memory leaks by T. Koi 
+// 080612 Delete unnecessary dependency and unused functions
+//        Change criterion of reaction by T. Koi
+// 081107 Add UnUseGEM (then use the default channel of G4Evaporation)
+//            UseFrag (chage criterion of a inelastic reaction)
+//        Fix bug in nucleon projectiles  by T. Koi    
+//
 #include "G4QMDReaction.hh"
 #include "G4QMDNucleus.hh"
-
 #include "G4QMDGroundStateNucleus.hh"
-#include "G4Fancy3DNucleus.hh"
 
 #include "G4NistManager.hh"
 
 G4QMDReaction::G4QMDReaction()
-:system ( 0 )
+: system ( NULL )
 , deltaT ( 1 ) // in fsec
 , maxTime ( 100 ) // will have maxTime-th time step
+, gem ( true )
+, frag ( false )
 {
    meanField = new G4QMDMeanField();
    collision = new G4QMDCollision();
 
+   excitationHandler = new G4ExcitationHandler;
    evaporation = new G4Evaporation;
-   evaporation->SetGEMChannel();
-   excitationHandler = new G4ExcitationHandler;
    excitationHandler->SetEvaporation( evaporation );
-//   preco = new G4PreCompoundModel( excitationHandler );
-
+   setEvaporationCh();
 }
 
@@ -53 +59 @@
 {
    delete evaporation; 
+   delete excitationHandler;
    delete collision;
    delete meanField;
-}
-
-
-
-void G4QMDReaction::setInitialCondition( G4QMDSystem* , G4QMDSystem* )
-{
-   ;
-}
-
-
-
-void G4QMDReaction::doPropagation()
-{
-   ;
 }
 
@@ -108 +101 @@
      //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*2;
+     G4double bmax_0 = std::sqrt( xs_0 / pi );
      //std::cout << "bmax_0 in fm (fermi) " <<  bmax_0/fermi << std::endl;
 
@@ -148 +141 @@
 // impact parameter 
       G4double bmax = 1.05*(bmax_0/fermi);  // 10% for Peripheral reactions
-      G4double b = bmax * G4UniformRand();
+      G4double b = bmax * std::sqrt ( G4UniformRand() );
 //071112
       //G4double b = 0;
@@ -164 +157 @@
       G4LorentzVector proj4pLAB = projectile.Get4Momentum()/GeV;
 
-      G4QMDNucleus* proj(NULL); 
-      if ( projectile.GetDefinition()->GetParticleType() == "nucleus" )
-      {
-         proj = new G4QMDNucleus; 
+      G4QMDGroundStateNucleus* proj(NULL); 
+      if ( projectile.GetDefinition()->GetParticleType() == "nucleus" 
+        || projectile.GetDefinition()->GetParticleName() == "proton"
+        || projectile.GetDefinition()->GetParticleName() == "neutron" )
+      {
 
          proj_Z = proj_pd->GetAtomicNumber();
@@ -185 +179 @@
       G4int ia = int ( target.GetN() );
 
-      //G4QMDNucleus* targ = new G4QMDNucleus;
-      G4QMDNucleus* targ = new G4QMDGroundStateNucleus( iz , ia );
+      G4QMDGroundStateNucleus* targ = new G4QMDGroundStateNucleus( iz , ia );
 
       meanField->SetSystem (targ );
@@ -283 +276 @@
       meanField->DoPropagation( deltaT );
       //system->ShowParticipants();
-      collision->CalKinematicsOfBinaryCollisions();
+      collision->CalKinematicsOfBinaryCollisions( deltaT );
 
       if ( i / 10 * 10 == i ) 
@@ -385 +378 @@
          //if ( elasticLike_energy == true && withCollision == true ) elastic = true;   // ielst = 1 
 //       InelasticLike without Collision 
-         if ( elasticLike_energy == false ) elastic = false;                          // ielst = 2                
+         //if ( elasticLike_energy == false ) elastic = false;                          // ielst = 2                
+         if ( frag == true )
+            if ( elasticLike_energy == false ) elastic = false;
 //       InelasticLike with Collision 
-         //if ( elasticLike_energy == false && withCollision == true ) elastic = false; // ielst = 3
+         if ( elasticLike_energy == false && withCollision == true ) elastic = false; // ielst = 3
 
       }
@@ -403 +398 @@
       //G4cout << elastic << G4endl;
       // if elastic is true try again from sampling of impact parameter 
+
+      if ( elastic == true )
+      {
+         // delete this nucleues
+         for ( std::vector< G4QMDNucleus* >::iterator
+               it = nucleuses.begin() ; it != nucleuses.end() ; it++ )
+         {
+            delete *it;
+         }
+         nucleuses.clear();
+      }
    } 
 
@@ -436 +442 @@
          for ( G4int i = 0 ; i < (*it)->GetTotalNumberOfParticipant() ; i++ )
          {
-            system->SetParticipant ( (*it)->GetParticipant( i ) );  
+            G4QMDParticipant* aP = new G4QMDParticipant( ( (*it)->GetParticipant( i ) )->GetDefinition() , ( (*it)->GetParticipant( i ) )->GetMomentum() , ( (*it)->GetParticipant( i ) )->GetPosition() );  
+            system->SetParticipant ( aP );  
          } 
          continue;  
@@ -508 +515 @@
       }
 
+      for ( G4ReactionProductVector::iterator itt
+          = rv->begin() ; itt != rv->end() ; itt++ )
+      {
+          delete *itt;
+      }
+      delete rv;
+
       delete aFragment;
 
-      delete *it;  // delete nulceuse 
-
    }
+
 
 
@@ -534 +547 @@
 
    }
+
+   for ( std::vector< G4QMDNucleus* >::iterator it
+       = nucleuses.begin() ; it != nucleuses.end() ; it++ )
+   {
+      delete *it;  // delete nulceuse 
+   }
+   nucleuses.clear();
 
    system->Clear();
@@ -666 +686 @@
 
 }
+
+
+
+void G4QMDReaction::setEvaporationCh()
+{
+
+   if ( gem == true ) 
+      evaporation->SetGEMChannel();
+   else
+      evaporation->SetDefaultChannel();
+
+}
+

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

@@ -24 +24 @@
 // ********************************************************************
 //
+// 081120 Add InsertParticipant Method by T. Koi
+
 #include "G4QMDSystem.hh"
 #include <iomanip>
@@ -97 +99 @@
    G4cout << "Sum upped Momentum and mag " << p_sum << " " << p_sum.mag() << G4endl;
 }
+
+
+
+void G4QMDSystem::InsertParticipant ( G4QMDParticipant* particle , G4int n )
+{
+
+   if ( (size_t) n > participants.size()+1 )
+      G4cout << "G4QMDSystem::InsertParticipant size error" << G4endl;
+
+   std::vector< G4QMDParticipant* >::iterator it; 
+   it = participants.begin();
+
+   for ( G4int i = 0; i < n ; i++ )
+      it++;
+
+   participants.insert( it, particle );
+}