#include <G4QMDCollision.hh>
Public Member Functions | |
G4QMDCollision () | |
~G4QMDCollision () | |
void | CalKinematicsOfBinaryCollisions (G4double) |
G4bool | CalFinalStateOfTheBinaryCollision (G4int, G4int) |
G4bool | CalFinalStateOfTheBinaryCollisionJQMD (G4double, G4double, G4ThreeVector, G4double, G4double, G4ThreeVector, G4double, G4int, G4int) |
void | SetMeanField (G4QMDMeanField *meanfield) |
Definition at line 47 of file G4QMDCollision.hh.
G4QMDCollision::G4QMDCollision | ( | ) |
Definition at line 39 of file G4QMDCollision.cc.
00040 : deltar ( 4 ) 00041 , bcmax0 ( 1.323142 ) // NN maximum impact parameter 00042 , bcmax1 ( 2.523 ) // others maximum impact parameter 00043 // , sig0 ( 55 ) // NN cross section 00044 //110617 fix for gcc 4.6 compilation warnings 00045 //, sig1 ( 200 ) // others cross section 00046 , epse ( 0.0001 ) 00047 { 00048 theScatterer = new G4Scatterer(); 00049 }
G4QMDCollision::~G4QMDCollision | ( | ) |
Definition at line 497 of file G4QMDCollision.cc.
References G4QMDMeanField::Cal2BodyQuantities(), G4QMDSystem::EraseParticipant(), G4cout, G4endl, G4QMDParticipant::Get4Momentum(), G4QMDParticipant::GetDefinition(), G4QMDSystem::GetParticipant(), G4QMDParticipant::GetPosition(), G4QMDMeanField::GetTotalPotential(), G4QMDSystem::InsertParticipant(), G4QMDMeanField::IsPauliBlocked(), G4Scatterer::Scatter(), G4QMDParticipant::SetDefinition(), G4QMDParticipant::SetMomentum(), and G4QMDMeanField::Update().
Referenced by CalKinematicsOfBinaryCollisions().
00498 { 00499 00500 //081103 00501 //G4cout << "CalFinalStateOfTheBinaryCollision " << i << " " << j << " " << theSystem->GetTotalNumberOfParticipant() << G4endl; 00502 00503 G4bool result = false; 00504 G4bool energyOK = false; 00505 G4bool pauliOK = false; 00506 G4bool abs = false; 00507 G4QMDParticipant* absorbed = NULL; 00508 00509 G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum(); 00510 G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum(); 00511 00512 //071031 00513 00514 G4double epot = theMeanField->GetTotalPotential(); 00515 00516 G4double eini = epot + p4i.e() + p4j.e(); 00517 00518 //071031 00519 // will use KineticTrack 00520 G4ParticleDefinition* pdi0 =theSystem->GetParticipant( i )->GetDefinition(); 00521 G4ParticleDefinition* pdj0 =theSystem->GetParticipant( j )->GetDefinition(); 00522 G4LorentzVector p4i0 = p4i*GeV; 00523 G4LorentzVector p4j0 = p4j*GeV; 00524 G4ThreeVector ri0 = ( theSystem->GetParticipant( i )->GetPosition() )*fermi; 00525 G4ThreeVector rj0 = ( theSystem->GetParticipant( j )->GetPosition() )*fermi; 00526 00527 for ( G4int iitry = 0 ; iitry < 4 ; iitry++ ) 00528 { 00529 00530 abs = false; 00531 00532 G4KineticTrack kt1( pdi0 , 0.0 , ri0 , p4i0 ); 00533 G4KineticTrack kt2( pdj0 , 0.0 , rj0 , p4j0 ); 00534 00535 G4LorentzVector p4ix_new; 00536 G4LorentzVector p4jx_new; 00537 G4KineticTrackVector* secs = NULL; 00538 secs = theScatterer->Scatter( kt1 , kt2 ); 00539 00540 //std::cout << "G4QMDSCATTERER BEFORE " << kt1.GetDefinition()->GetParticleName() << " " << kt1.Get4Momentum()/GeV << " " << kt1.GetPosition()/fermi << std::endl; 00541 //std::cout << "G4QMDSCATTERER BEFORE " << kt2.GetDefinition()->GetParticleName() << " " << kt2.Get4Momentum()/GeV << " " << kt2.GetPosition()/fermi << std::endl; 00542 //std::cout << "THESCATTERER " << theScatterer->GetCrossSection ( kt1 , kt2 )/millibarn << " " << elastic << " " << sig << std::endl; 00543 00544 00545 if ( secs ) 00546 { 00547 G4int iti = 0; 00548 if ( secs->size() == 2 ) 00549 { 00550 for ( G4KineticTrackVector::iterator it 00551 = secs->begin() ; it != secs->end() ; it++ ) 00552 { 00553 if ( iti == 0 ) 00554 { 00555 theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() ); 00556 p4ix_new = (*it)->Get4Momentum()/GeV; 00557 //std::cout << "THESCATTERER " << (*it)->GetDefinition()->GetParticleName() << std::endl; 00558 theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() ); 00559 } 00560 if ( iti == 1 ) 00561 { 00562 theSystem->GetParticipant( j )->SetDefinition( (*it)->GetDefinition() ); 00563 p4jx_new = (*it)->Get4Momentum()/GeV; 00564 //std::cout << "THESCATTERER " << p4jx_new.e()-p4jx_new.m() << std::endl; 00565 theSystem->GetParticipant( j )->SetMomentum( p4jx_new.v() ); 00566 } 00567 //std::cout << "G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() << " " << (*it)->Get4Momentum()/GeV << std::endl; 00568 iti++; 00569 } 00570 } 00571 else if ( secs->size() == 1 ) 00572 { 00573 //081118 00574 abs = true; 00575 //G4cout << "G4QMDCollision pion absrorption " << secs->front()->GetDefinition()->GetParticleName() << G4endl; 00576 //secs->front()->Decay(); 00577 theSystem->GetParticipant( i )->SetDefinition( secs->front()->GetDefinition() ); 00578 p4ix_new = secs->front()->Get4Momentum()/GeV; 00579 theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() ); 00580 00581 } 00582 00583 //081118 00584 if ( secs->size() > 2 ) 00585 { 00586 00587 G4cout << "G4QMDCollision secs size > 2; " << secs->size() << G4endl; 00588 00589 for ( G4KineticTrackVector::iterator it 00590 = secs->begin() ; it != secs->end() ; it++ ) 00591 { 00592 G4cout << "G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() << " " << (*it)->Get4Momentum()/GeV << G4endl; 00593 } 00594 00595 } 00596 00597 // deleteing KineticTrack 00598 for ( G4KineticTrackVector::iterator it 00599 = secs->begin() ; it != secs->end() ; it++ ) 00600 { 00601 delete *it; 00602 } 00603 00604 delete secs; 00605 } 00606 //071031 00607 00608 if ( !abs ) 00609 { 00610 theMeanField->Cal2BodyQuantities( i ); 00611 theMeanField->Cal2BodyQuantities( j ); 00612 } 00613 else 00614 { 00615 absorbed = theSystem->EraseParticipant( j ); 00616 theMeanField->Update(); 00617 } 00618 00619 epot = theMeanField->GetTotalPotential(); 00620 00621 G4double efin = epot + p4ix_new.e() + p4jx_new.e(); 00622 00623 //std::cout << "Collision NEW epot " << i << " " << j << " " << epot << " " << std::abs ( eini - efin ) - epse << std::endl; 00624 00625 /* 00626 std::cout << "Collision efin " << i << " " << j << " " << efin << std::endl; 00627 std::cout << "Collision " << i << " " << j << " " << std::abs ( eini - efin ) << " " << epse << std::endl; 00628 std::cout << "Collision " << std::abs ( eini - efin ) << " " << epse << std::endl; 00629 */ 00630 00631 //071031 00632 if ( std::abs ( eini - efin ) < epse ) 00633 { 00634 // Collison OK 00635 //std::cout << "collisions6" << std::endl; 00636 //std::cout << "collisions before " << p4i << " " << p4j << std::endl; 00637 //std::cout << "collisions after " << theSystem->GetParticipant( i )->Get4Momentum() << " " << theSystem->GetParticipant( j )->Get4Momentum() << std::endl; 00638 //std::cout << "collisions dif " << ( p4i + p4j ) - ( theSystem->GetParticipant( i )->Get4Momentum() + theSystem->GetParticipant( j )->Get4Momentum() ) << std::endl; 00639 //std::cout << "collisions before " << ri0/fermi << " " << rj0/fermi << std::endl; 00640 //std::cout << "collisions after " << theSystem->GetParticipant( i )->GetPosition() << " " << theSystem->GetParticipant( j )->GetPosition() << std::endl; 00641 energyOK = true; 00642 break; 00643 } 00644 else 00645 { 00646 //G4cout << "Energy Not OK " << G4endl; 00647 if ( abs ) 00648 { 00649 //G4cout << "TKDB reinsert j " << G4endl; 00650 theSystem->InsertParticipant( absorbed , j ); 00651 theMeanField->Update(); 00652 } 00653 // do not need reinsert in no absroption case 00654 } 00655 //071031 00656 } 00657 00658 // Energetically forbidden collision 00659 00660 if ( energyOK ) 00661 { 00662 // Pauli Check 00663 //G4cout << "Pauli Checking " << theSystem->GetTotalNumberOfParticipant() << G4endl; 00664 if ( !abs ) 00665 { 00666 if ( !( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) ) 00667 { 00668 //G4cout << "Binary Collision Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl; 00669 pauliOK = true; 00670 } 00671 } 00672 else 00673 { 00674 //if ( theMeanField->IsPauliBlocked ( i ) == false ) 00675 //090126 i-1 cause jth is erased 00676 if ( theMeanField->IsPauliBlocked ( i-1 ) == false ) 00677 { 00678 //G4cout << "Absorption Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl; 00679 delete absorbed; 00680 pauliOK = true; 00681 } 00682 } 00683 00684 00685 if ( pauliOK ) 00686 { 00687 result = true; 00688 } 00689 else 00690 { 00691 //G4cout << "Pauli Blocked" << G4endl; 00692 if ( abs ) 00693 { 00694 //G4cout << "TKDB reinsert j pauli block" << G4endl; 00695 theSystem->InsertParticipant( absorbed , j ); 00696 theMeanField->Update(); 00697 } 00698 } 00699 } 00700 00701 return result; 00702 00703 }
G4bool G4QMDCollision::CalFinalStateOfTheBinaryCollisionJQMD | ( | G4double | , | |
G4double | , | |||
G4ThreeVector | , | |||
G4double | , | |||
G4double | , | |||
G4ThreeVector | , | |||
G4double | , | |||
G4int | , | |||
G4int | ||||
) |
Definition at line 707 of file G4QMDCollision.cc.
References G4QMDMeanField::Cal2BodyQuantities(), G4UniformRand, G4QMDParticipant::Get4Momentum(), G4QMDParticipant::GetChargeInUnitOfEplus(), G4QMDParticipant::GetMass(), G4QMDSystem::GetParticipant(), G4QMDMeanField::GetTotalPotential(), G4INCL::Math::pi, and G4QMDParticipant::SetMomentum().
00708 { 00709 00710 //G4cout << "CalFinalStateOfTheBinaryCollisionJQMD" << G4endl; 00711 00712 G4bool result = true; 00713 00714 G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum(); 00715 G4double rmi = theSystem->GetParticipant( i )->GetMass(); 00716 G4int zi = theSystem->GetParticipant( i )->GetChargeInUnitOfEplus(); 00717 00718 G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum(); 00719 G4double rmj = theSystem->GetParticipant( j )->GetMass(); 00720 G4int zj = theSystem->GetParticipant( j )->GetChargeInUnitOfEplus(); 00721 00722 G4double pr = prcm; 00723 00724 G4double c2 = pcm.z()/pr; 00725 00726 G4double csrt = srt - cutoff; 00727 00728 //G4double pri = prcm; 00729 //G4double prf = sqrt ( 0.25 * srt*srt -rm2 ); 00730 00731 G4double asrt = srt - rmi - rmj; 00732 G4double pra = prcm; 00733 00734 00735 00736 G4double elastic = 0.0; 00737 00738 if ( zi == zj ) 00739 { 00740 if ( csrt < 0.4286 ) 00741 { 00742 elastic = 35.0 / ( 1. + csrt * 100.0 ) + 20.0; 00743 } 00744 else 00745 { 00746 elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 ) 00747 * 2. / pi + 1.0 ) * 9.65 + 7.0; 00748 } 00749 } 00750 else 00751 { 00752 if ( csrt < 0.4286 ) 00753 { 00754 elastic = 28.0 / ( 1. + csrt * 100.0 ) + 27.0; 00755 } 00756 else 00757 { 00758 elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 ) 00759 * 2. / pi + 1.0 ) * 12.34 + 10.0; 00760 } 00761 } 00762 00763 // std::cout << "Collision csrt " << i << " " << j << " " << csrt << std::endl; 00764 // std::cout << "Collision elstic " << i << " " << j << " " << elastic << std::endl; 00765 00766 00767 // std::cout << "Collision sig " << i << " " << j << " " << sig << std::endl; 00768 if ( G4UniformRand() > elastic / sig ) 00769 { 00770 //std::cout << "Inelastic " << std::endl; 00771 //std::cout << "elastic/sig " << elastic/sig << std::endl; 00772 return result; 00773 } 00774 else 00775 { 00776 //std::cout << "Elastic " << std::endl; 00777 } 00778 // std::cout << "Collision ELSTIC " << i << " " << j << std::endl; 00779 00780 00781 G4double as = std::pow ( 3.65 * asrt , 6 ); 00782 G4double a = 6.0 * as / (1.0 + as); 00783 G4double ta = -2.0 * pra*pra; 00784 G4double x = G4UniformRand(); 00785 G4double t1 = std::log( (1-x) * std::exp(2.*a*ta) + x ) / a; 00786 G4double c1 = 1.0 - t1/ta; 00787 00788 if( std::abs(c1) > 1.0 ) c1 = 2.0 * x - 1.0; 00789 00790 /* 00791 std::cout << "Collision as " << i << " " << j << " " << as << std::endl; 00792 std::cout << "Collision a " << i << " " << j << " " << a << std::endl; 00793 std::cout << "Collision ta " << i << " " << j << " " << ta << std::endl; 00794 std::cout << "Collision x " << i << " " << j << " " << x << std::endl; 00795 std::cout << "Collision t1 " << i << " " << j << " " << t1 << std::endl; 00796 std::cout << "Collision c1 " << i << " " << j << " " << c1 << std::endl; 00797 */ 00798 t1 = 2.0*pi*G4UniformRand(); 00799 // std::cout << "Collision t1 " << i << " " << j << " " << t1 << std::endl; 00800 G4double t2 = 0.0; 00801 if ( pcm.x() == 0.0 && pcm.y() == 0 ) 00802 { 00803 t2 = 0.0; 00804 } 00805 else 00806 { 00807 t2 = std::atan2( pcm.y() , pcm.x() ); 00808 } 00809 // std::cout << "Collision t2 " << i << " " << j << " " << t2 << std::endl; 00810 00811 G4double s1 = std::sqrt ( 1.0 - c1*c1 ); 00812 G4double s2 = std::sqrt ( 1.0 - c2*c2 ); 00813 00814 G4double ct1 = std::cos(t1); 00815 G4double st1 = std::sin(t1); 00816 00817 G4double ct2 = std::cos(t2); 00818 G4double st2 = std::sin(t2); 00819 00820 G4double ss = c2*s1*ct1 + s2*c1; 00821 00822 pcm.setX( pr * ( ss*ct2 - s1*st1*st2) ); 00823 pcm.setY( pr * ( ss*st2 + s1*st1*ct2) ); 00824 pcm.setZ( pr * ( c1*c2 - s1*s2*ct1) ); 00825 00826 // std::cout << "Collision pcm " << i << " " << j << " " << pcm << std::endl; 00827 00828 G4double epot = theMeanField->GetTotalPotential(); 00829 00830 G4double eini = epot + p4i.e() + p4j.e(); 00831 G4double etwo = p4i.e() + p4j.e(); 00832 00833 /* 00834 std::cout << "Collision epot " << i << " " << j << " " << epot << std::endl; 00835 std::cout << "Collision eini " << i << " " << j << " " << eini << std::endl; 00836 std::cout << "Collision etwo " << i << " " << j << " " << etwo << std::endl; 00837 */ 00838 00839 00840 for ( G4int itry = 0 ; itry < 4 ; itry++ ) 00841 { 00842 00843 G4double eicm = std::sqrt ( rmi*rmi + pcm*pcm ); 00844 G4double pibeta = pcm*beta; 00845 00846 G4double trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + eicm ); 00847 00848 G4ThreeVector pi_new = beta*trans + pcm; 00849 00850 G4double ejcm = std::sqrt ( rmj*rmj + pcm*pcm ); 00851 trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + ejcm ); 00852 00853 G4ThreeVector pj_new = beta*trans - pcm; 00854 00855 // 00856 // Delete old 00857 // Add new Particitipants 00858 // 00859 // Now only change momentum ( Beacuse we only have elastic sctter of nucleon 00860 // In future Definition also will be change 00861 // 00862 00863 theSystem->GetParticipant( i )->SetMomentum( pi_new ); 00864 theSystem->GetParticipant( j )->SetMomentum( pj_new ); 00865 00866 G4double pi_new_e = (theSystem->GetParticipant( i )->Get4Momentum()).e(); 00867 G4double pj_new_e = (theSystem->GetParticipant( j )->Get4Momentum()).e(); 00868 00869 theMeanField->Cal2BodyQuantities( i ); 00870 theMeanField->Cal2BodyQuantities( j ); 00871 00872 epot = theMeanField->GetTotalPotential(); 00873 00874 G4double efin = epot + pi_new_e + pj_new_e ; 00875 00876 //std::cout << "Collision NEW epot " << i << " " << j << " " << epot << " " << std::abs ( eini - efin ) - epse << std::endl; 00877 /* 00878 std::cout << "Collision efin " << i << " " << j << " " << efin << std::endl; 00879 std::cout << "Collision " << i << " " << j << " " << std::abs ( eini - efin ) << " " << epse << std::endl; 00880 std::cout << "Collision " << std::abs ( eini - efin ) << " " << epse << std::endl; 00881 */ 00882 00883 //071031 00884 if ( std::abs ( eini - efin ) < epse ) 00885 { 00886 // Collison OK 00887 //std::cout << "collisions6" << std::endl; 00888 //std::cout << "collisions before " << p4i << " " << p4j << std::endl; 00889 //std::cout << "collisions after " << theSystem->GetParticipant( i )->Get4Momentum() << " " << theSystem->GetParticipant( j )->Get4Momentum() << std::endl; 00890 //std::cout << "collisions dif " << ( p4i + p4j ) - ( theSystem->GetParticipant( i )->Get4Momentum() + theSystem->GetParticipant( j )->Get4Momentum() ) << std::endl; 00891 //std::cout << "collisions before " << rix/fermi << " " << rjx/fermi << std::endl; 00892 //std::cout << "collisions after " << theSystem->GetParticipant( i )->GetPosition() << " " << theSystem->GetParticipant( j )->GetPosition() << std::endl; 00893 } 00894 //071031 00895 00896 if ( std::abs ( eini - efin ) < epse ) return result; // Collison OK 00897 00898 G4double cona = ( eini - efin + etwo ) / gamma; 00899 G4double fac2 = 1.0 / ( 4.0 * cona*cona * pr*pr ) * 00900 ( ( cona*cona - ( rmi*rmi + rmj*rmj ) )*( cona*cona - ( rmi*rmi + rmj*rmj ) ) 00901 - 4.0 * rmi*rmi * rmj*rmj ); 00902 00903 if ( fac2 > 0 ) 00904 { 00905 G4double fact = std::sqrt ( fac2 ); 00906 pcm = fact*pcm; 00907 } 00908 00909 00910 } 00911 00912 // Energetically forbidden collision 00913 result = false; 00914 00915 return result; 00916 00917 }
void G4QMDCollision::CalKinematicsOfBinaryCollisions | ( | G4double | ) |
Definition at line 59 of file G4QMDCollision.cc.
References G4QMDMeanField::Cal2BodyQuantities(), CalFinalStateOfTheBinaryCollision(), G4KineticTrack::Decay(), G4QMDSystem::DeleteParticipant(), G4QMDParticipant::Get4Momentum(), G4QMDParticipant::GetDefinition(), G4QMDParticipant::GetMass(), G4QMDParticipant::GetMomentum(), G4QMDSystem::GetParticipant(), G4ParticleDefinition::GetPDGMass(), G4QMDParticipant::GetPosition(), G4QMDMeanField::GetRR2(), G4QMDSystem::GetTotalNumberOfParticipant(), G4QMDMeanField::GetTotalPotential(), G4QMDSystem::IncrementCollisionCounter(), G4QMDMeanField::IsPauliBlocked(), G4ParticleDefinition::IsShortLived(), G4QMDParticipant::IsThisHit(), G4QMDParticipant::IsThisProjectile(), G4QMDParticipant::IsThisTarget(), CLHEP::detail::n, G4QMDParticipant::SetDefinition(), G4QMDParticipant::SetHitMark(), G4QMDParticipant::SetMomentum(), G4QMDSystem::SetParticipant(), G4QMDParticipant::SetPosition(), G4QMDParticipant::UnsetHitMark(), G4QMDParticipant::UnsetInitialMark(), and G4QMDMeanField::Update().
Referenced by G4QMDReaction::ApplyYourself().
00060 { 00061 G4double deltaT = dt; 00062 00063 G4int n = theSystem->GetTotalNumberOfParticipant(); 00064 //081118 00065 //G4int nb = 0; 00066 for ( G4int i = 0 ; i < n ; i++ ) 00067 { 00068 theSystem->GetParticipant( i )->UnsetHitMark(); 00069 theSystem->GetParticipant( i )->UnsetHitMark(); 00070 //nb += theSystem->GetParticipant( i )->GetBaryonNumber(); 00071 } 00072 //G4cout << "nb = " << nb << " n = " << n << G4endl; 00073 00074 00075 //071101 00076 for ( G4int i = 0 ; i < n ; i++ ) 00077 { 00078 00079 //std::cout << i << " " << theSystem->GetParticipant( i )->GetDefinition()->GetParticleName() << " " << theSystem->GetParticipant( i )->GetPosition() << std::endl; 00080 00081 if ( theSystem->GetParticipant( i )->GetDefinition()->IsShortLived() ) 00082 { 00083 00084 G4bool decayed = false; 00085 00086 G4ParticleDefinition* pd0 = theSystem->GetParticipant( i )->GetDefinition(); 00087 G4ThreeVector p0 = theSystem->GetParticipant( i )->GetMomentum(); 00088 G4ThreeVector r0 = theSystem->GetParticipant( i )->GetPosition(); 00089 00090 G4LorentzVector p40 = theSystem->GetParticipant( i )->Get4Momentum(); 00091 00092 G4double eini = theMeanField->GetTotalPotential() + p40.e(); 00093 00094 G4int n0 = theSystem->GetTotalNumberOfParticipant(); 00095 G4int i0 = 0; 00096 00097 G4bool isThisEnergyOK = false; 00098 00099 for ( G4int ii = 0 ; ii < 4 ; ii++ ) 00100 { 00101 00102 //G4LorentzVector p4 = theSystem->GetParticipant( i )->Get4Momentum(); 00103 G4LorentzVector p400 = p40; 00104 00105 p400 *= GeV; 00106 //G4KineticTrack kt( theSystem->GetParticipant( i )->GetDefinition() , 0.0 , (theSystem->GetParticipant( i )->GetPosition())*fermi , p4 ); 00107 G4KineticTrack kt( pd0 , 0.0 , r0*fermi , p400 ); 00108 //std::cout << "G4KineticTrack " << i << " " << kt.GetDefinition()->GetParticleName() << kt.GetPosition() << std::endl; 00109 G4KineticTrackVector* secs = NULL; 00110 secs = kt.Decay(); 00111 G4int id = 0; 00112 G4double et = 0; 00113 if ( secs ) 00114 { 00115 for ( G4KineticTrackVector::iterator it 00116 = secs->begin() ; it != secs->end() ; it++ ) 00117 { 00118 /* 00119 G4cout << "G4KineticTrack" 00120 << " " << (*it)->GetDefinition()->GetParticleName() 00121 << " " << (*it)->Get4Momentum() 00122 << " " << (*it)->GetPosition()/fermi 00123 << G4endl; 00124 */ 00125 if ( id == 0 ) 00126 { 00127 theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() ); 00128 theSystem->GetParticipant( i )->SetMomentum( (*it)->Get4Momentum().v()/GeV ); 00129 theSystem->GetParticipant( i )->SetPosition( (*it)->GetPosition()/fermi ); 00130 //theMeanField->Cal2BodyQuantities( i ); 00131 et += (*it)->Get4Momentum().e()/GeV; 00132 } 00133 if ( id > 0 ) 00134 { 00135 // Append end; 00136 theSystem->SetParticipant ( new G4QMDParticipant ( (*it)->GetDefinition() , (*it)->Get4Momentum().v()/GeV , (*it)->GetPosition()/fermi ) ); 00137 et += (*it)->Get4Momentum().e()/GeV; 00138 if ( id > 1 ) 00139 { 00140 //081118 00141 //G4cout << "G4QMDCollision id >2; id= " << id << G4endl; 00142 } 00143 } 00144 id++; // number of daughter particles 00145 00146 delete *it; 00147 } 00148 00149 theMeanField->Update(); 00150 i0 = id-1; // 0 enter to i 00151 00152 delete secs; 00153 } 00154 00155 // EnergyCheck 00156 00157 G4double efin = theMeanField->GetTotalPotential() + et; 00158 //std::cout << std::abs ( eini - efin ) - epse << std::endl; 00159 // std::cout << std::abs ( eini - efin ) - epse*10 << std::endl; 00160 // *10 TK 00161 if ( std::abs ( eini - efin ) < epse*10 ) 00162 { 00163 // Energy OK 00164 isThisEnergyOK = true; 00165 break; 00166 } 00167 else 00168 { 00169 00170 theSystem->GetParticipant( i )->SetDefinition( pd0 ); 00171 theSystem->GetParticipant( i )->SetPosition( r0 ); 00172 theSystem->GetParticipant( i )->SetMomentum( p0 ); 00173 00174 for ( G4int i0i = 0 ; i0i < id-1 ; i0i++ ) 00175 { 00176 //081118 00177 //std::cout << "Decay Energitically Blocked deleteing " << i0i+n0 << std::endl; 00178 theSystem->DeleteParticipant( i0i+n0 ); 00179 } 00180 //081103 00181 theMeanField->Update(); 00182 } 00183 00184 } 00185 00186 00187 // Pauli Check 00188 if ( isThisEnergyOK == true ) 00189 { 00190 if ( theMeanField->IsPauliBlocked ( i ) != true ) 00191 { 00192 00193 G4bool allOK = true; 00194 for ( G4int i0i = 0 ; i0i < i0 ; i0i++ ) 00195 { 00196 if ( theMeanField->IsPauliBlocked ( i0i+n0 ) == true ) 00197 { 00198 allOK = false; 00199 break; 00200 } 00201 } 00202 00203 if ( allOK ) 00204 { 00205 decayed = true; //Decay Succeeded 00206 } 00207 } 00208 00209 } 00210 // 00211 00212 if ( decayed ) 00213 { 00214 //081119 00215 //G4cout << "Decay Suceeded! " << std::endl; 00216 theSystem->GetParticipant( i )->SetHitMark(); 00217 for ( G4int i0i = 0 ; i0i < i0 ; i0i++ ) 00218 { 00219 theSystem->GetParticipant( i0i+n0 )->SetHitMark(); 00220 } 00221 00222 } 00223 else 00224 { 00225 00226 // Decay Blocked and re-enter orginal participant; 00227 00228 if ( isThisEnergyOK == true ) // for false case already done 00229 { 00230 00231 theSystem->GetParticipant( i )->SetDefinition( pd0 ); 00232 theSystem->GetParticipant( i )->SetPosition( r0 ); 00233 theSystem->GetParticipant( i )->SetMomentum( p0 ); 00234 00235 for ( G4int i0i = 0 ; i0i < i0 ; i0i++ ) 00236 { 00237 //081118 00238 //std::cout << "Decay Blocked deleteing " << i0i+n0 << std::endl; 00239 theSystem->DeleteParticipant( i0i+n0 ); 00240 } 00241 //081103 00242 theMeanField->Update(); 00243 } 00244 00245 } 00246 00247 } //shortlive 00248 } // go next participant 00249 //071101 00250 00251 00252 n = theSystem->GetTotalNumberOfParticipant(); 00253 00254 //081118 00255 //for ( G4int i = 1 ; i < n ; i++ ) 00256 for ( G4int i = 1 ; i < theSystem->GetTotalNumberOfParticipant() ; i++ ) 00257 { 00258 00259 //std::cout << "Collision i " << i << std::endl; 00260 if ( theSystem->GetParticipant( i )->IsThisHit() ) continue; 00261 00262 G4ThreeVector ri = theSystem->GetParticipant( i )->GetPosition(); 00263 G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum(); 00264 G4double rmi = theSystem->GetParticipant( i )->GetMass(); 00265 G4ParticleDefinition* pdi = theSystem->GetParticipant( i )->GetDefinition(); 00266 //090331 gamma 00267 if ( pdi->GetPDGMass() == 0.0 ) continue; 00268 00269 //std::cout << " p4i00 " << p4i << std::endl; 00270 for ( G4int j = 0 ; j < i ; j++ ) 00271 { 00272 00273 00274 /* 00275 std::cout << "Collision " << i << " " << theSystem->GetParticipant( i )->IsThisProjectile() << std::endl; 00276 std::cout << "Collision " << j << " " << theSystem->GetParticipant( j )->IsThisProjectile() << std::endl; 00277 std::cout << "Collision " << i << " " << theSystem->GetParticipant( i )->IsThisTarget() << std::endl; 00278 std::cout << "Collision " << j << " " << theSystem->GetParticipant( j )->IsThisTarget() << std::endl; 00279 */ 00280 00281 // Only 1 Collision allowed for each particle in a time step. 00282 //081119 00283 if ( theSystem->GetParticipant( i )->IsThisHit() ) continue; 00284 if ( theSystem->GetParticipant( j )->IsThisHit() ) continue; 00285 00286 //std::cout << "Collision " << i << " " << j << std::endl; 00287 00288 // Do not allow collision between nucleons in target/projectile til its first collision. 00289 if ( theSystem->GetParticipant( i )->IsThisProjectile() ) 00290 { 00291 if ( theSystem->GetParticipant( j )->IsThisProjectile() ) continue; 00292 } 00293 else if ( theSystem->GetParticipant( i )->IsThisTarget() ) 00294 { 00295 if ( theSystem->GetParticipant( j )->IsThisTarget() ) continue; 00296 } 00297 00298 00299 G4ThreeVector rj = theSystem->GetParticipant( j )->GetPosition(); 00300 G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum(); 00301 G4double rmj = theSystem->GetParticipant( j )->GetMass(); 00302 G4ParticleDefinition* pdj = theSystem->GetParticipant( j )->GetDefinition(); 00303 //090331 gamma 00304 if ( pdj->GetPDGMass() == 0.0 ) continue; 00305 00306 G4double rr2 = theMeanField->GetRR2( i , j ); 00307 00308 // Here we assume elab (beam momentum less than 5 GeV/n ) 00309 if ( rr2 > deltar*deltar ) continue; 00310 00311 //G4double s = (p4i+p4j)*(p4i+p4j); 00312 //G4double srt = std::sqrt ( s ); 00313 00314 G4double srt = std::sqrt( (p4i+p4j)*(p4i+p4j) ); 00315 00316 G4double cutoff = 0.0; 00317 G4double bcmax = 0.0; 00318 //110617 fix for gcc 4.6 compilation warnings 00319 //G4double sig = 0.0; 00320 00321 if ( rmi < 0.94 && rmj < 0.94 ) 00322 { 00323 // nucleon or pion case 00324 cutoff = rmi + rmj + 0.02; 00325 bcmax = bcmax0; 00326 //110617 fix for gcc 4.6 compilation warnings 00327 //sig = sig0; 00328 } 00329 else 00330 { 00331 cutoff = rmi + rmj; 00332 bcmax = bcmax1; 00333 //110617 fix for gcc compilation warnings 00334 //sig = sig1; 00335 } 00336 00337 //std::cout << "Collision cutoff " << i << " " << j << " " << cutoff << std::endl; 00338 if ( srt < cutoff ) continue; 00339 00340 G4ThreeVector dr = ri - rj; 00341 G4double rsq = dr*dr; 00342 00343 G4double pij = p4i*p4j; 00344 G4double pidr = p4i.vect()*dr; 00345 G4double pjdr = p4j.vect()*dr; 00346 00347 G4double aij = 1.0 - ( rmi*rmj /pij ) * ( rmi*rmj /pij ); 00348 G4double bij = pidr / rmi - pjdr*rmi/pij; 00349 G4double cij = rsq + ( pidr / rmi ) * ( pidr / rmi ); 00350 G4double brel = std::sqrt ( std::abs ( cij - bij*bij/aij ) ); 00351 00352 if ( brel > bcmax ) continue; 00353 //std::cout << "collisions3 " << std::endl; 00354 00355 G4double bji = -pjdr/rmj + pidr * rmj /pij; 00356 00357 G4double ti = ( pidr/rmi - bij / aij ) * p4i.e() / rmi; 00358 G4double tj = (-pjdr/rmj - bji / aij ) * p4j.e() / rmj; 00359 00360 00361 /* 00362 std::cout << "collisions4 p4i " << p4i << std::endl; 00363 std::cout << "collisions4 ri " << ri << std::endl; 00364 std::cout << "collisions4 p4j " << p4j << std::endl; 00365 std::cout << "collisions4 rj " << rj << std::endl; 00366 std::cout << "collisions4 dr " << dr << std::endl; 00367 std::cout << "collisions4 pij " << pij << std::endl; 00368 std::cout << "collisions4 aij " << aij << std::endl; 00369 std::cout << "collisions4 bij bji " << bij << " " << bji << std::endl; 00370 std::cout << "collisions4 pidr pjdr " << pidr << " " << pjdr << std::endl; 00371 std::cout << "collisions4 p4i.e() p4j.e() " << p4i.e() << " " << p4j.e() << std::endl; 00372 std::cout << "collisions4 rmi rmj " << rmi << " " << rmj << std::endl; 00373 std::cout << "collisions4 " << ti << " " << tj << std::endl; 00374 */ 00375 if ( std::abs ( ti + tj ) > deltaT ) continue; 00376 //std::cout << "collisions4 " << std::endl; 00377 00378 G4ThreeVector beta = ( p4i + p4j ).boostVector(); 00379 00380 G4LorentzVector p = p4i; 00381 G4LorentzVector p4icm = p.boost( p.findBoostToCM ( p4j ) ); 00382 G4ThreeVector pcm = p4icm.vect(); 00383 00384 G4double prcm = pcm.mag(); 00385 00386 if ( prcm <= 0.00001 ) continue; 00387 //std::cout << "collisions5 " << std::endl; 00388 00389 G4bool energetically_forbidden = !( CalFinalStateOfTheBinaryCollision ( i , j ) ); // Use Geant4 Collision Library 00390 //G4bool energetically_forbidden = !( CalFinalStateOfTheBinaryCollisionJQMD ( sig , cutoff , pcm , prcm , srt, beta , gamma , i , j ) ); // JQMD Elastic 00391 00392 /* 00393 G4bool pauli_blocked = false; 00394 if ( energetically_forbidden == false ) // result true 00395 { 00396 if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) 00397 { 00398 pauli_blocked = true; 00399 //std::cout << "G4QMDRESULT Collsion Pauli Blocked " << std::endl; 00400 } 00401 } 00402 else 00403 { 00404 if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) 00405 pauli_blocked = false; 00406 //std::cout << "G4QMDRESULT Collsion Blocked " << std::endl; 00407 } 00408 */ 00409 00410 /* 00411 std::cout << "G4QMDRESULT Collsion initial p4 i and j " 00412 << p4i << " " << p4j 00413 << std::endl; 00414 */ 00415 // 081118 00416 //if ( energetically_forbidden == true || pauli_blocked == true ) 00417 if ( energetically_forbidden == true ) 00418 { 00419 00420 //G4cout << " energetically_forbidden " << G4endl; 00421 // Collsion not allowed then re enter orginal participants 00422 // Now only momentum, becasuse we only consider elastic scattering of nucleons 00423 00424 theSystem->GetParticipant( i )->SetMomentum( p4i.vect() ); 00425 theSystem->GetParticipant( i )->SetDefinition( pdi ); 00426 theSystem->GetParticipant( i )->SetPosition( ri ); 00427 00428 theSystem->GetParticipant( j )->SetMomentum( p4j.vect() ); 00429 theSystem->GetParticipant( j )->SetDefinition( pdj ); 00430 theSystem->GetParticipant( j )->SetPosition( rj ); 00431 00432 theMeanField->Cal2BodyQuantities( i ); 00433 theMeanField->Cal2BodyQuantities( j ); 00434 00435 } 00436 else 00437 { 00438 00439 00440 G4bool absorption = false; 00441 if ( n == theSystem->GetTotalNumberOfParticipant()+1 ) absorption = true; 00442 if ( absorption ) 00443 { 00444 //G4cout << "Absorption happend " << G4endl; 00445 i = i-1; 00446 n = n-1; 00447 } 00448 00449 // Collsion allowed (really happened) 00450 00451 // Unset Projectile/Target flag 00452 theSystem->GetParticipant( i )->UnsetInitialMark(); 00453 if ( !absorption ) theSystem->GetParticipant( j )->UnsetInitialMark(); 00454 00455 theSystem->GetParticipant( i )->SetHitMark(); 00456 if ( !absorption ) theSystem->GetParticipant( j )->SetHitMark(); 00457 00458 theSystem->IncrementCollisionCounter(); 00459 00460 /* 00461 std::cout << "G4QMDRESULT Collsion Really Happened between " 00462 << i << " and " << j 00463 << std::endl; 00464 std::cout << "G4QMDRESULT Collsion initial p4 i and j " 00465 << p4i << " " << p4j 00466 << std::endl; 00467 std::cout << "G4QMDRESULT Collsion after p4 i and j " 00468 << theSystem->GetParticipant( i )->Get4Momentum() 00469 << " " 00470 << theSystem->GetParticipant( j )->Get4Momentum() 00471 << std::endl; 00472 std::cout << "G4QMDRESULT Collsion Diff " 00473 << p4i + p4j - theSystem->GetParticipant( i )->Get4Momentum() - theSystem->GetParticipant( j )->Get4Momentum() 00474 << std::endl; 00475 std::cout << "G4QMDRESULT Collsion initial r i and j " 00476 << ri << " " << rj 00477 << std::endl; 00478 std::cout << "G4QMDRESULT Collsion after r i and j " 00479 << theSystem->GetParticipant( i )->GetPosition() 00480 << " " 00481 << theSystem->GetParticipant( j )->GetPosition() 00482 << std::endl; 00483 */ 00484 00485 00486 } 00487 00488 } 00489 00490 } 00491 00492 00493 }
void G4QMDCollision::SetMeanField | ( | G4QMDMeanField * | meanfield | ) | [inline] |
Definition at line 58 of file G4QMDCollision.hh.
References G4QMDMeanField::GetSystem().
Referenced by G4QMDReaction::ApplyYourself().
00058 { theMeanField = meanfield; theSystem = meanfield->GetSystem(); }