Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
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 ( )

Definition at line 39 of file G4QMDCollision.cc.

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

Definition at line 107 of file G4QMDCollision.cc.

108 {
109  //if ( theSystem != NULL ) delete theSystem;
110  //if ( theMeanField != NULL ) delete theMeanField;
111  delete theScatterer;
112 }

Member Function Documentation

void G4QMDCollision::bcmax0 ( G4double  x)
inline

Definition at line 65 of file G4QMDCollision.hh.

References test::x.

65 { fbcmax0 = x; };
G4double G4QMDCollision::bcmax0 ( )
inline

Definition at line 66 of file G4QMDCollision.hh.

66 { return fbcmax0; };
void G4QMDCollision::bcmax1 ( G4double  x)
inline

Definition at line 67 of file G4QMDCollision.hh.

References test::x.

67 { fbcmax1 = x; };
G4double G4QMDCollision::bcmax1 ( )
inline

Definition at line 68 of file G4QMDCollision.hh.

68 { return fbcmax1; };
G4bool G4QMDCollision::CalFinalStateOfTheBinaryCollision ( G4int  i,
G4int  j 
)

Definition at line 553 of file G4QMDCollision.cc.

References absorbed, G4QMDMeanField::Cal2BodyQuantities(), CLHEP::HepLorentzVector::e(), G4QMDSystem::EraseParticipant(), python.hepunit::fermi, G4cout, G4endl, G4QMDParticipant::Get4Momentum(), G4QMDParticipant::GetDefinition(), G4QMDSystem::GetParticipant(), G4QMDParticipant::GetPosition(), G4QMDMeanField::GetTotalPotential(), python.hepunit::GeV, G4QMDSystem::InsertParticipant(), G4QMDMeanField::IsPauliBlocked(), G4Scatterer::Scatter(), G4QMDParticipant::SetDefinition(), G4QMDParticipant::SetMomentum(), G4QMDMeanField::Update(), and CLHEP::HepLorentzVector::v().

Referenced by CalKinematicsOfBinaryCollisions().

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

Definition at line 763 of file G4QMDCollision.cc.

References test::a, plottest35::c1, G4QMDMeanField::Cal2BodyQuantities(), CLHEP::HepLorentzVector::e(), G4INCL::CrossSections::elastic(), G4UniformRand, G4QMDParticipant::Get4Momentum(), G4QMDParticipant::GetChargeInUnitOfEplus(), G4QMDParticipant::GetMass(), G4QMDSystem::GetParticipant(), G4QMDMeanField::GetTotalPotential(), python.hepunit::pi, gammaraytel::pr, G4QMDParticipant::SetMomentum(), CLHEP::Hep3Vector::setX(), CLHEP::Hep3Vector::setY(), CLHEP::Hep3Vector::setZ(), plottest35::t1, test::x, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

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

Definition at line 115 of file G4QMDCollision.cc.

References CLHEP::HepLorentzVector::boost(), G4QMDMeanField::Cal2BodyQuantities(), CalFinalStateOfTheBinaryCollision(), G4KineticTrack::Decay(), G4QMDSystem::DeleteParticipant(), CLHEP::HepLorentzVector::e(), python.hepunit::fermi, CLHEP::HepLorentzVector::findBoostToCM(), G4QMDParticipant::Get4Momentum(), G4QMDParticipant::GetDefinition(), G4QMDParticipant::GetMass(), G4QMDParticipant::GetMomentum(), G4QMDSystem::GetParticipant(), G4ParticleDefinition::GetPDGMass(), G4QMDParticipant::GetPosition(), G4QMDMeanField::GetRR2(), G4QMDSystem::GetTotalNumberOfParticipant(), G4QMDMeanField::GetTotalPotential(), python.hepunit::GeV, G4QMDSystem::IncrementCollisionCounter(), G4QMDMeanField::IsPauliBlocked(), G4ParticleDefinition::IsShortLived(), G4QMDParticipant::IsThisHit(), G4QMDParticipant::IsThisProjectile(), G4QMDParticipant::IsThisTarget(), CLHEP::Hep3Vector::mag(), n, G4QMDParticipant::SetDefinition(), G4QMDParticipant::SetHitMark(), G4QMDParticipant::SetMomentum(), G4QMDSystem::SetParticipant(), G4QMDParticipant::SetPosition(), G4QMDParticipant::UnsetHitMark(), G4QMDParticipant::UnsetInitialMark(), G4QMDMeanField::Update(), and CLHEP::HepLorentzVector::vect().

Referenced by G4QMDReaction::ApplyYourself().

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

Definition at line 63 of file G4QMDCollision.hh.

References test::x.

63 { fdeltar = x; };
G4double G4QMDCollision::deltar ( )
inline

Definition at line 64 of file G4QMDCollision.hh.

64 { return fdeltar; };
void G4QMDCollision::epse ( G4double  x)
inline

Definition at line 69 of file G4QMDCollision.hh.

References test::x.

69 { fepse = x; };
G4double G4QMDCollision::epse ( )
inline

Definition at line 70 of file G4QMDCollision.hh.

70 { return fepse; };
void G4QMDCollision::SetMeanField ( G4QMDMeanField *  meanfield)
inline

Definition at line 60 of file G4QMDCollision.hh.

References G4QMDMeanField::GetSystem().

Referenced by G4QMDReaction::ApplyYourself().

60 { theMeanField = meanfield; theSystem = meanfield->GetSystem(); }
G4QMDMeanField::GetSystem
G4QMDSystem * GetSystem()
Definition: G4QMDMeanField.hh:52
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:59 for Geant4.10 by   doxygen 1.8.7