Ignore:
Timestamp:
Apr 6, 2009, 12:30:29 PM (15 years ago)
Author:
garnier
Message:

update processes

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r819 r962  
    2424// ********************************************************************
    2525//
     26// 080602 Fix memory leaks by T. Koi
     27// 081120 Add deltaT in signature of CalKinematicsOfBinaryCollisions
     28//        Add several required updating of Mean Filed
     29//        Modified handling of absorption case by T. Koi
     30//
    2631#include "G4QMDCollision.hh"
    2732#include "G4ParticleDefinition.hh"
    2833#include "G4Scatterer.hh"
    2934#include "Randomize.hh"
    30 #include "G4PionZero.hh"
    3135
    3236G4QMDCollision::G4QMDCollision()
    … …  
    4953
    5054
    51 void G4QMDCollision::CalKinematicsOfBinaryCollisions()
     55void G4QMDCollision::CalKinematicsOfBinaryCollisions( G4double dt )
    5256{
    53 
     57   G4double deltaT = dt;
    5458
    5559   G4int n = theSystem->GetTotalNumberOfParticipant();
     60//081118
     61   //G4int nb = 0;
     62   for ( G4int i = 0 ; i < n ; i++ )
     63   {
     64      theSystem->GetParticipant( i )->UnsetHitMark();
     65      theSystem->GetParticipant( i )->UnsetHitMark();
     66      //nb += theSystem->GetParticipant( i )->GetBaryonNumber();
     67   }
     68   //G4cout << "nb = " << nb << " n = " << n << G4endl;
     69
    5670
    5771//071101
    5872   for ( G4int i = 0 ; i < n ; i++ )
    5973   {
     74
    6075      //std::cout << i << " " << theSystem->GetParticipant( i )->GetDefinition()->GetParticleName() << " " << theSystem->GetParticipant( i )->GetPosition() << std::endl;
     76
    6177      if ( theSystem->GetParticipant( i )->GetDefinition()->IsShortLived() )
    6278      {
     79
     80         G4bool decayed = false;
     81
    6382         G4ParticleDefinition* pd0 = theSystem->GetParticipant( i )->GetDefinition();
    6483         G4ThreeVector p0 = theSystem->GetParticipant( i )->GetMomentum();
    … …  
    6786         G4LorentzVector p40 = theSystem->GetParticipant( i )->Get4Momentum();
    6887
    69          G4double epot = theMeanField->GetTotalPotential();
    70          G4double eini = epot + p40.e();
     88         G4double eini = theMeanField->GetTotalPotential() + p40.e();
    7189
    7290         G4int n0 = theSystem->GetTotalNumberOfParticipant();
    7391         G4int i0 = 0;
    74 G4bool isThisEnergyOK = false;
     92
     93         G4bool isThisEnergyOK = false;
     94
    7595         for ( G4int ii = 0 ; ii < 4 ; ii++ )
    76 {
    77 
    78          //G4LorentzVector p4 = theSystem->GetParticipant( i )->Get4Momentum();
    79          G4LorentzVector p400 = p40;
    80 
    81          p400 *= GeV;
    82          //G4KineticTrack kt( theSystem->GetParticipant( i )->GetDefinition() , 0.0 , (theSystem->GetParticipant( i )->GetPosition())*fermi , p4 );
    83          G4KineticTrack kt( pd0 , 0.0 , r0*fermi , p400 );
    84 //         std::cout << "G4KineticTrack " << i << " " <<  kt.GetDefinition()->GetParticleName() <<  kt.GetPosition() << std::endl;
    85          G4KineticTrackVector* secs = NULL;
    86          secs = kt.Decay();
    87          G4int id = 0;
    88          G4double et = 0;
    89          if ( secs )
    90          {
    91             for ( G4KineticTrackVector::iterator it
    92                   = secs->begin() ; it != secs->end() ; it++ )
     96         {
     97
     98            //G4LorentzVector p4 = theSystem->GetParticipant( i )->Get4Momentum();
     99            G4LorentzVector p400 = p40;
     100
     101            p400 *= GeV;
     102            //G4KineticTrack kt( theSystem->GetParticipant( i )->GetDefinition() , 0.0 , (theSystem->GetParticipant( i )->GetPosition())*fermi , p4 );
     103            G4KineticTrack kt( pd0 , 0.0 , r0*fermi , p400 );
     104            //std::cout << "G4KineticTrack " << i << " " <<  kt.GetDefinition()->GetParticleName() <<  kt.GetPosition() << std::endl;
     105            G4KineticTrackVector* secs = NULL;
     106            secs = kt.Decay();
     107            G4int id = 0;
     108            G4double et = 0;
     109            if ( secs )
    93110            {
    94 //              std::cout << "G4KineticTrack"
    95 //                 << " " << (*it)->GetDefinition()->GetParticleName()
    96 //                 << " " << (*it)->Get4Momentum()
    97 //                 << " " << (*it)->GetPosition()/fermi
    98 //                         << std::endl;
    99                 if ( id == 0 )
    100                 {
    101                    theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() );
    102                    theSystem->GetParticipant( i )->SetMomentum( (*it)->Get4Momentum().v()/GeV );
    103                    theSystem->GetParticipant( i )->SetPosition( (*it)->GetPosition()/fermi );
    104                    //theMeanField->Cal2BodyQuantities( i );
    105                    et += (*it)->Get4Momentum().e()/GeV;
    106                 }
    107                 if ( id > 0 )
    108                 {
    109                    // Append end;
    110                    theSystem->SetParticipant ( new G4QMDParticipant ( (*it)->GetDefinition() , (*it)->Get4Momentum().v()/GeV , (*it)->GetPosition()/fermi ) );
    111                    et += (*it)->Get4Momentum().e()/GeV;
    112                    if ( id > 1 )
     111               for ( G4KineticTrackVector::iterator it
     112                     = secs->begin() ; it != secs->end() ; it++ )
     113               {
     114/*
     115                  G4cout << "G4KineticTrack"
     116                  << " " << (*it)->GetDefinition()->GetParticleName()
     117                  << " " << (*it)->Get4Momentum()
     118                  << " " << (*it)->GetPosition()/fermi
     119                  << G4endl;
     120*/
     121                   if ( id == 0 )
    113122                   {
    114  //                     std::cout << "NAGISA id >2; id= " << id  << std::endl;
     123                      theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() );
     124                      theSystem->GetParticipant( i )->SetMomentum( (*it)->Get4Momentum().v()/GeV );
     125                      theSystem->GetParticipant( i )->SetPosition( (*it)->GetPosition()/fermi );
     126                      //theMeanField->Cal2BodyQuantities( i );
     127                      et += (*it)->Get4Momentum().e()/GeV;
    115128                   }
    116                 }
    117                 id++;
    118 
    119                 delete *it;
     129                   if ( id > 0 )
     130                   {
     131                      // Append end;
     132                      theSystem->SetParticipant ( new G4QMDParticipant ( (*it)->GetDefinition() , (*it)->Get4Momentum().v()/GeV , (*it)->GetPosition()/fermi ) );
     133                      et += (*it)->Get4Momentum().e()/GeV;
     134                      if ( id > 1 )
     135                      {
     136                         //081118
     137                         //G4cout << "G4QMDCollision id >2; id= " << id  << G4endl;
     138                      }
     139                   }
     140                   id++; // number of daughter particles
     141
     142                   delete *it;
     143               }
     144
     145               theMeanField->Update();
     146               i0 = id-1; // 0 enter to i
     147
     148               delete secs;
    120149            }
    121             theMeanField->SetSystem ( theSystem );
    122             i0 = id-1; // 0 enter to i
    123          }
    124 
    125 //       EnergyCheck 
    126 
    127          G4double epot = theMeanField->GetTotalPotential();
    128          G4double efin = epot + et;
    129          //std::cout <<  std::abs ( eini - efin ) - epse << std::endl;
    130 //         std::cout <<  std::abs ( eini - efin ) - epse*10 << std::endl;
    131 //071031
    132 //                                        *10 TK 
    133          if ( std::abs ( eini - efin ) < epse*10 )
    134          {
    135             // Energy OK
    136 //            std::cout << "Decay Succeeded Energy OK" << std::endl;
    137             isThisEnergyOK = true;
    138             break;
    139          }
    140          else
    141          {
    142             for ( G4int i0i = 0 ; i0i < id-1 ; i0i++ )
     150
     151//          EnergyCheck 
     152
     153            G4double efin = theMeanField->GetTotalPotential() + et;
     154            //std::cout <<  std::abs ( eini - efin ) - epse << std::endl;
     155//            std::cout <<  std::abs ( eini - efin ) - epse*10 << std::endl;
     156//                                           *10 TK 
     157            if ( std::abs ( eini - efin ) < epse*10 )
    143158            {
    144 //               std::cout << "Decay Energitically Blocked deleteing " << i0i+n0  << std::endl;
    145                theSystem->DeleteParticipant( i0i+n0 );
     159               // Energy OK
     160               isThisEnergyOK = true;
     161               break;
    146162            }
    147          }
    148 }
     163            else
     164            {
     165
     166               theSystem->GetParticipant( i )->SetDefinition( pd0 );
     167               theSystem->GetParticipant( i )->SetPosition( r0 );
     168               theSystem->GetParticipant( i )->SetMomentum( p0 );
     169
     170               for ( G4int i0i = 0 ; i0i < id-1 ; i0i++ )
     171               {
     172                  //081118
     173                  //std::cout << "Decay Energitically Blocked deleteing " << i0i+n0  << std::endl;
     174                  theSystem->DeleteParticipant( i0i+n0 );
     175               }
     176               //081103
     177               theMeanField->Update();
     178            }
     179
     180         }
     181
    149182
    150183//       Pauli Check
    151184         if ( isThisEnergyOK == true )
    152185         {
    153 //          if ( theMeanField->IsPauliBlocked ( i ) != true )
     186            if ( theMeanField->IsPauliBlocked ( i ) != true )
    154187            {
    155                bool allOK = true;
     188
     189               G4bool allOK = true;
    156190               for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
    157191               {
    … …  
    163197               }
    164198
    165 //               if ( allOK ) std::cout << "Decay Succeeded" << std::endl;
    166                if ( allOK ) continue; //Do not Pauli Blocked
     199               if ( allOK )
     200               {
     201                  decayed = true; //Decay Succeeded
     202               }
    167203            }
     204
    168205         }
    169206//       
    170207
    171 //         std::cout << "Decay Blocked" << std::endl;
    172          theSystem->GetParticipant( i )->SetDefinition( pd0 );
    173          theSystem->GetParticipant( i )->SetPosition( r0 );
    174          theSystem->GetParticipant( i )->SetMomentum( p0 );
    175 
    176          if ( isThisEnergyOK == true )
    177          {
    178          for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
    179          {
    180 //            std::cout << "Decay Blocked deleteing " << i0i+n0  << std::endl;
    181             theSystem->DeleteParticipant( i0i+n0 );
    182          }
    183          }
    184          
    185       }
    186    }
     208         if ( decayed )
     209         {
     210            //081119
     211            //G4cout << "Decay Suceeded! " << std::endl;
     212            theSystem->GetParticipant( i )->SetHitMark();
     213            for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
     214            {
     215                theSystem->GetParticipant( i0i+n0 )->SetHitMark();
     216            }
     217
     218         }
     219         else
     220         {
     221
     222//          Decay Blocked and re-enter orginal participant;
     223
     224            if ( isThisEnergyOK == true )  // for false case already done
     225            {
     226
     227               theSystem->GetParticipant( i )->SetDefinition( pd0 );
     228               theSystem->GetParticipant( i )->SetPosition( r0 );
     229               theSystem->GetParticipant( i )->SetMomentum( p0 );
     230
     231               for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
     232               {
     233                  //081118
     234                  //std::cout << "Decay Blocked deleteing " << i0i+n0  << std::endl;
     235                  theSystem->DeleteParticipant( i0i+n0 );
     236               }
     237               //081103
     238               theMeanField->Update();
     239            }
     240
     241         }
     242
     243      }  //shortlive
     244   }  // go next participant
    187245//071101
    188246
    189247
    190248   n = theSystem->GetTotalNumberOfParticipant();
    191    //std::cout << "Collision n " << n << std::endl;
    192 
    193    std::vector< G4bool > isCollided ( n , false );
    194 
    195    for ( G4int i = 1 ; i < n ; i++ )
     249
     250//081118
     251   //for ( G4int i = 1 ; i < n ; i++ )
     252   for ( G4int i = 1 ; i < theSystem->GetTotalNumberOfParticipant() ; i++ )
    196253   {
    197254
    198255      //std::cout << "Collision i " << i << std::endl;
     256      if ( theSystem->GetParticipant( i )->IsThisHit() ) continue;
    199257
    200258      G4ThreeVector ri =  theSystem->GetParticipant( i )->GetPosition();
    … …  
    206264      for ( G4int j = 0 ; j < i ; j++ )
    207265      {
    208 //         std::cout << "Collision " << i << " " << j << std::endl;
     266
    209267
    210268/*
    … …  
    216274
    217275         // Only 1 Collision allowed for each particle in a time step.
    218          if ( isCollided[ i ] == true ) continue;
    219          if ( isCollided[ j ] == true ) continue;
     276         //081119
     277         if ( theSystem->GetParticipant( i )->IsThisHit() ) continue;
     278         if ( theSystem->GetParticipant( j )->IsThisHit() ) continue;
     279
     280         //std::cout << "Collision " << i << " " << j << std::endl;
    220281
    221282         // Do not allow collision between nucleons in target/projectile til its first collision.
    … …  
    285346         G4double tj = (-pjdr/rmj - bji / aij ) * p4j.e() / rmj;
    286347
    287          G4double deltaT = 0.0;
    288          deltaT = 1.0; // TK 
    289348
    290349/*
    … …  
    319378         //G4bool energetically_forbidden = !( CalFinalStateOfTheBinaryCollisionJQMD ( sig , cutoff , pcm , prcm , srt, beta , gamma , i , j ) ); // JQMD Elastic
    320379
     380/*
    321381         G4bool pauli_blocked = false;
    322          if ( energetically_forbidden != true )
     382         if ( energetically_forbidden == false ) // result true
    323383         {
    324384            if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true )
    … …  
    330390         else
    331391         {
     392            if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true )
     393               pauli_blocked = false;
    332394            //std::cout << "G4QMDRESULT Collsion Blocked " << std::endl;
    333395         }
     396*/
    334397
    335398/*
    … …  
    338401                      << std::endl;
    339402*/
    340 
    341  
    342          if ( energetically_forbidden == true ||  pauli_blocked == true )
    343          {
    344 //       Collsion not allowed then re enter orginal participants
    345 //       Now only momentum, becasuse we only consider elastic scattering of nucleons
     403//       081118
     404         //if ( energetically_forbidden == true || pauli_blocked == true )
     405         if ( energetically_forbidden == true )
     406         {
     407
     408            //G4cout << " energetically_forbidden  " << G4endl;
     409//          Collsion not allowed then re enter orginal participants
     410//          Now only momentum, becasuse we only consider elastic scattering of nucleons
    346411
    347412            theSystem->GetParticipant( i )->SetMomentum( p4i.vect() );
    348413            theSystem->GetParticipant( i )->SetDefinition( pdi );
    349414            theSystem->GetParticipant( i )->SetPosition( ri );
     415
    350416            theSystem->GetParticipant( j )->SetMomentum( p4j.vect() );
    351417            theSystem->GetParticipant( j )->SetDefinition( pdj );
    352418            theSystem->GetParticipant( j )->SetPosition( rj );
    353419
     420            theMeanField->Cal2BodyQuantities( i );
     421            theMeanField->Cal2BodyQuantities( j );
     422
    354423         }
    355424         else
    356425         {
    357 //       Collsion allowed (really happened)
     426
     427           
     428           G4bool absorption = false;
     429           if ( n == theSystem->GetTotalNumberOfParticipant()+1 ) absorption = true;
     430           if ( absorption )
     431           {
     432              //G4cout << "Absorption happend " << G4endl;
     433              i = i-1;
     434              n = n-1;
     435           }
     436             
     437//          Collsion allowed (really happened)
    358438
    359439            // Unset Projectile/Target flag
    360440            theSystem->GetParticipant( i )->UnsetInitialMark();
    361             theSystem->GetParticipant( j )->UnsetInitialMark();
    362 
    363             isCollided[ i ] = true;
    364             isCollided[ j ] = true;
     441            if ( !absorption ) theSystem->GetParticipant( j )->UnsetInitialMark();
     442
     443            theSystem->GetParticipant( i )->SetHitMark();
     444            if ( !absorption ) theSystem->GetParticipant( j )->SetHitMark();
    365445
    366446            theSystem->IncrementCollisionCounter();
    … …  
    390470                      << std::endl;
    391471*/
    392 
    393          }
    394 
    395 //         theMeanField
    396 
    397       }
     472             
     473
     474         }
     475
     476      }
     477
    398478   }
    399479
    400 //071106
    401    n = theSystem->GetTotalNumberOfParticipant();
    402    G4bool isThisModefied = false;
    403    for ( G4int i = 0 ; i < n ; i++ )
    404    {
    405       if ( theSystem->GetParticipant( i )->GetDefinition() == G4PionZero::PionZero() )
    406       {
    407          if ( theSystem->GetParticipant( i )->GetPosition().mag() > 1.0e9 )
    408          {
    409 //            std::cout << "Deleting " << i << " " << theSystem->GetParticipant( i )->GetPosition().mag() << std::endl;
    410             theSystem->DeleteParticipant( i );
    411             isThisModefied = true;
    412          }
    413       }
    414    }
    415    if ( isThisModefied == true ) theMeanField->SetSystem ( theSystem );
    416 //071106
    417480
    418481}
    … …  
    423486{
    424487
    425    //G4cout << "CalFinalStateOfTheBinaryCollision " << G4endl;
    426 
    427    G4bool result = true;
    428 
    429    G4LorentzVector p4i =  theSystem->GetParticipant( i )->Get4Momentum();
    430    G4LorentzVector p4j =  theSystem->GetParticipant( j )->Get4Momentum();
     488//081103
     489   //G4cout << "CalFinalStateOfTheBinaryCollision " << i << " " << j << " " << theSystem->GetTotalNumberOfParticipant() << G4endl;
     490
     491   G4bool result = false;
     492   G4bool energyOK = false;
     493   G4bool pauliOK = false;
     494   G4bool abs = false;
     495   G4QMDParticipant* absorbed = NULL; 
     496
     497   G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum();
     498   G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum();
     499
     500//071031
     501
     502   G4double epot = theMeanField->GetTotalPotential();
     503
     504   G4double eini = epot + p4i.e() + p4j.e();
    431505
    432506//071031
    433507   // will use KineticTrack
    434    G4LorentzVector p4ix = p4i*GeV;
    435    G4LorentzVector p4jx = p4j*GeV;
    436    G4ThreeVector rix = (theSystem->GetParticipant( i )->GetPosition())*fermi;
    437    G4ThreeVector rjx = (theSystem->GetParticipant( j )->GetPosition())*fermi;
    438 //071031
    439 
    440    G4double epot = theMeanField->GetTotalPotential();
    441 
    442    G4double eini = epot + p4i.e() + p4j.e();
    443 
    444 
    445 //071031
    446508   G4ParticleDefinition* pdi0 =theSystem->GetParticipant( i )->GetDefinition();
    447509   G4ParticleDefinition* pdj0 =theSystem->GetParticipant( j )->GetDefinition();
    448    G4ThreeVector ri0 =(theSystem->GetParticipant( i )->GetPosition())*fermi;
    449    G4ThreeVector rj0 =(theSystem->GetParticipant( j )->GetPosition())*fermi;
     510   G4LorentzVector p4i0 = p4i*GeV;
     511   G4LorentzVector p4j0 = p4j*GeV;
     512   G4ThreeVector ri0 = ( theSystem->GetParticipant( i )->GetPosition() )*fermi;
     513   G4ThreeVector rj0 = ( theSystem->GetParticipant( j )->GetPosition() )*fermi;
    450514
    451515   for ( G4int iitry = 0 ; iitry < 4 ; iitry++ )
    452516   {
    453517
    454       G4KineticTrack kt1( pdi0 , 0.0 , ri0 , p4ix );
    455       G4KineticTrack kt2( pdj0 , 0.0 , rj0 , p4jx );
     518      abs = false;
     519
     520      G4KineticTrack kt1( pdi0 , 0.0 , ri0 , p4i0 );
     521      G4KineticTrack kt2( pdj0 , 0.0 , rj0 , p4j0 );
     522
    456523      G4LorentzVector p4ix_new;
    457524      G4LorentzVector p4jx_new;
    … …  
    462529      //std::cout << "G4QMDSCATTERER BEFORE " << kt2.GetDefinition()->GetParticleName() << " " << kt2.Get4Momentum()/GeV << " " << kt2.GetPosition()/fermi << std::endl;
    463530      //std::cout << "THESCATTERER " << theScatterer->GetCrossSection ( kt1 , kt2 )/millibarn << " " << elastic << " " << sig << std::endl;
     531
     532
    464533      if ( secs )
    465534      {
    … …  
    488557            }
    489558         }
    490          else
    491          {
    492             //std::cout << "NAGISA pion absrorption " << secs->front()->GetDefinition()->GetParticleName() << std::endl;
     559         else if ( secs->size() == 1 )
     560         {
     561//081118
     562            abs = true;
     563            //G4cout << "G4QMDCollision pion absrorption " << secs->front()->GetDefinition()->GetParticleName() << G4endl;
    493564            //secs->front()->Decay();
    494565            theSystem->GetParticipant( i )->SetDefinition( secs->front()->GetDefinition() );
    495566            p4ix_new = secs->front()->Get4Momentum()/GeV;
    496567            theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() );
    497            
    498               //std::cout << "THESCATTERER " << (*it)->GetDefinition()->GetParticleName() << std::endl;
    499               p4jx_new( 0 );
    500               //theSystem->GetParticipant( j )->SetDefinition( G4Gamma::Gamma() );
    501               //theSystem->GetParticipant( j )->SetDefinition( G4Neutron::Neutron() );
    502               theSystem->GetParticipant( j )->SetDefinition( G4PionZero::PionZero() );
    503               theSystem->GetParticipant( j )->SetMomentum( G4ThreeVector( G4UniformRand() )*eV );
    504               theSystem->GetParticipant( j )->SetPosition( G4ThreeVector( 1000, 1000, 1000 )*km );
    505568
    506569         }
    507570
    508          if ( secs->size() > 2 ) std::cout << "NAGISA secs size > 2;  " << secs->size() << std::endl;
     571//081118
     572         if ( secs->size() > 2 )
     573         {
     574
     575            G4cout << "G4QMDCollision secs size > 2;  " << secs->size() << G4endl;
     576
     577            for ( G4KineticTrackVector::iterator it
     578                = secs->begin() ; it != secs->end() ; it++ )
     579            {
     580               G4cout << "G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() << " " << (*it)->Get4Momentum()/GeV << G4endl;
     581            }
     582
     583         }
    509584
    510585         // deleteing KineticTrack
    … …  
    514589            delete *it;
    515590         }
     591
     592         delete secs;
    516593      }
    517594//071031
    518595
    519       theMeanField->Cal2BodyQuantities( i );
    520       theMeanField->Cal2BodyQuantities( j );
     596      if ( !abs )
     597      {
     598         theMeanField->Cal2BodyQuantities( i );
     599         theMeanField->Cal2BodyQuantities( j );
     600      }
     601      else
     602      {
     603         absorbed = theSystem->EraseParticipant( j );
     604         theMeanField->Update();
     605      }
    521606
    522607      epot = theMeanField->GetTotalPotential();
    … …  
    540625         //std::cout << "collisions after " << theSystem->GetParticipant( i )->Get4Momentum() << " " << theSystem->GetParticipant( j )->Get4Momentum() << std::endl;
    541626         //std::cout << "collisions dif " << ( p4i + p4j ) - ( theSystem->GetParticipant( i )->Get4Momentum() + theSystem->GetParticipant( j )->Get4Momentum() ) << std::endl;
    542          //std::cout << "collisions before " << rix/fermi << " " << rjx/fermi << std::endl;
     627         //std::cout << "collisions before " << ri0/fermi << " " << rj0/fermi << std::endl;
    543628         //std::cout << "collisions after " << theSystem->GetParticipant( i )->GetPosition() << " " << theSystem->GetParticipant( j )->GetPosition() << std::endl;
     629         energyOK = true;
     630         break;
     631      }
     632      else
     633      {
     634         //G4cout << "Energy Not OK " << G4endl;
     635         if ( abs )
     636         {
     637            //G4cout << "TKDB reinsert j " << G4endl;
     638            theSystem->InsertParticipant( absorbed , j );   
     639            theMeanField->Update();
     640         }
     641         // do not need reinsert in no absroption case
    544642      }
    545643//071031
    546 
    547       if ( std::abs ( eini - efin ) < epse ) return result;  // Collison OK
    548 
    549644   }
    550645
    551 //  Energetically forbidden collision
    552     result = false;
    553 
    554     return result;
     646// Energetically forbidden collision
     647
     648   if ( energyOK )
     649   {
     650      // Pauli Check
     651      //G4cout << "Pauli Checking " << theSystem->GetTotalNumberOfParticipant() << G4endl;
     652      if ( !abs )
     653      {
     654         if ( !( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) )
     655         {
     656            //G4cout << "Binary Collision Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl;
     657            pauliOK = true;
     658         }
     659      }
     660      else
     661      {
     662         if ( theMeanField->IsPauliBlocked ( i ) == false )
     663         {
     664            //G4cout << "Absorption Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl;
     665            delete absorbed;
     666            pauliOK = true;
     667         }
     668      }
     669     
     670
     671      if ( pauliOK )
     672      {
     673         result = true;
     674      }
     675      else
     676      {
     677         //G4cout << "Pauli Blocked" << G4endl;
     678         if ( abs )
     679         {
     680            //G4cout << "TKDB reinsert j pauli block" << G4endl;
     681            theSystem->InsertParticipant( absorbed , j );   
     682            theMeanField->Update();
     683         }
     684      }
     685   }
     686
     687   return result;
    555688
    556689}
Note: See TracChangeset for help on using the changeset viewer.