Remove Particle interface and use ParticleRecord interface instead.

projects/dataclasses/private/InteractionRecord.cxx

@@ -607,20 +607,6 @@ double const & CrossSectionDistributionRecord::GetPrimaryHelicity() const {
     return primary_helicity;
 }
 
-Particle CrossSectionDistributionRecord::GetPrimaryParticle() const {
-    Particle p;
-    p.id = primary_id;
-    p.type = primary_type;
-    p.mass = primary_mass;
-    p.momentum = primary_momentum;
-    p.position = primary_initial_position;
-    p.length = std::sqrt((interaction_vertex.at(0) - primary_initial_position.at(0))*(interaction_vertex.at(0) - primary_initial_position.at(0)) +
-                         (interaction_vertex.at(1) - primary_initial_position.at(1))*(interaction_vertex.at(1) - primary_initial_position.at(1)) +
-                         (interaction_vertex.at(2) - primary_initial_position.at(2))*(interaction_vertex.at(2) - primary_initial_position.at(2)));
-    p.helicity = primary_helicity;
-    return p;
-}
-
 std::array<double, 3> const & CrossSectionDistributionRecord::GetInteractionVertex() const {
     return interaction_vertex;
 }
@@ -669,25 +655,12 @@ void CrossSectionDistributionRecord::SetInteractionParameter(std::string const &
     interaction_parameters[name] = value;
 }
 
-Particle CrossSectionDistributionRecord::GetTargetParticle() const {
-    Particle p;
-    p.id = target_id;
-    p.type = target_type;
-    p.mass = target_mass;
-    p.helicity = target_helicity;
-    return p;
+std::vector<SecondaryParticleRecord> & CrossSectionDistributionRecord::GetSecondaryParticleRecords() {
+    return secondary_particles;
 }
 
-void CrossSectionDistributionRecord::SetTargetParticle(Particle const & particle) {
-    if(particle.id != target_id) {
-        throw std::runtime_error("Cannot set particle with different ID!");
-    }
-    if(particle.type != target_type) {
-        throw std::runtime_error("Cannot set particle with different type!");
-    }
-
-    target_mass = particle.mass;
-    target_helicity = particle.helicity;
+std::vector<SecondaryParticleRecord> const & CrossSectionDistributionRecord::GetSecondaryParticleRecords() const {
+    return secondary_particles;
 }
 
 SecondaryParticleRecord & CrossSectionDistributionRecord::GetSecondaryParticleRecord(size_t index) {
@@ -698,31 +671,6 @@ SecondaryParticleRecord const & CrossSectionDistributionRecord::GetSecondaryPart
     return secondary_particles.at(index);
 }
 
-Particle CrossSectionDistributionRecord::GetSecondaryParticle(size_t index) const {
-    return secondary_particles.at(index).GetParticle();
-}
-
-std::vector<Particle> CrossSectionDistributionRecord::GetSecondaryParticles() const {
-    std::vector<Particle> particles;
-    for(SecondaryParticleRecord const & secondary : secondary_particles) {
-        particles.push_back(secondary.GetParticle());
-    }
-    return particles;
-}
-
-void CrossSectionDistributionRecord::SetSecondaryParticle(size_t index, Particle const & particle) {
-    secondary_particles.at(index).SetParticle(particle);
-}
-
-void CrossSectionDistributionRecord::SetSecondaryParticles(std::vector<Particle> const & particles) {
-    if(particles.size() != secondary_particles.size()) {
-        throw std::runtime_error("Cannot set particles with different size!");
-    }
-    for(size_t i = 0; i < particles.size(); ++i) {
-        secondary_particles.at(i).SetParticle(particles.at(i));
-    }
-}
-
 void CrossSectionDistributionRecord::Finalize(InteractionRecord & record) const {
     record.target_id = target_id;
     record.target_mass = target_mass;

projects/dataclasses/private/pybindings/dataclasses.cxx

@@ -111,14 +111,12 @@ PYBIND11_MODULE(dataclasses,m) {
         .def_property("target_mass", ((double const & (LI::dataclasses::CrossSectionDistributionRecord::*)() const)(&LI::dataclasses::CrossSectionDistributionRecord::GetTargetMass)), &LI::dataclasses::CrossSectionDistributionRecord::SetTargetMass)
         .def_property("target_helicity", ((double const & (LI::dataclasses::CrossSectionDistributionRecord::*)() const)(&LI::dataclasses::CrossSectionDistributionRecord::GetTargetHelicity)), &LI::dataclasses::CrossSectionDistributionRecord::SetTargetHelicity)
         .def_property("interaction_parameters", ((std::map<std::string, double> const & (LI::dataclasses::CrossSectionDistributionRecord::*)())(&LI::dataclasses::CrossSectionDistributionRecord::GetInteractionParameters)), &LI::dataclasses::CrossSectionDistributionRecord::SetInteractionParameters)
-        .def("GetPrimaryParticle", &CrossSectionDistributionRecord::GetPrimaryParticle)
-        .def("GetTargetParticle", &CrossSectionDistributionRecord::GetTargetParticle)
         .def("GetSecondaryParticleRecord",
-                [](LI::dataclasses::CrossSectionDistributionRecord const & cdr, size_t i) {return cdr.GetSecondaryParticleRecord(i);})
-        .def("GetSecondaryParticle", &CrossSectionDistributionRecord::GetSecondaryParticle)
-        .def("GetSecondaryParticles", &CrossSectionDistributionRecord::GetSecondaryParticles)
-        .def("SetSecondaryParticle", &CrossSectionDistributionRecord::SetSecondaryParticle)
-        .def("SetSecondaryParticles", &CrossSectionDistributionRecord::SetSecondaryParticles)
+                [](LI::dataclasses::CrossSectionDistributionRecord & cdr, size_t i) -> LI::dataclasses::SecondaryParticleRecord & {return cdr.GetSecondaryParticleRecord(i);},
+                return_value_policy::reference_internal)
+        .def("GetSecondaryParticleRecords",
+                [](LI::dataclasses::CrossSectionDistributionRecord & cdr) -> std::vector<LI::dataclasses::SecondaryParticleRecord> & {return cdr.GetSecondaryParticleRecords();},
+                return_value_policy::reference_internal)
         .def("Finalize", &CrossSectionDistributionRecord::Finalize);
 
 

projects/dataclasses/public/LeptonInjector/dataclasses/InteractionRecord.h

@@ -69,7 +69,8 @@ class PrimaryDistributionRecord {
 public:
     friend std::ostream& ::operator<<(std::ostream& os, PrimaryDistributionRecord const& record);
 
-    PrimaryDistributionRecord(PrimaryDistributionRecord const & other) = default;
+    PrimaryDistributionRecord(PrimaryDistributionRecord const & other) = delete;
+    PrimaryDistributionRecord(PrimaryDistributionRecord && other) = default;
     PrimaryDistributionRecord & operator=(PrimaryDistributionRecord const & other) = delete;
     PrimaryDistributionRecord & operator=(PrimaryDistributionRecord && other) = delete;
 
@@ -139,7 +140,8 @@ class SecondaryParticleRecord {
 public:
     friend std::ostream& ::operator<<(std::ostream& os, SecondaryParticleRecord const& record);
 
-    SecondaryParticleRecord(SecondaryParticleRecord const & other) = default;
+    SecondaryParticleRecord(SecondaryParticleRecord const & other) = delete;
+    SecondaryParticleRecord(SecondaryParticleRecord && other) = default;
     SecondaryParticleRecord & operator=(SecondaryParticleRecord const & other) = delete;
     SecondaryParticleRecord & operator=(SecondaryParticleRecord && other) = delete;
 
@@ -202,7 +204,8 @@ class CrossSectionDistributionRecord {
 public:
     friend std::ostream& ::operator<<(std::ostream& os, CrossSectionDistributionRecord const& record);
 
-    CrossSectionDistributionRecord(CrossSectionDistributionRecord const & other) = default;
+    CrossSectionDistributionRecord(CrossSectionDistributionRecord const & other) = delete;
+    CrossSectionDistributionRecord(CrossSectionDistributionRecord && other) = default;
     CrossSectionDistributionRecord & operator=(CrossSectionDistributionRecord const & other) = delete;
     CrossSectionDistributionRecord & operator=(CrossSectionDistributionRecord && other) = delete;
 
@@ -215,7 +218,6 @@ class CrossSectionDistributionRecord {
     double const & GetPrimaryMass() const;
     std::array<double, 4> const & GetPrimaryMomentum() const;
     double const & GetPrimaryHelicity() const;
-    Particle GetPrimaryParticle() const;
 
     std::array<double, 3> const & GetInteractionVertex() const;
     ParticleID const & GetTargetID() const;
@@ -232,18 +234,11 @@ class CrossSectionDistributionRecord {
     void SetInteractionParameters(std::map<std::string, double> const & parameters);
     void SetInteractionParameter(std::string const & name, double value);
 
-    Particle GetTargetParticle() const;
-    void SetTargetParticle(Particle const & particle);
-
+    std::vector<SecondaryParticleRecord> & GetSecondaryParticleRecords();
+    std::vector<SecondaryParticleRecord> const & GetSecondaryParticleRecords() const;
     SecondaryParticleRecord & GetSecondaryParticleRecord(size_t index);
     SecondaryParticleRecord const & GetSecondaryParticleRecord(size_t index) const;
 
-    Particle GetSecondaryParticle(size_t index) const;
-    std::vector<Particle> GetSecondaryParticles() const;
-
-    void SetSecondaryParticle(size_t index, Particle const & particle);
-    void SetSecondaryParticles(std::vector<Particle> const & particles);
-
     void Finalize(InteractionRecord & record) const;
 };
 

projects/interactions/private/DISFromSpline.cxx

@@ -544,25 +544,18 @@ void DISFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionRecord
     p3 = p3_lab;
     p4 = p4_lab;
 
-    std::vector<LI::dataclasses::Particle> secondaries = record.GetSecondaryParticles();
-    LI::dataclasses::Particle & lepton = secondaries[lepton_index];
-    LI::dataclasses::Particle & other = secondaries[other_index];
-
-    lepton.momentum[0] = p3.e(); // p3_energy
-    lepton.momentum[1] = p3.px(); // p3_x
-    lepton.momentum[2] = p3.py(); // p3_y
-    lepton.momentum[3] = p3.pz(); // p3_z
-    lepton.mass = p3.m();
-    lepton.helicity = record.primary_helicity;
-
-    other.momentum[0] = p4.e(); // p4_energy
-    other.momentum[1] = p4.px(); // p4_x
-    other.momentum[2] = p4.py(); // p4_y
-    other.momentum[3] = p4.pz(); // p4_z
-    other.mass = p4.m();
-    other.helicity = record.target_helicity;
-
-    record.SetSecondaryParticles(secondaries);
+    std::vector<LI::dataclasses::SecondaryParticleRecord> & secondaries = record.GetSecondaryParticleRecords();
+    LI::dataclasses::SecondaryParticleRecord & lepton = secondaries[lepton_index];
+    LI::dataclasses::SecondaryParticleRecord & other = secondaries[other_index];
+
+
+    lepton.SetFourMomentum({p3.e(), p3.px(), p3.py(), p3.pz()});
+    lepton.SetMass(p3.m());
+    lepton.SetHelicity(record.primary_helicity);
+
+    other.SetFourMomentum({p4.e(), p4.px(), p4.py(), p4.pz()});
+    other.SetMass(p4.m());
+    other.SetHelicity(record.target_helicity);
 }
 
 double DISFromSpline::FinalStateProbability(dataclasses::InteractionRecord const & interaction) const {

projects/interactions/private/DarkNewsCrossSection.cxx

@@ -130,31 +130,22 @@ void DarkNewsCrossSection::SampleFinalState(dataclasses::CrossSectionDistributio
     std::vector<double> secondary_masses = SecondaryMasses(record.signature.secondary_types);
     std::vector<double> secondary_helicities = SecondaryHelicities(record.record);
 
-    std::vector<LI::dataclasses::Particle> secondaries = record.GetSecondaryParticles();
-    for(size_t i = 0; i < secondaries.size(); i++) {
-        secondaries[i].mass = secondary_masses[i];
-        secondaries[i].helicity = secondary_helicities[i];
-    }
-
-    LI::dataclasses::Particle primary = record.GetPrimaryParticle();
-    LI::dataclasses::Particle target = record.GetTargetParticle();
-
     // Uses Metropolis-Hastings Algorithm
     // Assumes we have the differential xsec v.s. Q^2
 
-    rk::P4 p1(geom3::Vector3(primary.momentum[1], primary.momentum[2], primary.momentum[3]), primary.mass);
-    rk::P4 p2(geom3::Vector3(0, 0, 0), target.mass);
+    rk::P4 p1(geom3::Vector3(record.primary_momentum[1], record.primary_momentum[2], record.primary_momentum[3]), record.primary_mass);
+    rk::P4 p2(geom3::Vector3(0, 0, 0), record.target_mass);
 
     // we assume that:
     // the target is stationary so its energy is just its mass
     // the incoming neutrino is massless, so its kinetic energy is its total energy
     // double s = std::pow(rk::invMass(p1, p2), 2);
 
     // define masses that we will use
-    double m1 = primary.mass;
-    double m2 = target.mass;
-    double m3 = secondaries[0].mass;
-    double m4 = secondaries[1].mass;
+    double m1 = record.primary_mass;
+    double m2 = record.target_mass;
+    double m3 = secondary_masses.at(0);
+    double m4 = secondary_masses.at(1);
 
     double primary_energy;
     rk::P4 p1_lab;
@@ -262,18 +253,18 @@ void DarkNewsCrossSection::SampleFinalState(dataclasses::CrossSectionDistributio
     p3 = p3_lab;
     p4 = p4_lab;
 
-    // TODO: helicity update for secondary particles
-    secondaries[0].momentum[0] = p3.e(); // p3_energy
-    secondaries[0].momentum[1] = p3.px(); // p3_x
-    secondaries[0].momentum[2] = p3.py(); // p3_y
-    secondaries[0].momentum[3] = p3.pz(); // p3_z
+    std::vector<LI::dataclasses::SecondaryParticleRecord> & secondaries = record.GetSecondaryParticleRecords();
+
+    LI::dataclasses::SecondaryParticleRecord & p3_record = secondaries[0];
+    LI::dataclasses::SecondaryParticleRecord & p4_record = secondaries[1];
 
-    secondaries[1].momentum[0] = p4.e(); // p4_energy
-    secondaries[1].momentum[1] = p4.px(); // p4_x
-    secondaries[1].momentum[2] = p4.py(); // p4_y
-    secondaries[1].momentum[3] = p4.pz(); // p4_z
+    p3_record.SetFourMomentum({p3.e(), p3.px(), p3.py(), p3.pz()});
+    p3_record.SetMass(secondary_masses.at(0));
+    p3_record.SetHelicity(secondary_helicities.at(0));
 
-    record.SetSecondaryParticles(secondaries);
+    p4_record.SetFourMomentum({p4.e(), p4.px(), p4.py(), p4.pz()});
+    p4_record.SetMass(secondary_masses.at(1));
+    p4_record.SetHelicity(secondary_helicities.at(1));
 }
 
 double DarkNewsCrossSection::FinalStateProbability(dataclasses::InteractionRecord const & record) const {

projects/interactions/private/DipoleFromTable.cxx

@@ -442,20 +442,21 @@ void DipoleFromTable::SampleFinalState(dataclasses::CrossSectionDistributionReco
     else if(channel == Flipping)
         helicity_mul = -1.0;
 
-    std::vector<LI::dataclasses::Particle> secondaries = record.GetSecondaryParticles();
-    secondaries[lepton_index].type = secondary_types[lepton_index];
-    secondaries[lepton_index].mass = p3.m();
-    secondaries[lepton_index].momentum = {p3.e(), p3.px(), p3.py(), p3.pz()};
-    secondaries[lepton_index].length = 0;
-    secondaries[lepton_index].helicity = std::copysign(0.5, record.GetPrimaryHelicity() * helicity_mul);
-
-    secondaries[other_index].type = secondary_types[other_index];
-    secondaries[other_index].mass = p4.m();
-    secondaries[other_index].momentum = {p4.e(), p4.px(), p4.py(), p4.pz()};
-    secondaries[other_index].length = 0;
-    secondaries[other_index].helicity = std::copysign(0.5, record.GetPrimaryHelicity() * helicity_mul);
-
-    record.SetSecondaryParticles(secondaries);
+    std::vector<LI::dataclasses::SecondaryParticleRecord> & secondaries = record.GetSecondaryParticleRecords();
+
+    LI::dataclasses::SecondaryParticleRecord & lepton = secondaries[lepton_index];
+    LI::dataclasses::SecondaryParticleRecord & other = secondaries[other_index];
+
+    assert(lepton.type == secondary_types[lepton_index]);
+    assert(other.type == secondary_types[other_index]);
+
+    lepton.SetMass(p3.m());
+    lepton.SetFourMomentum({p3.e(), p3.px(), p3.py(), p3.pz()});
+    lepton.SetHelicity(std::copysign(0.5, record.GetPrimaryHelicity() * helicity_mul));
+
+    other.SetMass(p4.m());
+    other.SetFourMomentum({p4.e(), p4.px(), p4.py(), p4.pz()});
+    other.SetHelicity(std::copysign(0.5, record.GetPrimaryHelicity() * helicity_mul));
 }
 
 double DipoleFromTable::FinalStateProbability(dataclasses::InteractionRecord const & interaction) const {

projects/interactions/private/DummyCrossSection.cxx

@@ -118,25 +118,17 @@ void DummyCrossSection::SampleFinalState(dataclasses::CrossSectionDistributionRe
     size_t lepton_index = 0;
     size_t other_index = 1;
 
-    std::vector<LI::dataclasses::Particle> secondaries = record.GetSecondaryParticles();
-    LI::dataclasses::Particle & lepton = secondaries[lepton_index];
-    LI::dataclasses::Particle & other = secondaries[other_index];
-
-    lepton.momentum[0] = p3.e(); // p3_energy
-    lepton.momentum[1] = p3.px(); // p3_x
-    lepton.momentum[2] = p3.py(); // p3_y
-    lepton.momentum[3] = p3.pz(); // p3_z
-    lepton.mass = p3.m();
-    lepton.helicity = record.primary_helicity;
-
-    other.momentum[0] = p4.e(); // p4_energy
-    other.momentum[1] = p4.px(); // p4_x
-    other.momentum[2] = p4.py(); // p4_y
-    other.momentum[3] = p4.pz(); // p4_z
-    other.mass = p4.m();
-    other.helicity = record.target_helicity;
-
-    record.SetSecondaryParticles(secondaries);
+    std::vector<LI::dataclasses::SecondaryParticleRecord> & secondaries = record.GetSecondaryParticleRecords();
+    LI::dataclasses::SecondaryParticleRecord & lepton = secondaries[lepton_index];
+    LI::dataclasses::SecondaryParticleRecord & other = secondaries[other_index];
+
+    lepton.SetFourMomentum({p3.e(), p3.px(), p3.py(), p3.pz()});
+    lepton.SetMass(p3.m());
+    lepton.SetHelicity(record.primary_helicity);
+
+    other.SetFourMomentum({p4.e(), p4.px(), p4.py(), p4.pz()});
+    other.SetMass(p4.m());
+    other.SetHelicity(record.target_helicity);
 }
 
 double DummyCrossSection::FinalStateProbability(dataclasses::InteractionRecord const & interaction) const {

projects/interactions/private/ElasticScattering.cxx

@@ -203,25 +203,16 @@ void ElasticScattering::SampleFinalState(dataclasses::CrossSectionDistributionRe
     // doing something dumb, ignore outgoing neutrino
     rk::P4 p4_lab = p1_lab;//p2_lab + (p1_lab - p3_lab);
 
-    std::vector<LI::dataclasses::Particle> secondaries = record.GetSecondaryParticles();
-    LI::dataclasses::Particle & electron = secondaries[electron_index];
-    LI::dataclasses::Particle & neutrino = secondaries[nu_index];
-
-    electron.momentum[0] = p3_lab.e(); // p3_energy
-    electron.momentum[1] = p3_lab.px(); // p3_x
-    electron.momentum[2] = p3_lab.py(); // p3_y
-    electron.momentum[3] = p3_lab.pz(); // p3_z
-    electron.mass = p3_lab.m();
-    electron.helicity = record.target_helicity;
-
-    neutrino.momentum[0] = p4_lab.e(); // p4_energy
-    neutrino.momentum[1] = p4_lab.px(); // p4_x
-    neutrino.momentum[2] = p4_lab.py(); // p4_y
-    neutrino.momentum[3] = p4_lab.pz(); // p4_z
-    neutrino.mass = p4_lab.m();
-    neutrino.helicity = record.primary_helicity;
-
-    record.SetSecondaryParticles(secondaries);
+    LI::dataclasses::SecondaryParticleRecord & electron = record.GetSecondaryParticleRecord(electron_index);
+    LI::dataclasses::SecondaryParticleRecord & neutrino = record.GetSecondaryParticleRecord(nu_index);
+
+    electron.SetFourMomentum({p3_lab.e(), p3_lab.px(), p3_lab.py(), p3_lab.pz()});
+    electron.SetMass(p3_lab.m());
+    electron.SetHelicity(record.target_helicity);
+
+    neutrino.SetFourMomentum({p4_lab.e(), p4_lab.px(), p4_lab.py(), p4_lab.pz()});
+    neutrino.SetMass(p4_lab.m());
+    neutrino.SetHelicity(record.primary_helicity);
 }
 
 std::vector<LI::dataclasses::Particle::ParticleType> ElasticScattering::GetPossibleTargets() const {