188 const G4DynamicParticle* aDynamicParticle = aTrack.GetDynamicParticle();
195 G4Material* aMaterial = aTrack.GetMaterial();
196 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
197 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
199 G4double NumberOfProtons = 0;
200 G4double NumberOfNucleons = 0;
202 for (
size_t elm = 0; elm < aMaterial->GetNumberOfElements(); elm++) {
204 NumberOfProtons += NbOfAtomsPerVolume[elm] * (*theElementVector)[elm]->GetZ();
206 NumberOfNucleons += NbOfAtomsPerVolume[elm] * (*theElementVector)[elm]->GetN();
209 if (G4UniformRand() < NumberOfProtons / NumberOfNucleons) {
218 G4int theIncidentPDG = aDynamicParticle->GetDefinition()->GetPDGEncoding();
221 aDynamicParticle->GetDefinition()->GetPDGCharge() == 0.) {
223 theIncidentPDG *= -1;
226 bool baryonise =
false;
247 std::vector<bool> theChargeChangeList;
249 for (ReactionProductList::iterator prod_it = aReactionProductList->begin(); prod_it != aReactionProductList->end();
251 G4int secondaries = prod_it->size();
258 theReactionProductList.push_back(*prod_it);
259 if (secondaries == 2) {
261 }
else if (secondaries == 3) {
264 G4cerr <<
"ReactionProduct has unsupported number of secondaries: " << secondaries << G4endl;
272 if (theReactionProductList.empty())
273 G4Exception(
"G4ProcessHelper",
276 "GetFinalState: No process could be selected from the given list.");
280 int n_rps = theReactionProductList.size();
281 int select = (int)(G4UniformRand() * n_rps);
283 return theReactionProductList[
select];
289 G4double p23 = 1 - p22;
291 std::vector<G4double> Probabilities;
292 std::vector<G4bool> TwotoThreeFlag;
294 G4double CumulatedProbability = 0;
298 for (
unsigned int i = 0;
i != theReactionProductList.size();
i++) {
299 if (theReactionProductList[
i].
size() == 2) {
300 CumulatedProbability += p22 / N22;
301 TwotoThreeFlag.push_back(
false);
303 CumulatedProbability += p23 / N23;
304 TwotoThreeFlag.push_back(
true);
306 Probabilities.push_back(CumulatedProbability);
312 for (std::vector<G4double>::iterator it = Probabilities.begin(); it != Probabilities.end(); it++) {
313 *it /= CumulatedProbability;
320 G4bool selected =
false;
326 while (!selected && tries < 100) {
328 G4double dice = G4UniformRand();
330 while (dice > Probabilities[i] && i < theReactionProductList.size()) {
337 if (!TwotoThreeFlag[i]) {
342 if (
PhaseSpace(theReactionProductList[i], aDynamicParticle) > G4UniformRand())
347 if (selected &&
particleTable->FindParticle(theReactionProductList[i][0])->GetPDGCharge() !=
348 aDynamicParticle->GetDefinition()->GetPDGCharge()) {
362 G4cerr <<
"Could not select process!!!!" << G4endl;
368 if (theReactionProductList[i].
size() == 2) {
378 return theReactionProductList[
i];
G4ParticleTable * particleTable
G4ParticleDefinition * theProton
static bool s_isSbaryon(int pdg)
G4ParticleDefinition * theNeutron
static bool s_isRMeson(int pdg)
G4double PhaseSpace(const ReactionProduct &aReaction, const G4DynamicParticle *aDynamicParticle)
static bool s_isMesonino(int pdg)
static bool s_isRBaryon(int pdg)
ReactionMap * theReactionMap
G4ParticleDefinition * theTarget
G4bool ReactionIsPossible(const ReactionProduct &aReaction, const G4DynamicParticle *aDynamicParticle)
std::vector< ReactionProduct > ReactionProductList
tuple size
Write out results.
G4bool ReactionGivesBaryon(const ReactionProduct &aReaction)