153 std::vector< G4QMDNucleus* > nucleuses;
162 G4double e1 = std::sqrt( p1*p1 + m1*m1 );
168 G4double beta_nncm = ( - boostLABtoCM.beta() + boostLABtoNN.
beta() ) / ( 1 - boostLABtoCM.beta() * boostLABtoNN.
beta() ) ;
176 boostToReac = boostLABtoNN;
177 boostBackToLAB = -boostLABtoNN;
182 G4int icounter_max = 1024;
186 if ( icounter > icounter_max ) {
187 G4cout <<
"Loop-counter exceeded the threshold value at " << __LINE__ <<
"th line of " << __FILE__ <<
"." <<
G4endl;
240 targ->CalEnergyAndAngularMomentumInCM();
247 for (
G4int i = 0 ; i < targ->GetTotalNumberOfParticipant() ; i++ )
303 G4int i = targ->GetTotalNumberOfParticipant();
340 if ( i / 10 * 10 == i )
365 if ( numberOfSecondary == 2 )
368 G4bool elasticLike_system =
false;
369 if ( nucleuses.size() == 2 )
373 sec_a_A = nucleuses[0]->GetMassNumber();
374 sec_b_Z = nucleuses[1]->GetAtomicNumber();
375 sec_b_A = nucleuses[1]->GetMassNumber();
377 if ( ( sec_a_Z == proj_Z && sec_a_A == proj_A && sec_b_Z == targ_Z && sec_b_A == targ_A )
378 || ( sec_a_Z == targ_Z && sec_a_A == targ_A && sec_b_Z == proj_Z && sec_b_A == proj_A ) )
380 elasticLike_system =
true;
384 else if ( nucleuses.size() == 1 )
387 sec_a_Z = nucleuses[0]->GetAtomicNumber();
388 sec_a_A = nucleuses[0]->GetMassNumber();
391 if ( ( sec_a_Z == proj_Z && sec_a_A == proj_A && sec_b_pd == targ_pd )
392 || ( sec_a_Z == targ_Z && sec_a_A == targ_A && sec_b_pd == proj_pd ) )
394 elasticLike_system =
true;
404 if ( ( sec_a_pd == proj_pd && sec_b_pd == targ_pd )
405 || ( sec_a_pd == targ_pd && sec_b_pd == proj_pd ) )
407 elasticLike_system =
true;
412 if ( elasticLike_system ==
true )
415 G4bool elasticLike_energy =
true;
417 for (
G4int i = 0 ; i <
int ( nucleuses.size() ) ; i++ )
425 if ( nucleuses[i]->GetExcitationEnergy()*
GeV > 1.0*
MeV ) elasticLike_energy =
false;
430 G4bool withCollision =
true;
442 if ( elasticLike_energy ==
false ) elastic =
false;
444 if ( elasticLike_energy ==
false && withCollision ==
true ) elastic =
false;
461 if ( elastic ==
true )
464 for ( std::vector< G4QMDNucleus* >::iterator
465 it = nucleuses.begin() ; it != nucleuses.end() ; it++ )
477 for ( std::vector< G4QMDNucleus* >::iterator it
478 = nucleuses.begin() ; it != nucleuses.end() ; it++ )
499 if ( (*it)->GetAtomicNumber() == 0
500 || (*it)->GetAtomicNumber() == (*it)->GetMassNumber() )
503 for (
G4int i = 0 ; i < (*it)->GetTotalNumberOfParticipant() ; i++ )
505 G4QMDParticipant* aP =
new G4QMDParticipant( ( (*it)->GetParticipant( i ) )->GetDefinition() , ( (*it)->GetParticipant( i ) )->GetMomentum() , ( (*it)->GetParticipant( i ) )->GetPosition() );
512 G4LorentzVector nucleus_p4CM ( (*it)->Get4Momentum().vect() , nucleus_e );
516 G4int ia = (*it)->GetMassNumber();
517 G4int iz = (*it)->GetAtomicNumber();
526 for ( G4ReactionProductVector::iterator itt
527 = rv->begin() ; itt != rv->end() ; itt++ )
550 randomized_direction = randomized_direction.
unit();
599 if ( notBreak ==
true )
610 for ( G4ReactionProductVector::iterator itt
611 = rv->begin() ; itt != rv->end() ; itt++ )
643 for ( std::vector< G4QMDNucleus* >::iterator it
644 = nucleuses.begin() ; it != nucleuses.end() ; it++ )
670 G4double stot = std::sqrt ( etot*etot - ptot*ptot );
672 G4double pstt = std::sqrt ( ( stot*stot - ( mass_proj + mass_targ ) * ( mass_proj + mass_targ )
673 ) * ( stot*stot - ( mass_proj - mass_targ ) * ( mass_proj - mass_targ ) ) )
677 G4double eccm = stot - ( mass_proj + mass_targ );
700 G4double rmax = std::sqrt( rmax0*rmax0 + b*b );
703 G4double pcca = 1.0 - double ( zp * zt ) * ccoul / eccm / rmax - ( b / rmax )*( b / rmax );
713 G4double aas = 2.0 * eccm * b / double ( zp * zt ) / ccoul;
714 bbs = 1.0 / std::sqrt ( 1.0 + aas*aas );
715 aas1 = ( 1.0 + aas * b / rmax ) * bbs;
722 if ( 1.0 - aas1*aas1 <= 0 || 1.0 - bbs*bbs <= 0.0 )
729 G4double aat1 = aas1 / std::sqrt ( 1.0 - aas1*aas1 );
730 G4double aat2 = bbs / std::sqrt ( 1.0 - bbs*bbs );
732 thet1 = std::atan ( aat1 );
733 thet2 = std::atan ( aat2 );
737 cost = std::cos( theta );
738 sint = std::sin( theta );
741 G4double rzpr = -rmax * cost * ( mass_targ ) / ( mass_proj + mass_targ );
742 G4double rzta = rmax * cost * ( mass_proj ) / ( mass_proj + mass_targ );
749 G4double pzpc = pzcc * ( cost * pccf + sint * b / rmax );
750 G4double pxpr = pzcc * ( -sint * pccf + cost * b / rmax );
755 G4double epc = std::sqrt ( pzpc*pzpc + pxpr*pxpr + mass_proj*mass_proj );
756 G4double etc = std::sqrt ( pztc*pztc + pxta*pxta + mass_targ*mass_targ );
767 pzpr = pzpc + betacm * gammacm * ( gammacm / ( 1. + gammacm ) * pzpc * betacm + epc );
768 pzta = pztc + betacm * gammacm * ( gammacm / ( 1. + gammacm ) * pztc * betacm + etc );
769 epr = gammacm * ( epc + betacm * pzpc );
770 eta = gammacm * ( etc + betacm * pztc );
775 G4double gammpr = epr / ( mass_proj );
776 G4double gammta = eta / ( mass_targ );
778 pzta = pzta / double ( at );
779 pxta = pxta / double ( at );
781 pzpr = pzpr / double ( ap );
782 pxpr = pxpr / double ( ap );
818 outFile <<
"Lorentz covarianted Quantum Molecular Dynamics model for nucleus (particle) vs nucleus reactions\n";
HepLorentzVector boostOf(const HepLorentzVector &vec, const Hep3Vector &betaVector)
CLHEP::Hep3Vector G4ThreeVector
G4int GetAtomicNumber() const
system("rm -rf microbeam.root")
static constexpr double MeV
void SetMeanField(G4QMDMeanField *meanfield)
static G4IonTable * GetIonTable()
G4double GetTotalPotential()
G4ExcitationHandler * excitationHandler
const G4ParticleDefinition * GetDefinition()
void CalEnergyAndAngularMomentumInCM()
G4int GetTotalNumberOfParticipant()
G4double elastic(Particle const *const p1, Particle const *const p2)
G4LorentzVector Get4Momentum()
Hep3Vector findBoostToCM() const
virtual void ModelDescription(std::ostream &outFile) const
const G4String & GetParticleName() const
const G4String & GetParticleType() const
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState)
void SetTotalPotential(G4double x)
G4double GetPDGCharge() const
G4double GetTotalMomentum() const
static G4CrossSectionDataSetRegistry * Instance()
G4Evaporation * evaporation
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
G4double GetPDGMass() const
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
G4double coulomb_collision_pz_proj
G4double powN(G4double x, G4int n) const
G4double coulomb_collision_rz_targ
void CalKinematicsOfBinaryCollisions(G4double)
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
void SetEvaporation(G4VEvaporation *ptr, G4bool isLocal=false)
G4VCrossSectionDataSet * GetCrossSectionDataSet(const G4String &name, G4bool warning=true)
G4PiNuclearCrossSection * piNucXS
G4double coulomb_collision_gamma_proj
static constexpr double fermi
static G4Pow * GetInstance()
G4ParticleDefinition * GetDefinition() const
G4double coulomb_collision_rx_proj
G4double coulomb_collision_px_proj
G4double coulomb_collision_gamma_targ
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
G4LorentzVector Get4Momentum() const
void DoPropagation(G4double)
G4IonsShenCrossSection * shenXS
std::vector< G4ReactionProduct * > G4ReactionProductVector
void calcOffSetOfCollision(G4double, const G4ParticleDefinition *, const G4ParticleDefinition *, G4double, G4double, G4double, G4ThreeVector)
const G4LorentzVector & Get4Momentum() const
static const char * Default_Name()
G4double GetKineticEnergy() const
G4QMDCollision * collision
G4ThreeVector GetMomentum() const
void SetParticipant(G4QMDParticipant *particle)
G4ThreeVector GetMomentum()
static constexpr double eplus
const G4ParticleDefinition * GetDefinition() const
G4QMDMeanField * meanField
void SetSystem(G4QMDSystem *aSystem)
G4int GetAtomicMass() const
std::vector< G4QMDNucleus * > DoClusterJudgment()
virtual G4double GetIsoCrossSection(const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso=nullptr, const G4Element *elm=nullptr, const G4Material *mat=nullptr)
G4GLOB_DLL std::ostream G4cout
G4double coulomb_collision_pz_targ
static constexpr double pi
G4QMDParticipant * GetParticipant(G4int i)
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
G4double coulomb_collision_rz_proj
G4ThreeVector GetPosition()
void SetNucleus(G4QMDNucleus *aSystem)
void SetDeexChannelsType(G4DeexChannelType val)
static constexpr double GeV
G4double coulomb_collision_px_targ
G4HadFinalState theParticleChange
G4double coulomb_collision_rx_targ
void SetStatusChange(G4HadFinalStateStatus aS)