44 , fbcmax0 ( 1.323142 )
125 for (
G4int i = 0 ; i <
n ; i++ )
135 for (
G4int i = 0 ; i <
n ; i++ )
156 G4bool isThisEnergyOK =
false;
158 G4int maximumNumberOfTrial=4;
159 for (
G4int ii = 0 ; ii < maximumNumberOfTrial ; ii++ )
175 for ( G4KineticTrackVector::iterator it
176 = secs->begin() ; it != secs->end() ; it++ )
191 et += (*it)->Get4Momentum().e()/
GeV;
197 et += (*it)->Get4Momentum().e()/
GeV;
221 if ( std::abs ( eini - efin ) <
fepse*10 )
224 isThisEnergyOK =
true;
236 for (
G4int i0i =
id-2 ; 0 <= i0i ; i0i-- ) {
249 if ( isThisEnergyOK ==
true )
255 for (
G4int i0i = 0 ; i0i < i0 ; i0i++ )
278 for (
G4int i0i = 0 ; i0i < i0 ; i0i++ )
289 if ( isThisEnergyOK ==
true )
296 for (
G4int i0i = 0 ; i0i < i0 ; i0i++ )
332 for (
G4int j = 0 ; j < i ; j++ )
376 G4double srt = std::sqrt( (p4i+p4j)*(p4i+p4j) );
383 if ( rmi < 0.94 && rmj < 0.94 )
386 cutoff = rmi + rmj + 0.02;
400 if ( srt < cutoff )
continue;
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 ) );
414 if ( brel > fbcmax )
continue;
417 G4double bji = -pjdr/rmj + pidr * rmj /pij;
419 G4double ti = ( pidr/rmi - bij / aij ) * p4i.e() / rmi;
420 G4double tj = (-pjdr/rmj - bji / aij ) * p4j.e() / rmj;
437 if ( std::abs ( ti + tj ) > deltaT )
continue;
448 if ( prcm <= 0.00001 )
continue;
479 if ( energetically_forbidden ==
true )
502 G4bool absorption =
false;
578 G4double eini = epot + p4i.e() + p4j.e();
589 for (
G4int iitry = 0 ; iitry < 4 ; iitry++ )
610 if ( secs->size() == 2 )
612 for ( G4KineticTrackVector::iterator it
613 = secs->begin() ; it != secs->end() ; it++ )
618 p4ix_new = (*it)->Get4Momentum()/
GeV;
625 p4jx_new = (*it)->Get4Momentum()/
GeV;
633 else if ( secs->size() == 1 )
640 p4ix_new = secs->front()->Get4Momentum()/
GeV;
646 if ( secs->size() > 2 )
649 G4cout <<
"G4QMDCollision secs size > 2; " << secs->size() <<
G4endl;
651 for ( G4KineticTrackVector::iterator it
652 = secs->begin() ; it != secs->end() ; it++ )
654 G4cout <<
"G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() <<
" " << (*it)->Get4Momentum()/GeV <<
G4endl;
660 for ( G4KineticTrackVector::iterator it
661 = secs->begin() ; it != secs->end() ; it++ )
683 G4double efin = epot + p4ix_new.e() + p4jx_new.e();
694 if ( std::abs ( eini - efin ) <
fepse )
804 elastic = 35.0 / ( 1. + csrt * 100.0 ) + 20.0;
808 elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 )
809 * 2. /
pi + 1.0 ) * 9.65 + 7.0;
816 elastic = 28.0 / ( 1. + csrt * 100.0 ) + 27.0;
820 elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 )
821 * 2. /
pi + 1.0 ) * 12.34 + 10.0;
850 if( std::abs(c1) > 1.0 ) c1 = 2.0 * x - 1.0;
863 if ( pcm.x() == 0.0 && pcm.y() == 0 )
869 t2 = std::atan2( pcm.y() , pcm.x() );
873 G4double s1 = std::sqrt ( 1.0 - c1*c1 );
874 G4double s2 = std::sqrt ( 1.0 - c2*c2 );
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) );
892 G4double eini = epot + p4i.e() + p4j.e();
902 for (
G4int itry = 0 ; itry < 4 ; itry++ )
905 G4double eicm = std::sqrt ( rmi*rmi + pcm*pcm );
908 G4double trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + eicm );
912 G4double ejcm = std::sqrt ( rmj*rmj + pcm*pcm );
913 trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + ejcm );
936 G4double efin = epot + pi_new_e + pj_new_e ;
946 if ( std::abs ( eini - efin ) <
fepse )
958 if ( std::abs ( eini - efin ) <
fepse )
return result;
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 );
static G4Pow * GetInstance()
G4ThreeVector GetPosition()
static c2_factory< G4double > c2
G4QMDMeanField * theMeanField
G4double powN(G4double x, G4int n) const
void SetParticipant(G4QMDParticipant *particle)
CLHEP::Hep3Vector G4ThreeVector
G4KineticTrackVector * Decay()
void SetPosition(G4ThreeVector r)
G4int GetChargeInUnitOfEplus()
const G4ParticleDefinition * GetDefinition()
G4QMDParticipant * EraseParticipant(G4int i)
void Cal2BodyQuantities()
G4bool IsPauliBlocked(G4int)
G4ThreeVector GetMomentum()
void SetDefinition(const G4ParticleDefinition *pd)
G4GLOB_DLL std::ostream G4cout
G4Scatterer * theScatterer
void CalKinematicsOfBinaryCollisions(G4double)
G4int GetTotalNumberOfParticipant()
void SetMomentum(G4ThreeVector p)
G4double GetTotalPotential()
G4QMDParticipant * GetParticipant(G4int i)
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4bool IsShortLived() const
G4bool IsThisProjectile()
G4double GetPDGMass() const
void DeleteParticipant(G4int i)
G4LorentzVector Get4Momentum()
G4double GetRR2(G4int i, G4int j)
virtual G4KineticTrackVector * Scatter(const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
G4bool CalFinalStateOfTheBinaryCollisionJQMD(G4double, G4double, G4ThreeVector, G4double, G4double, G4ThreeVector, G4double, G4int, G4int)
const G4double x[NPOINTSGL]
G4bool CalFinalStateOfTheBinaryCollision(G4int, G4int)
void InsertParticipant(G4QMDParticipant *particle, G4int j)
void IncrementCollisionCounter()
static const double fermi
G4double elastic(Particle const *const p1, Particle const *const p2)
CLHEP::HepLorentzVector G4LorentzVector