73{
74 auto & tA = const_cast<G4Track&>(trackA);
75 auto & tB = const_cast<G4Track&>(trackB);
77
78 std::unique_ptr<G4ITReactionChange> pChanges(new G4ITReactionChange());
79 pChanges->Initialize(trackA, trackB);
80
83
84 const auto pReactionData =
fMolReactionTable->GetReactionData(pMoleculeA, pMoleculeB);
85 const G4int nbProducts = pReactionData->GetNbProducts();
86 if (nbProducts != 0)
87 {
88 const G4double D1 = pMoleculeA->GetDiffusionCoefficient();
89 const G4double D2 = pMoleculeB->GetDiffusionCoefficient();
90 const G4double sqrD1 = D1 == 0. ? 0. : std::sqrt(D1);
91 const G4double sqrD2 = D2 == 0. ? 0. : std::sqrt(D2);
92 const G4double inv_numerator = 1./(sqrD1 + sqrD2);
93 const G4ThreeVector reactionSite = sqrD2 * inv_numerator * tA.GetPosition()
94 + sqrD1 * inv_numerator * tB.GetPosition();
95
97 auto randP = (1-u) * tA.GetPosition() + u * tB.GetPosition();
98
99 for (
G4int j = 0; j < nbProducts; ++j)
100 {
101 auto pProduct = new G4Molecule(pReactionData->GetProduct(j));
102 auto pProductTrack = pProduct->BuildTrack(trackA.
GetGlobalTime(), (reactionSite + randP)/2);
103 pProductTrack->SetTrackStatus(
fAlive);
105 pChanges->AddSecondary(pProductTrack);
106 }
107 }
108 pChanges->KillParents(true);
109 return pChanges;
110}
G4Molecule * GetMolecule(const G4Track &track)
CLHEP::Hep3Vector G4ThreeVector
void UpdatePositionForReaction(G4Track &, G4Track &)
void Push(G4Track *) override
static G4ITTrackHolder * Instance()
const G4MolecularConfiguration * GetMolecularConfiguration() const
G4double GetGlobalTime() const