Merge pull request #484 from LeanderSchlegel/fix_mass_inheritance

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;
 	}
 }

test/testCore.cpp

@@ -149,6 +149,16 @@ TEST(ParticleID, isNucleus) {
 	EXPECT_FALSE(isNucleus(11));
 }
 
+TEST(ParticleMass, particleMass) {
+	//particleMass(int id) interfaces nuclearMass for nuclei
+	EXPECT_DOUBLE_EQ(nuclearMass(nucleusId(1,1)), particleMass(nucleusId(1,1)));
+	//particleMass(int id) for electron/positron, photon and neutrino
+	EXPECT_DOUBLE_EQ(mass_electron,particleMass(11));
+	EXPECT_DOUBLE_EQ(mass_electron,particleMass(-11));
+	EXPECT_DOUBLE_EQ(0.0,particleMass(22));
+	EXPECT_DOUBLE_EQ(0.0,particleMass(14));
+}
+
 TEST(HepPID, consistencyWithReferenceImplementation) {
 	// Tests the performance improved version against the default one
 	unsigned long testPID = rand() % 1000000000 + 1000000000;
@@ -219,8 +229,14 @@ TEST(Candidate, addSecondary) {
 
 	c.addSecondary(nucleusId(1,1), 200);
 	c.addSecondary(nucleusId(1,1), 200, 5.);
-	Candidate s1 = *c.secondaries[0];
-	Candidate s2 = *c.secondaries[1];
+	c.addSecondary(11, 200);
+	c.addSecondary(14, 200);
+	c.addSecondary(22, 200);
+	Candidate s1 = *c.secondaries[0]; //proton
+	Candidate s2 = *c.secondaries[1]; //proton
+	Candidate s3 = *c.secondaries[2]; //electron
+	Candidate s4 = *c.secondaries[3]; //neutrino
+	Candidate s5 = *c.secondaries[4]; //photon
 
 	EXPECT_EQ(nucleusId(1,1), s1.current.getId());
 	EXPECT_EQ(200, s1.current.getEnergy());
@@ -231,6 +247,9 @@ TEST(Candidate, addSecondary) {
 	EXPECT_TRUE(Vector3d(1,2,3) == s1.created.getPosition());
 	EXPECT_TRUE(Vector3d(0,0,1) == s1.created.getDirection());
 	EXPECT_TRUE(s1.getTagOrigin() == "SEC");
+	EXPECT_EQ(mass_electron,s3.current.getMass());
+	EXPECT_EQ(0.0,s4.current.getMass());
+	EXPECT_EQ(0.0,s5.current.getMass());
 
 	EXPECT_EQ(15., s2.getWeight());
 }