Possible fix for wrong secondary masses of photons and neutrinos

CHANGELOG.md

@@ -3,8 +3,10 @@
 ### Bug fixes:
  * Fixed sign for exponential decay of magn. field strength with Galactic height in LogarithmicSpiralField 
  * Fixed r term in source distribution for SNR and Pulsar 
+ * Fixed wrong mass inheritance for secondaries other than nuclei or electron/positron
 
 ### New features:
+ * Added new backwards-compatible function particleMass that returns particle mass also for non-nuclei 
 
 ### Interface changes:
 

include/crpropa/ParticleMass.h

@@ -6,7 +6,16 @@ namespace crpropa {
  * \addtogroup PhysicsDefinitions
  * @{
  */
-
+
+ /** Get the particle mass by lookup from a table.
+ For nuclei, the function nuclearMass is called, for the case of
+ electrons or positrons the mass_electron is returned and for all
+ other cases like photons and also neutrinos, zero mass is returned.
+ @param id		id of the particle following the PDG numbering scheme
+ @returns The mass of a the particle
+ */
+ double particleMass(int id);
+
 /** Get the nucleus mass by lookup from a table.
  The masses are the atomic masses from the NIST database: 
    http://www.nist.gov/pml/data/comp.cfm
@@ -17,6 +26,7 @@ namespace crpropa {
  @returns The mass of a the nucleus
  */
 double nuclearMass(int id);
+
 /** Get the nucleus mass by lookup from a table.
  The masses are the atomic masses from the NIST database: 
    http://www.nist.gov/pml/data/comp.cfm

src/ParticleMass.cpp

@@ -53,6 +53,14 @@ struct NuclearMassTable {
 
 static NuclearMassTable nuclearMassTable;
 
+double particleMass(int id) {
+	if (isNucleus(id))
+		return nuclearMass(id);
+	if (abs(id) == 11)
+		return mass_electron;
+	return 0.0;
+}
+
 double nuclearMass(int id) {
 	int A = massNumber(id);
 	int Z = chargeNumber(id);

src/ParticleState.cpp

@@ -49,14 +49,12 @@ double ParticleState::getRigidity() const {
 
 void ParticleState::setId(int newId) {
 	id = newId;
+	pmass = particleMass(id);
 	if (isNucleus(id)) {
-		pmass = nuclearMass(id);
 		charge = chargeNumber(id) * eplus;
 		if (id < 0)
 			charge *= -1; // anti-nucleus
 	} else {
-		if (abs(id) == 11)
-			pmass = mass_electron;
 		charge = HepPID::charge(id) * eplus;
 	}
 }