diff --git a/CHANGELOG.md b/CHANGELOG.md index 4332cf998..4d7d0ccc3 100644 --- a/CHANGELOG.md +++ b/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: diff --git a/include/crpropa/ParticleMass.h b/include/crpropa/ParticleMass.h index f4e5a6d6c..06175f20c 100644 --- a/include/crpropa/ParticleMass.h +++ b/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 diff --git a/src/ParticleMass.cpp b/src/ParticleMass.cpp index a2a5c5090..075668487 100644 --- a/src/ParticleMass.cpp +++ b/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); diff --git a/src/ParticleState.cpp b/src/ParticleState.cpp index 9f867c384..85736a856 100644 --- a/src/ParticleState.cpp +++ b/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; } } diff --git a/test/testCore.cpp b/test/testCore.cpp index 7aebff63a..b6738ac2a 100644 --- a/test/testCore.cpp +++ b/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()); }