Pion-absorption bug in 10.2.3 Geant4 Bertini cascade

[EinarElen]

While working on some of the Bertini cascade internals I saw reactions like the following
$$\pi^+ + (nn) \rightarrow n + n$$
$$\pi^+ + (pn) \rightarrow p + n$$
$$\pi^- + (pn) \rightarrow p + n$$
$$\pi^- + (pp) \rightarrow p + p$$
etc as part of the cascade which clearly violates charge conservation. This particular reaction happens really often. I made some quick bookkeeping by running the event generator directly for 3 and 6 GeV photons on W.

Bookkeeping for 3 GeV photons 
For 5000 events
1091 had 1 non-charge conserving interactions -> ~21.8 %
319 had 2 non-charge conserving interactions -> ~6.36 %
59 had 3 non-charge conserving interactions -> ~1.18 %
8 had 4 non-charge conserving interactions -> ~0.16 %
2 had 5 non-charge conserving interactions -> ~0.0399 %
0 had 6 or more non-charge conserving interactions -> ~0 %
In total 1479 events had at least 1 non-charge conserving interaction -> ~29.5 %

and

Bookkeeping for 6 GeV photons 
For 5000 events
1207 had 1 non-charge conserving interactions -> ~24.1 %
462 had 2 non-charge conserving interactions -> ~9.21 %
180 had 3 non-charge conserving interactions -> ~3.59 %
43 had 4 non-charge conserving interactions -> ~0.857 %
17 had 5 non-charge conserving interactions -> ~0.339 %
6 had 6 non-charge conserving interactions -> ~0.12 %
0 had 7 non-charge conserving interactions -> ~0 %

Note: This takes the number of attempts for each event into account!

The way that quasi-deutron absorption is modelled for pions is to "explode" the quasi-nucleon pair. If the interaction is allowed, and selected the collision output is just the two nucleons. This of course only can happen if the pion and the quasi-deutron matches up with the final state.

After some debugging and making sure that I wasn't missing something, this happens because the charge conservation check in

geant4/source/processes/hadronic/models/cascade/cascade/src/G4ElementaryParticleCollider.cc

Line 401 in 6271441

G4ElementaryParticleCollider::generateSCMpionAbsorption(G4double etot_scm,

is wrong

geant4/source/processes/hadronic/models/cascade/cascade/src/G4ElementaryParticleCollider.cc

Lines 417 to 425 in 6271441

	if (!G4NucleiModel::useQuasiDeuteron(type1, type2)) {
	G4cerr << " pion absorption: "
	<< particle1->getDefinition()->GetParticleName() << " + "
	<< particle2->getDefinition()->GetParticleName() << " -> ?"
	<< G4endl;
	return;
	}
	if (!splitQuasiDeuteron(type2)) return; // Get constituents of [NN]

The check is here

geant4/source/processes/hadronic/models/cascade/cascade/src/G4NucleiModel.cc

Lines 1047 to 1058 in 6271441

	// Test if particle is suitable for absorption with dibaryons
	G4bool G4NucleiModel::useQuasiDeuteron(G4int ptype, G4int qdtype) {
	if (qdtype==pn || qdtype==0) // All absorptive particles
	return (ptype==pi0 || ptype==pip || ptype==pim || ptype==gam || ptype==mum);
	else if (qdtype==pp) // Negative or neutral only
	return (ptype==pi0 || ptype==pim || ptype==gam || ptype==mum);
	else if (qdtype==nn) // Positive or neutral only
	return (ptype==pi0 || ptype==pip || ptype==gam);
	return false; // Not a quasideuteron target
	}

This was patched in Geant4 10.5.0

geant4/source/processes/hadronic/models/cascade/cascade/src/G4ElementaryParticleCollider.cc

Lines 569 to 579 in db49709

	G4int type1 = particle1->type();
	G4int type2 = particle2->type();
	// Ensure that single-nucleon absportion is valid (charge exchangeable)
	if ((type1*type2 != pim*pro && type1*type2 != pip*neu)) {
	G4cerr << " pion-nucleon absorption: "
	<< particle1->getDefinition()->GetParticleName() << " + "
	<< particle2->getDefinition()->GetParticleName() << " -> ?"
	<< G4endl;
	return;
	}

Please don't check what the local Friday evening time in Lund is right now...

NB: I don't think this affects photo-absorption by quasi-deutrons but I haven't checked explicitly

[EinarElen]

On second thought, this obviously doesn't affect photo-absorption