Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4QMDCollision Class Reference

#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)
 
void deltar (G4double x)
 
G4double deltar ()
 
void bcmax0 (G4double x)
 
G4double bcmax0 ()
 
void bcmax1 (G4double x)
 
G4double bcmax1 ()
 
void epse (G4double x)
 
G4double epse ()
 

Detailed Description

Definition at line 47 of file G4QMDCollision.hh.

Constructor & Destructor Documentation

◆ G4QMDCollision()

G4QMDCollision::G4QMDCollision ( )

Definition at line 42 of file G4QMDCollision.cc.

43: fdeltar ( 4.0 )
44, fbcmax0 ( 1.323142 ) // NN maximum impact parameter
45, fbcmax1 ( 2.523 ) // others maximum impact parameter
46// , sig0 ( 55 ) // NN cross section
47//110617 fix for gcc 4.6 compilation warnings
48//, sig1 ( 200 ) // others cross section
49, fepse ( 0.0001 )
50{
51 //These two pointers will be set through SetMeanField method
52 theSystem=NULL;
53 theMeanField=NULL;
54 theScatterer = new G4Scatterer();
55}

◆ ~G4QMDCollision()

G4QMDCollision::~G4QMDCollision ( )

Definition at line 110 of file G4QMDCollision.cc.

111{
112 //if ( theSystem != NULL ) delete theSystem;
113 //if ( theMeanField != NULL ) delete theMeanField;
114 delete theScatterer;
115}

Member Function Documentation

◆ bcmax0() [1/2]

G4double G4QMDCollision::bcmax0 ( )
inline

Definition at line 64 of file G4QMDCollision.hh.

64{ return fbcmax0; };

◆ bcmax0() [2/2]

void G4QMDCollision::bcmax0 ( G4double x)
inline

Definition at line 63 of file G4QMDCollision.hh.

63{ fbcmax0 = x; };

◆ bcmax1() [1/2]

G4double G4QMDCollision::bcmax1 ( )
inline

Definition at line 66 of file G4QMDCollision.hh.

66{ return fbcmax1; };

◆ bcmax1() [2/2]

void G4QMDCollision::bcmax1 ( G4double x)
inline

Definition at line 65 of file G4QMDCollision.hh.

65{ fbcmax1 = x; };

◆ CalFinalStateOfTheBinaryCollision()

G4bool G4QMDCollision::CalFinalStateOfTheBinaryCollision ( G4int i,
G4int j )

Definition at line 559 of file G4QMDCollision.cc.

560{
561
562//081103
563 //G4cout << "CalFinalStateOfTheBinaryCollision " << i << " " << j << " " << theSystem->GetTotalNumberOfParticipant() << G4endl;
564
565 G4bool result = false;
566 G4bool energyOK = false;
567 G4bool pauliOK = false;
568 G4bool abs = false;
569 G4QMDParticipant* absorbed = NULL;
570
571 G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum();
572 G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum();
573
574//071031
575
576 G4double epot = theMeanField->GetTotalPotential();
577
578 G4double eini = epot + p4i.e() + p4j.e();
579
580//071031
581 // will use KineticTrack
582 const G4ParticleDefinition* pdi0 =theSystem->GetParticipant( i )->GetDefinition();
583 const G4ParticleDefinition* pdj0 =theSystem->GetParticipant( j )->GetDefinition();
584 G4LorentzVector p4i0 = p4i*GeV;
585 G4LorentzVector p4j0 = p4j*GeV;
586 G4ThreeVector ri0 = ( theSystem->GetParticipant( i )->GetPosition() )*fermi;
587 G4ThreeVector rj0 = ( theSystem->GetParticipant( j )->GetPosition() )*fermi;
588
589 for ( G4int iitry = 0 ; iitry < 4 ; iitry++ )
590 {
591
592 abs = false;
593
594 G4KineticTrack kt1( pdi0 , 0.0 , ri0 , p4i0 );
595 G4KineticTrack kt2( pdj0 , 0.0 , rj0 , p4j0 );
596
597 G4LorentzVector p4ix_new;
598 G4LorentzVector p4jx_new;
599 G4KineticTrackVector* secs = NULL;
600 secs = theScatterer->Scatter( kt1 , kt2 );
601
602 //std::cout << "G4QMDSCATTERER BEFORE " << kt1.GetDefinition()->GetParticleName() << " " << kt1.Get4Momentum()/GeV << " " << kt1.GetPosition()/fermi << std::endl;
603 //std::cout << "G4QMDSCATTERER BEFORE " << kt2.GetDefinition()->GetParticleName() << " " << kt2.Get4Momentum()/GeV << " " << kt2.GetPosition()/fermi << std::endl;
604 //std::cout << "THESCATTERER " << theScatterer->GetCrossSection ( kt1 , kt2 )/millibarn << " " << elastic << " " << sig << std::endl;
605
606
607 if ( secs )
608 {
609 G4int iti = 0;
610 if ( secs->size() == 2 )
611 {
612 for ( G4KineticTrackVector::iterator it
613 = secs->begin() ; it != secs->end() ; it++ )
614 {
615 if ( iti == 0 )
616 {
617 theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() );
618 p4ix_new = (*it)->Get4Momentum()/GeV;
619 //std::cout << "THESCATTERER " << (*it)->GetDefinition()->GetParticleName() << std::endl;
620 theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() );
621 }
622 if ( iti == 1 )
623 {
624 theSystem->GetParticipant( j )->SetDefinition( (*it)->GetDefinition() );
625 p4jx_new = (*it)->Get4Momentum()/GeV;
626 //std::cout << "THESCATTERER " << p4jx_new.e()-p4jx_new.m() << std::endl;
627 theSystem->GetParticipant( j )->SetMomentum( p4jx_new.v() );
628 }
629 //std::cout << "G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() << " " << (*it)->Get4Momentum()/GeV << std::endl;
630 iti++;
631 }
632 }
633 else if ( secs->size() == 1 )
634 {
635//081118
636 abs = true;
637 //G4cout << "G4QMDCollision pion absrorption " << secs->front()->GetDefinition()->GetParticleName() << G4endl;
638 //secs->front()->Decay();
639 theSystem->GetParticipant( i )->SetDefinition( secs->front()->GetDefinition() );
640 p4ix_new = secs->front()->Get4Momentum()/GeV;
641 theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() );
642
643 }
644
645//081118
646 if ( secs->size() > 2 )
647 {
648
649 G4cout << "G4QMDCollision secs size > 2; " << secs->size() << G4endl;
650
651 for ( G4KineticTrackVector::iterator it
652 = secs->begin() ; it != secs->end() ; it++ )
653 {
654 G4cout << "G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() << " " << (*it)->Get4Momentum()/GeV << G4endl;
655 }
656
657 }
658
659 // deleteing KineticTrack
660 for ( G4KineticTrackVector::iterator it
661 = secs->begin() ; it != secs->end() ; it++ )
662 {
663 delete *it;
664 }
665
666 delete secs;
667 }
668//071031
669
670 if ( !abs )
671 {
672 theMeanField->Cal2BodyQuantities( i );
673 theMeanField->Cal2BodyQuantities( j );
674 }
675 else
676 {
677 absorbed = theSystem->EraseParticipant( j );
678 theMeanField->Update();
679 }
680
681 epot = theMeanField->GetTotalPotential();
682
683 G4double efin = epot + p4ix_new.e() + p4jx_new.e();
684
685 //std::cout << "Collision NEW epot " << i << " " << j << " " << epot << " " << std::abs ( eini - efin ) - fepse << std::endl;
686
687/*
688 G4cout << "Collision efin " << i << " " << j << " " << efin << G4endl;
689 G4cout << "Collision " << i << " " << j << " " << std::abs ( eini - efin ) << " " << fepse << G4endl;
690 G4cout << "Collision " << std::abs ( eini - efin ) << " " << fepse << G4endl;
691*/
692
693//071031
694 if ( std::abs ( eini - efin ) < fepse )
695 {
696 // Collison OK
697 //std::cout << "collisions6" << std::endl;
698 //std::cout << "collisions before " << p4i << " " << p4j << std::endl;
699 //std::cout << "collisions after " << theSystem->GetParticipant( i )->Get4Momentum() << " " << theSystem->GetParticipant( j )->Get4Momentum() << std::endl;
700 //std::cout << "collisions dif " << ( p4i + p4j ) - ( theSystem->GetParticipant( i )->Get4Momentum() + theSystem->GetParticipant( j )->Get4Momentum() ) << std::endl;
701 //std::cout << "collisions before " << ri0/fermi << " " << rj0/fermi << std::endl;
702 //std::cout << "collisions after " << theSystem->GetParticipant( i )->GetPosition() << " " << theSystem->GetParticipant( j )->GetPosition() << std::endl;
703 energyOK = true;
704 break;
705 }
706 else
707 {
708 //G4cout << "Energy Not OK " << G4endl;
709 if ( abs )
710 {
711 //G4cout << "TKDB reinsert j " << G4endl;
712 theSystem->InsertParticipant( absorbed , j );
713 theMeanField->Update();
714 }
715 // do not need reinsert in no absroption case
716 }
717//071031
718 }
719
720// Energetically forbidden collision
721
722 if ( energyOK )
723 {
724 // Pauli Check
725 //G4cout << "Pauli Checking " << theSystem->GetTotalNumberOfParticipant() << G4endl;
726 if ( !abs )
727 {
728 if ( !( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) )
729 {
730 //G4cout << "Binary Collision Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl;
731 pauliOK = true;
732 }
733 }
734 else
735 {
736 //if ( theMeanField->IsPauliBlocked ( i ) == false )
737 //090126 i-1 cause jth is erased
738 if ( theMeanField->IsPauliBlocked ( i-1 ) == false )
739 {
740 //G4cout << "Absorption Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl;
741 delete absorbed;
742 pauliOK = true;
743 }
744 }
745
746
747 if ( pauliOK )
748 {
749 result = true;
750 }
751 else
752 {
753 //G4cout << "Pauli Blocked" << G4endl;
754 if ( abs )
755 {
756 //G4cout << "TKDB reinsert j pauli block" << G4endl;
757 theSystem->InsertParticipant( absorbed , j );
758 theMeanField->Update();
759 }
760 }
761 }
762
763 return result;
764
765}
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
CLHEP::HepLorentzVector::v
Hep3Vector v() const
G4QMDMeanField::GetTotalPotential
G4double GetTotalPotential()
Definition G4QMDMeanField.cc:405
G4QMDMeanField::IsPauliBlocked
G4bool IsPauliBlocked(G4int)
Definition G4QMDMeanField.cc:469
G4QMDParticipant::GetDefinition
const G4ParticleDefinition * GetDefinition()
Definition G4QMDParticipant.hh:56
G4QMDParticipant::Get4Momentum
G4LorentzVector Get4Momentum()
Definition G4QMDParticipant.cc:50
G4QMDParticipant::SetDefinition
void SetDefinition(const G4ParticleDefinition *pd)
Definition G4QMDParticipant.hh:55
G4QMDParticipant::SetMomentum
void SetMomentum(G4ThreeVector p)
Definition G4QMDParticipant.hh:61
G4QMDSystem::InsertParticipant
void InsertParticipant(G4QMDParticipant *particle, G4int j)
Definition G4QMDSystem.cc:110
G4QMDSystem::GetParticipant
G4QMDParticipant * GetParticipant(G4int i)
Definition G4QMDSystem.hh:83
G4QMDSystem::EraseParticipant
G4QMDParticipant * EraseParticipant(G4int i)
Definition G4QMDSystem.hh:61
G4Scatterer::Scatter
virtual G4KineticTrackVector * Scatter(const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
Definition G4Scatterer.cc:283

Referenced by CalKinematicsOfBinaryCollisions().

◆ CalFinalStateOfTheBinaryCollisionJQMD()

G4bool G4QMDCollision::CalFinalStateOfTheBinaryCollisionJQMD ( G4double sig,
G4double cutoff,
G4ThreeVector pcm,
G4double prcm,
G4double srt,
G4ThreeVector beta,
G4double gamma,
G4int i,
G4int j )

Definition at line 769 of file G4QMDCollision.cc.

770{
771
772 //G4cout << "CalFinalStateOfTheBinaryCollisionJQMD" << G4endl;
773
774 G4bool result = true;
775
776 G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum();
777 G4double rmi = theSystem->GetParticipant( i )->GetMass();
778 G4int zi = theSystem->GetParticipant( i )->GetChargeInUnitOfEplus();
779
780 G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum();
781 G4double rmj = theSystem->GetParticipant( j )->GetMass();
782 G4int zj = theSystem->GetParticipant( j )->GetChargeInUnitOfEplus();
783
784 G4double pr = prcm;
785
786 G4double c2 = pcm.z()/pr;
787
788 G4double csrt = srt - cutoff;
789
790 //G4double pri = prcm;
791 //G4double prf = sqrt ( 0.25 * srt*srt -rm2 );
792
793 G4double asrt = srt - rmi - rmj;
794 G4double pra = prcm;
795
796
797
798 G4double elastic = 0.0;
799
800 if ( zi == zj )
801 {
802 if ( csrt < 0.4286 )
803 {
804 elastic = 35.0 / ( 1. + csrt * 100.0 ) + 20.0;
805 }
806 else
807 {
808 elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 )
809 * 2. / pi + 1.0 ) * 9.65 + 7.0;
810 }
811 }
812 else
813 {
814 if ( csrt < 0.4286 )
815 {
816 elastic = 28.0 / ( 1. + csrt * 100.0 ) + 27.0;
817 }
818 else
819 {
820 elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 )
821 * 2. / pi + 1.0 ) * 12.34 + 10.0;
822 }
823 }
824
825// std::cout << "Collision csrt " << i << " " << j << " " << csrt << std::endl;
826// std::cout << "Collision elstic " << i << " " << j << " " << elastic << std::endl;
827
828
829// std::cout << "Collision sig " << i << " " << j << " " << sig << std::endl;
830 if ( G4UniformRand() > elastic / sig )
831 {
832 //std::cout << "Inelastic " << std::endl;
833 //std::cout << "elastic/sig " << elastic/sig << std::endl;
834 return result;
835 }
836 else
837 {
838 //std::cout << "Elastic " << std::endl;
839 }
840// std::cout << "Collision ELSTIC " << i << " " << j << std::endl;
841
842
843 G4double as = G4Pow::GetInstance()->powN ( 3.65 * asrt , 6 );
844 G4double a = 6.0 * as / (1.0 + as);
845 G4double ta = -2.0 * pra*pra;
846 G4double x = G4UniformRand();
847 G4double t1 = G4Log( (1-x) * G4Exp(2.*a*ta) + x ) / a;
848 G4double c1 = 1.0 - t1/ta;
849
850 if( std::abs(c1) > 1.0 ) c1 = 2.0 * x - 1.0;
851
852/*
853 G4cout << "Collision as " << i << " " << j << " " << as << G4endl;
854 G4cout << "Collision a " << i << " " << j << " " << a << G4endl;
855 G4cout << "Collision ta " << i << " " << j << " " << ta << G4endl;
856 G4cout << "Collision x " << i << " " << j << " " << x << G4endl;
857 G4cout << "Collision t1 " << i << " " << j << " " << t1 << G4endl;
858 G4cout << "Collision c1 " << i << " " << j << " " << c1 << G4endl;
859*/
860 t1 = 2.0*pi*G4UniformRand();
861// std::cout << "Collision t1 " << i << " " << j << " " << t1 << std::endl;
862 G4double t2 = 0.0;
863 if ( pcm.x() == 0.0 && pcm.y() == 0 )
864 {
865 t2 = 0.0;
866 }
867 else
868 {
869 t2 = std::atan2( pcm.y() , pcm.x() );
870 }
871// std::cout << "Collision t2 " << i << " " << j << " " << t2 << std::endl;
872
873 G4double s1 = std::sqrt ( 1.0 - c1*c1 );
874 G4double s2 = std::sqrt ( 1.0 - c2*c2 );
875
876 G4double ct1 = std::cos(t1);
877 G4double st1 = std::sin(t1);
878
879 G4double ct2 = std::cos(t2);
880 G4double st2 = std::sin(t2);
881
882 G4double ss = c2*s1*ct1 + s2*c1;
883
884 pcm.setX( pr * ( ss*ct2 - s1*st1*st2) );
885 pcm.setY( pr * ( ss*st2 + s1*st1*ct2) );
886 pcm.setZ( pr * ( c1*c2 - s1*s2*ct1) );
887
888// std::cout << "Collision pcm " << i << " " << j << " " << pcm << std::endl;
889
890 G4double epot = theMeanField->GetTotalPotential();
891
892 G4double eini = epot + p4i.e() + p4j.e();
893 G4double etwo = p4i.e() + p4j.e();
894
895/*
896 G4cout << "Collision epot " << i << " " << j << " " << epot << G4endl;
897 G4cout << "Collision eini " << i << " " << j << " " << eini << G4endl;
898 G4cout << "Collision etwo " << i << " " << j << " " << etwo << G4endl;
899*/
900
901
902 for ( G4int itry = 0 ; itry < 4 ; itry++ )
903 {
904
905 G4double eicm = std::sqrt ( rmi*rmi + pcm*pcm );
906 G4double pibeta = pcm*beta;
907
908 G4double trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + eicm );
909
910 G4ThreeVector pi_new = beta*trans + pcm;
911
912 G4double ejcm = std::sqrt ( rmj*rmj + pcm*pcm );
913 trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + ejcm );
914
915 G4ThreeVector pj_new = beta*trans - pcm;
916
917//
918// Delete old
919// Add new Particitipants
920//
921// Now only change momentum ( Beacuse we only have elastic sctter of nucleon
922// In future Definition also will be change
923//
924
925 theSystem->GetParticipant( i )->SetMomentum( pi_new );
926 theSystem->GetParticipant( j )->SetMomentum( pj_new );
927
928 G4double pi_new_e = (theSystem->GetParticipant( i )->Get4Momentum()).e();
929 G4double pj_new_e = (theSystem->GetParticipant( j )->Get4Momentum()).e();
930
931 theMeanField->Cal2BodyQuantities( i );
932 theMeanField->Cal2BodyQuantities( j );
933
934 epot = theMeanField->GetTotalPotential();
935
936 G4double efin = epot + pi_new_e + pj_new_e ;
937
938 //std::cout << "Collision NEW epot " << i << " " << j << " " << epot << " " << std::abs ( eini - efin ) - fepse << std::endl;
939/*
940 G4cout << "Collision efin " << i << " " << j << " " << efin << G4endl;
941 G4cout << "Collision " << i << " " << j << " " << std::abs ( eini - efin ) << " " << fepse << G4endl;
942 G4cout << "Collision " << std::abs ( eini - efin ) << " " << fepse << G4endl;
943*/
944
945//071031
946 if ( std::abs ( eini - efin ) < fepse )
947 {
948 // Collison OK
949 //std::cout << "collisions6" << std::endl;
950 //std::cout << "collisions before " << p4i << " " << p4j << std::endl;
951 //std::cout << "collisions after " << theSystem->GetParticipant( i )->Get4Momentum() << " " << theSystem->GetParticipant( j )->Get4Momentum() << std::endl;
952 //std::cout << "collisions dif " << ( p4i + p4j ) - ( theSystem->GetParticipant( i )->Get4Momentum() + theSystem->GetParticipant( j )->Get4Momentum() ) << std::endl;
953 //std::cout << "collisions before " << rix/fermi << " " << rjx/fermi << std::endl;
954 //std::cout << "collisions after " << theSystem->GetParticipant( i )->GetPosition() << " " << theSystem->GetParticipant( j )->GetPosition() << std::endl;
955 }
956//071031
957
958 if ( std::abs ( eini - efin ) < fepse ) return result; // Collison OK
959
960 G4double cona = ( eini - efin + etwo ) / gamma;
961 G4double fac2 = 1.0 / ( 4.0 * cona*cona * pr*pr ) *
962 ( ( cona*cona - ( rmi*rmi + rmj*rmj ) )*( cona*cona - ( rmi*rmi + rmj*rmj ) )
963 - 4.0 * rmi*rmi * rmj*rmj );
964
965 if ( fac2 > 0 )
966 {
967 G4double fact = std::sqrt ( fac2 );
968 pcm = fact*pcm;
969 }
970
971
972 }
973
974// Energetically forbidden collision
975 result = false;
976
977 return result;
978
979}
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:180
G4Log
G4double G4Log(G4double x)
Definition G4Log.hh:227
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::setY
void setY(double)
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::setZ
void setZ(double)
CLHEP::Hep3Vector::setX
void setX(double)
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition G4Pow.cc:41
G4Pow::powN
G4double powN(G4double x, G4int n) const
Definition G4Pow.cc:162
G4INCL::CrossSections::elastic
G4double elastic(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:57

◆ CalKinematicsOfBinaryCollisions()

void G4QMDCollision::CalKinematicsOfBinaryCollisions ( G4double dt)

Definition at line 118 of file G4QMDCollision.cc.

119{
120 G4double deltaT = dt;
121
122 G4int n = theSystem->GetTotalNumberOfParticipant();
123//081118
124 //G4int nb = 0;
125 for ( G4int i = 0 ; i < n ; i++ )
126 {
127 theSystem->GetParticipant( i )->UnsetHitMark();
128 theSystem->GetParticipant( i )->UnsetHitMark();
129 //nb += theSystem->GetParticipant( i )->GetBaryonNumber();
130 }
131 //G4cout << "nb = " << nb << " n = " << n << G4endl;
132
133
134//071101
135 for ( G4int i = 0 ; i < n ; i++ )
136 {
137
138 //std::cout << i << " " << theSystem->GetParticipant( i )->GetDefinition()->GetParticleName() << " " << theSystem->GetParticipant( i )->GetPosition() << std::endl;
139
140 if ( theSystem->GetParticipant( i )->GetDefinition()->IsShortLived() )
141 {
142
143 G4bool decayed = false;
144
145 const G4ParticleDefinition* pd0 = theSystem->GetParticipant( i )->GetDefinition();
146 G4ThreeVector p0 = theSystem->GetParticipant( i )->GetMomentum();
147 G4ThreeVector r0 = theSystem->GetParticipant( i )->GetPosition();
148
149 G4LorentzVector p40 = theSystem->GetParticipant( i )->Get4Momentum();
150
151 G4double eini = theMeanField->GetTotalPotential() + p40.e();
152
153 G4int n0 = theSystem->GetTotalNumberOfParticipant();
154 G4int i0 = 0;
155
156 G4bool isThisEnergyOK = false;
157
158 G4int maximumNumberOfTrial=4;
159 for ( G4int ii = 0 ; ii < maximumNumberOfTrial ; ii++ )
160 {
161
162 //G4LorentzVector p4 = theSystem->GetParticipant( i )->Get4Momentum();
163 G4LorentzVector p400 = p40;
164
165 p400 *= GeV;
166 //G4KineticTrack kt( theSystem->GetParticipant( i )->GetDefinition() , 0.0 , (theSystem->GetParticipant( i )->GetPosition())*fermi , p4 );
167 G4KineticTrack kt( pd0 , 0.0 , r0*fermi , p400 );
168 //std::cout << "G4KineticTrack " << i << " " << kt.GetDefinition()->GetParticleName() << kt.GetPosition() << std::endl;
169 G4KineticTrackVector* secs = NULL;
170 secs = kt.Decay();
171 G4int id = 0;
172 G4double et = 0;
173 if ( secs )
174 {
175 for ( G4KineticTrackVector::iterator it
176 = secs->begin() ; it != secs->end() ; it++ )
177 {
178/*
179 G4cout << "G4KineticTrack"
180 << " " << (*it)->GetDefinition()->GetParticleName()
181 << " " << (*it)->Get4Momentum()
182 << " " << (*it)->GetPosition()/fermi
183 << G4endl;
184*/
185 if ( id == 0 )
186 {
187 theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() );
188 theSystem->GetParticipant( i )->SetMomentum( (*it)->Get4Momentum().v()/GeV );
189 theSystem->GetParticipant( i )->SetPosition( (*it)->GetPosition()/fermi );
190 //theMeanField->Cal2BodyQuantities( i );
191 et += (*it)->Get4Momentum().e()/GeV;
192 }
193 if ( id > 0 )
194 {
195 // Append end;
196 theSystem->SetParticipant ( new G4QMDParticipant ( (*it)->GetDefinition() , (*it)->Get4Momentum().v()/GeV , (*it)->GetPosition()/fermi ) );
197 et += (*it)->Get4Momentum().e()/GeV;
198 if ( id > 1 )
199 {
200 //081118
201 //G4cout << "G4QMDCollision id >2; id= " << id << G4endl;
202 }
203 }
204 id++; // number of daughter particles
205
206 delete *it;
207 }
208
209 theMeanField->Update();
210 i0 = id-1; // 0 enter to i
211
212 delete secs;
213 }
214
215// EnergyCheck
216
217 G4double efin = theMeanField->GetTotalPotential() + et;
218 //std::cout << std::abs ( eini - efin ) - fepse << std::endl;
219// std::cout << std::abs ( eini - efin ) - fepse*10 << std::endl;
220// *10 TK
221 if ( std::abs ( eini - efin ) < fepse*10 )
222 {
223 // Energy OK
224 isThisEnergyOK = true;
225 break;
226 }
227 else
228 {
229
230 theSystem->GetParticipant( i )->SetDefinition( pd0 );
231 theSystem->GetParticipant( i )->SetPosition( r0 );
232 theSystem->GetParticipant( i )->SetMomentum( p0 );
233
234 //for ( G4int i0i = 0 ; i0i < id-1 ; i0i++ )
235 //160210 deletion must be done in descending order
236 for ( G4int i0i = id-2 ; 0 <= i0i ; i0i-- ) {
237 //081118
238 //std::cout << "Decay Energitically Blocked deleteing " << i0i+n0 << std::endl;
239 theSystem->DeleteParticipant( i0i+n0 );
240 }
241 //081103
242 theMeanField->Update();
243 }
244
245 }
246
247
248// Pauli Check
249 if ( isThisEnergyOK == true )
250 {
251 if ( theMeanField->IsPauliBlocked ( i ) != true )
252 {
253
254 G4bool allOK = true;
255 for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
256 {
257 if ( theMeanField->IsPauliBlocked ( i0i+n0 ) == true )
258 {
259 allOK = false;
260 break;
261 }
262 }
263
264 if ( allOK )
265 {
266 decayed = true; //Decay Succeeded
267 }
268 }
269
270 }
271//
272
273 if ( decayed )
274 {
275 //081119
276 //G4cout << "Decay Suceeded! " << std::endl;
277 theSystem->GetParticipant( i )->SetHitMark();
278 for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
279 {
280 theSystem->GetParticipant( i0i+n0 )->SetHitMark();
281 }
282
283 }
284 else
285 {
286
287// Decay Blocked and re-enter orginal participant;
288
289 if ( isThisEnergyOK == true ) // for false case already done
290 {
291
292 theSystem->GetParticipant( i )->SetDefinition( pd0 );
293 theSystem->GetParticipant( i )->SetPosition( r0 );
294 theSystem->GetParticipant( i )->SetMomentum( p0 );
295
296 for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
297 {
298 //081118
299 //std::cout << "Decay Blocked deleteing " << i0i+n0 << std::endl;
300 //160210 adding commnet: deletion must be done in descending order
301 theSystem->DeleteParticipant( i0+n0-i0i-1 );
302 }
303 //081103
304 theMeanField->Update();
305 }
306
307 }
308
309 } //shortlive
310 } // go next participant
311//071101
312
313
314 n = theSystem->GetTotalNumberOfParticipant();
315
316//081118
317 //for ( G4int i = 1 ; i < n ; i++ )
318 for ( G4int i = 1 ; i < theSystem->GetTotalNumberOfParticipant() ; i++ )
319 {
320
321 //std::cout << "Collision i " << i << std::endl;
322 if ( theSystem->GetParticipant( i )->IsThisHit() ) continue;
323
324 G4ThreeVector ri = theSystem->GetParticipant( i )->GetPosition();
325 G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum();
326 G4double rmi = theSystem->GetParticipant( i )->GetMass();
327 const G4ParticleDefinition* pdi = theSystem->GetParticipant( i )->GetDefinition();
328//090331 gamma
329 if ( pdi->GetPDGMass() == 0.0 ) continue;
330
331 //std::cout << " p4i00 " << p4i << std::endl;
332 for ( G4int j = 0 ; j < i ; j++ )
333 {
334
335
336/*
337 G4cout << "Collision " << i << " " << theSystem->GetParticipant( i )->IsThisProjectile() << G4endl;
338 G4cout << "Collision " << j << " " << theSystem->GetParticipant( j )->IsThisProjectile() << G4endl;
339 G4cout << "Collision " << i << " " << theSystem->GetParticipant( i )->IsThisTarget() << G4endl;
340 G4cout << "Collision " << j << " " << theSystem->GetParticipant( j )->IsThisTarget() << G4endl;
341*/
342
343 // Only 1 Collision allowed for each particle in a time step.
344 //081119
345 if ( theSystem->GetParticipant( i )->IsThisHit() ) continue;
346 if ( theSystem->GetParticipant( j )->IsThisHit() ) continue;
347
348 //std::cout << "Collision " << i << " " << j << std::endl;
349
350 // Do not allow collision between nucleons in target/projectile til its first collision.
351 if ( theSystem->GetParticipant( i )->IsThisProjectile() )
352 {
353 if ( theSystem->GetParticipant( j )->IsThisProjectile() ) continue;
354 }
355 else if ( theSystem->GetParticipant( i )->IsThisTarget() )
356 {
357 if ( theSystem->GetParticipant( j )->IsThisTarget() ) continue;
358 }
359
360
361 G4ThreeVector rj = theSystem->GetParticipant( j )->GetPosition();
362 G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum();
363 G4double rmj = theSystem->GetParticipant( j )->GetMass();
364 const G4ParticleDefinition* pdj = theSystem->GetParticipant( j )->GetDefinition();
365//090331 gamma
366 if ( pdj->GetPDGMass() == 0.0 ) continue;
367
368 G4double rr2 = theMeanField->GetRR2( i , j );
369
370// Here we assume elab (beam momentum less than 5 GeV/n )
371 if ( rr2 > fdeltar*fdeltar ) continue;
372
373 //G4double s = (p4i+p4j)*(p4i+p4j);
374 //G4double srt = std::sqrt ( s );
375
376 G4double srt = std::sqrt( (p4i+p4j)*(p4i+p4j) );
377
378 G4double cutoff = 0.0;
379 G4double fbcmax = 0.0;
380 //110617 fix for gcc 4.6 compilation warnings
381 //G4double sig = 0.0;
382
383 if ( rmi < 0.94 && rmj < 0.94 )
384 {
385// nucleon or pion case
386 cutoff = rmi + rmj + 0.02;
387 fbcmax = fbcmax0;
388 //110617 fix for gcc 4.6 compilation warnings
389 //sig = sig0;
390 }
391 else
392 {
393 cutoff = rmi + rmj;
394 fbcmax = fbcmax1;
395 //110617 fix for gcc compilation warnings
396 //sig = sig1;
397 }
398
399 //std::cout << "Collision cutoff " << i << " " << j << " " << cutoff << std::endl;
400 if ( srt < cutoff ) continue;
401
402 G4ThreeVector dr = ri - rj;
403 G4double rsq = dr*dr;
404
405 G4double pij = p4i*p4j;
406 G4double pidr = p4i.vect()*dr;
407 G4double pjdr = p4j.vect()*dr;
408
409 G4double aij = 1.0 - ( rmi*rmj /pij ) * ( rmi*rmj /pij );
410 G4double bij = pidr / rmi - pjdr*rmi/pij;
411 G4double cij = rsq + ( pidr / rmi ) * ( pidr / rmi );
412 G4double brel = std::sqrt ( std::abs ( cij - bij*bij/aij ) );
413
414 if ( brel > fbcmax ) continue;
415 //std::cout << "collisions3 " << std::endl;
416
417 G4double bji = -pjdr/rmj + pidr * rmj /pij;
418
419 G4double ti = ( pidr/rmi - bij / aij ) * p4i.e() / rmi;
420 G4double tj = (-pjdr/rmj - bji / aij ) * p4j.e() / rmj;
421
422
423/*
424 G4cout << "collisions4 p4i " << p4i << G4endl;
425 G4cout << "collisions4 ri " << ri << G4endl;
426 G4cout << "collisions4 p4j " << p4j << G4endl;
427 G4cout << "collisions4 rj " << rj << G4endl;
428 G4cout << "collisions4 dr " << dr << G4endl;
429 G4cout << "collisions4 pij " << pij << G4endl;
430 G4cout << "collisions4 aij " << aij << G4endl;
431 G4cout << "collisions4 bij bji " << bij << " " << bji << G4endl;
432 G4cout << "collisions4 pidr pjdr " << pidr << " " << pjdr << G4endl;
433 G4cout << "collisions4 p4i.e() p4j.e() " << p4i.e() << " " << p4j.e() << G4endl;
434 G4cout << "collisions4 rmi rmj " << rmi << " " << rmj << G4endl;
435 G4cout << "collisions4 " << ti << " " << tj << G4endl;
436*/
437 if ( std::abs ( ti + tj ) > deltaT ) continue;
438 //std::cout << "collisions4 " << std::endl;
439
440 G4ThreeVector beta = ( p4i + p4j ).boostVector();
441
442 G4LorentzVector p = p4i;
443 G4LorentzVector p4icm = p.boost( p.findBoostToCM ( p4j ) );
444 G4ThreeVector pcm = p4icm.vect();
445
446 G4double prcm = pcm.mag();
447
448 if ( prcm <= 0.00001 ) continue;
449 //std::cout << "collisions5 " << std::endl;
450
451 G4bool energetically_forbidden = !( CalFinalStateOfTheBinaryCollision ( i , j ) ); // Use Geant4 Collision Library
452 //G4bool energetically_forbidden = !( CalFinalStateOfTheBinaryCollisionJQMD ( sig , cutoff , pcm , prcm , srt, beta , gamma , i , j ) ); // JQMD Elastic
453
454/*
455 G4bool pauli_blocked = false;
456 if ( energetically_forbidden == false ) // result true
457 {
458 if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true )
459 {
460 pauli_blocked = true;
461 //std::cout << "G4QMDRESULT Collsion Pauli Blocked " << std::endl;
462 }
463 }
464 else
465 {
466 if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true )
467 pauli_blocked = false;
468 //std::cout << "G4QMDRESULT Collsion Blocked " << std::endl;
469 }
470*/
471
472/*
473 G4cout << "G4QMDRESULT Collsion initial p4 i and j "
474 << p4i << " " << p4j
475 << G4endl;
476*/
477// 081118
478 //if ( energetically_forbidden == true || pauli_blocked == true )
479 if ( energetically_forbidden == true )
480 {
481
482 //G4cout << " energetically_forbidden " << G4endl;
483// Collsion not allowed then re enter orginal participants
484// Now only momentum, becasuse we only consider elastic scattering of nucleons
485
486 theSystem->GetParticipant( i )->SetMomentum( p4i.vect() );
487 theSystem->GetParticipant( i )->SetDefinition( pdi );
488 theSystem->GetParticipant( i )->SetPosition( ri );
489
490 theSystem->GetParticipant( j )->SetMomentum( p4j.vect() );
491 theSystem->GetParticipant( j )->SetDefinition( pdj );
492 theSystem->GetParticipant( j )->SetPosition( rj );
493
494 theMeanField->Cal2BodyQuantities( i );
495 theMeanField->Cal2BodyQuantities( j );
496
497 }
498 else
499 {
500
501
502 G4bool absorption = false;
503 if ( n == theSystem->GetTotalNumberOfParticipant()+1 ) absorption = true;
504 if ( absorption )
505 {
506 //G4cout << "Absorption happend " << G4endl;
507 i = i-1;
508 n = n-1;
509 }
510
511// Collsion allowed (really happened)
512
513 // Unset Projectile/Target flag
514 theSystem->GetParticipant( i )->UnsetInitialMark();
515 if ( !absorption ) theSystem->GetParticipant( j )->UnsetInitialMark();
516
517 theSystem->GetParticipant( i )->SetHitMark();
518 if ( !absorption ) theSystem->GetParticipant( j )->SetHitMark();
519
520 theSystem->IncrementCollisionCounter();
521
522/*
523 G4cout << "G4QMDRESULT Collsion Really Happened between "
524 << i << " and " << j
525 << G4endl;
526 G4cout << "G4QMDRESULT Collsion initial p4 i and j "
527 << p4i << " " << p4j
528 << G4endl;
529 G4cout << "G4QMDRESULT Collsion after p4 i and j "
530 << theSystem->GetParticipant( i )->Get4Momentum()
531 << " "
532 << theSystem->GetParticipant( j )->Get4Momentum()
533 << G4endl;
534 G4cout << "G4QMDRESULT Collsion Diff "
535 << p4i + p4j - theSystem->GetParticipant( i )->Get4Momentum() - theSystem->GetParticipant( j )->Get4Momentum()
536 << G4endl;
537 G4cout << "G4QMDRESULT Collsion initial r i and j "
538 << ri << " " << rj
539 << G4endl;
540 G4cout << "G4QMDRESULT Collsion after r i and j "
541 << theSystem->GetParticipant( i )->GetPosition()
542 << " "
543 << theSystem->GetParticipant( j )->GetPosition()
544 << G4endl;
545*/
546
547
548 }
549
550 }
551
552 }
553
554
555}
CLHEP::Hep3Vector::mag
double mag() const
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition LorentzVector.cc:55
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
CLHEP::HepLorentzVector::findBoostToCM
Hep3Vector findBoostToCM() const
Definition LorentzVectorK.cc:220
G4QMDCollision::CalFinalStateOfTheBinaryCollision
G4bool CalFinalStateOfTheBinaryCollision(G4int, G4int)
Definition G4QMDCollision.cc:559
G4QMDMeanField::GetRR2
G4double GetRR2(G4int i, G4int j)
Definition G4QMDMeanField.hh:61
G4QMDParticipant::GetPosition
G4ThreeVector GetPosition()
Definition G4QMDParticipant.hh:59
G4QMDParticipant::SetPosition
void SetPosition(G4ThreeVector r)
Definition G4QMDParticipant.hh:58
G4QMDParticipant::GetMomentum
G4ThreeVector GetMomentum()
Definition G4QMDParticipant.hh:62
G4QMDSystem::GetTotalNumberOfParticipant
G4int GetTotalNumberOfParticipant()
Definition G4QMDSystem.hh:78
G4QMDSystem::DeleteParticipant
void DeleteParticipant(G4int i)
Definition G4QMDSystem.hh:69
G4QMDSystem::SetParticipant
void SetParticipant(G4QMDParticipant *particle)
Definition G4QMDSystem.hh:52
G4QMDSystem::IncrementCollisionCounter
void IncrementCollisionCounter()
Definition G4QMDSystem.hh:88

Referenced by G4QMDReaction::ApplyYourself().

◆ deltar() [1/2]

G4double G4QMDCollision::deltar ( )
inline

Definition at line 62 of file G4QMDCollision.hh.

62{ return fdeltar; };

◆ deltar() [2/2]

void G4QMDCollision::deltar ( G4double x)
inline

Definition at line 61 of file G4QMDCollision.hh.

61{ fdeltar = x; };

◆ epse() [1/2]

G4double G4QMDCollision::epse ( )
inline

Definition at line 68 of file G4QMDCollision.hh.

68{ return fepse; };

◆ epse() [2/2]

void G4QMDCollision::epse ( G4double x)
inline

Definition at line 67 of file G4QMDCollision.hh.

67{ fepse = x; };

◆ SetMeanField()

void G4QMDCollision::SetMeanField ( G4QMDMeanField * meanfield)
inline

Definition at line 58 of file G4QMDCollision.hh.

58{ theMeanField = meanfield; theSystem = meanfield->GetSystem(); }
G4QMDMeanField::GetSystem
G4QMDSystem * GetSystem()
Definition G4QMDMeanField.hh:45

Referenced by G4QMDReaction::ApplyYourself().

The documentation for this class was generated from the following files: