From 9a551d6f0e8a32b36704c0b3bf0e5d8dca45c516 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sat, 24 Feb 2024 09:10:05 -0700 Subject: [PATCH] Remove Particle interface and use ParticleRecord interface instead. --- .../dataclasses/private/InteractionRecord.cxx | 60 ++----------------- .../private/pybindings/dataclasses.cxx | 12 ++-- .../dataclasses/InteractionRecord.h | 21 +++---- .../interactions/private/DISFromSpline.cxx | 31 ++++------ .../private/DarkNewsCrossSection.cxx | 41 +++++-------- .../interactions/private/DipoleFromTable.cxx | 29 ++++----- .../private/DummyCrossSection.cxx | 30 ++++------ .../private/ElasticScattering.cxx | 29 ++++----- .../interactions/private/HNLFromSpline.cxx | 31 ++++------ .../interactions/private/NeutrissimoDecay.cxx | 39 ++++++------ 10 files changed, 115 insertions(+), 208 deletions(-) diff --git a/projects/dataclasses/private/InteractionRecord.cxx b/projects/dataclasses/private/InteractionRecord.cxx index 6ed87985..d2c791b6 100644 --- a/projects/dataclasses/private/InteractionRecord.cxx +++ b/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 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 & 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 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 CrossSectionDistributionRecord::GetSecondaryParticles() const { - std::vector 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 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; diff --git a/projects/dataclasses/private/pybindings/dataclasses.cxx b/projects/dataclasses/private/pybindings/dataclasses.cxx index 9a578d80..a5c89bd9 100644 --- a/projects/dataclasses/private/pybindings/dataclasses.cxx +++ b/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 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 & {return cdr.GetSecondaryParticleRecords();}, + return_value_policy::reference_internal) .def("Finalize", &CrossSectionDistributionRecord::Finalize); diff --git a/projects/dataclasses/public/LeptonInjector/dataclasses/InteractionRecord.h b/projects/dataclasses/public/LeptonInjector/dataclasses/InteractionRecord.h index eb97b870..cbf69d11 100644 --- a/projects/dataclasses/public/LeptonInjector/dataclasses/InteractionRecord.h +++ b/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 const & GetPrimaryMomentum() const; double const & GetPrimaryHelicity() const; - Particle GetPrimaryParticle() const; std::array const & GetInteractionVertex() const; ParticleID const & GetTargetID() const; @@ -232,18 +234,11 @@ class CrossSectionDistributionRecord { void SetInteractionParameters(std::map const & parameters); void SetInteractionParameter(std::string const & name, double value); - Particle GetTargetParticle() const; - void SetTargetParticle(Particle const & particle); - + std::vector & GetSecondaryParticleRecords(); + std::vector const & GetSecondaryParticleRecords() const; SecondaryParticleRecord & GetSecondaryParticleRecord(size_t index); SecondaryParticleRecord const & GetSecondaryParticleRecord(size_t index) const; - Particle GetSecondaryParticle(size_t index) const; - std::vector GetSecondaryParticles() const; - - void SetSecondaryParticle(size_t index, Particle const & particle); - void SetSecondaryParticles(std::vector const & particles); - void Finalize(InteractionRecord & record) const; }; diff --git a/projects/interactions/private/DISFromSpline.cxx b/projects/interactions/private/DISFromSpline.cxx index 07173ad6..7d69d1ca 100644 --- a/projects/interactions/private/DISFromSpline.cxx +++ b/projects/interactions/private/DISFromSpline.cxx @@ -544,25 +544,18 @@ void DISFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionRecord p3 = p3_lab; p4 = p4_lab; - std::vector 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 & 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 { diff --git a/projects/interactions/private/DarkNewsCrossSection.cxx b/projects/interactions/private/DarkNewsCrossSection.cxx index 727a6910..420192cb 100644 --- a/projects/interactions/private/DarkNewsCrossSection.cxx +++ b/projects/interactions/private/DarkNewsCrossSection.cxx @@ -130,20 +130,11 @@ void DarkNewsCrossSection::SampleFinalState(dataclasses::CrossSectionDistributio std::vector secondary_masses = SecondaryMasses(record.signature.secondary_types); std::vector secondary_helicities = SecondaryHelicities(record.record); - std::vector 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 @@ -151,10 +142,10 @@ void DarkNewsCrossSection::SampleFinalState(dataclasses::CrossSectionDistributio // 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 & 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 { diff --git a/projects/interactions/private/DipoleFromTable.cxx b/projects/interactions/private/DipoleFromTable.cxx index 7ba3adf0..2274bc18 100644 --- a/projects/interactions/private/DipoleFromTable.cxx +++ b/projects/interactions/private/DipoleFromTable.cxx @@ -442,20 +442,21 @@ void DipoleFromTable::SampleFinalState(dataclasses::CrossSectionDistributionReco else if(channel == Flipping) helicity_mul = -1.0; - std::vector 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 & 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 { diff --git a/projects/interactions/private/DummyCrossSection.cxx b/projects/interactions/private/DummyCrossSection.cxx index cf6fb9ec..e1b91b90 100644 --- a/projects/interactions/private/DummyCrossSection.cxx +++ b/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 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 & 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 { diff --git a/projects/interactions/private/ElasticScattering.cxx b/projects/interactions/private/ElasticScattering.cxx index 7f6b0a8b..d84f3edf 100644 --- a/projects/interactions/private/ElasticScattering.cxx +++ b/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 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 ElasticScattering::GetPossibleTargets() const { diff --git a/projects/interactions/private/HNLFromSpline.cxx b/projects/interactions/private/HNLFromSpline.cxx index a6dc7e96..b38421d3 100644 --- a/projects/interactions/private/HNLFromSpline.cxx +++ b/projects/interactions/private/HNLFromSpline.cxx @@ -542,25 +542,18 @@ void HNLFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionRecord p3 = p3_lab; p4 = p4_lab; - std::vector 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 & 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 HNLFromSpline::FinalStateProbability(dataclasses::InteractionRecord const & interaction) const { diff --git a/projects/interactions/private/NeutrissimoDecay.cxx b/projects/interactions/private/NeutrissimoDecay.cxx index 750c8383..d521d7ef 100644 --- a/projects/interactions/private/NeutrissimoDecay.cxx +++ b/projects/interactions/private/NeutrissimoDecay.cxx @@ -135,7 +135,6 @@ double NeutrissimoDecay::DifferentialDecayWidth(dataclasses::InteractionRecord c void NeutrissimoDecay::SampleFinalState(dataclasses::CrossSectionDistributionRecord & record, std::shared_ptr random) const { LI::dataclasses::InteractionSignature const & signature = record.GetSignature(); - std::vector secondaries = record.GetSecondaryParticles(); unsigned int gamma_index = (signature.secondary_types[0] == LI::dataclasses::Particle::ParticleType::Gamma) ? 0 : 1; unsigned int nu_index = 1 - gamma_index; @@ -143,12 +142,13 @@ void NeutrissimoDecay::SampleFinalState(dataclasses::CrossSectionDistributionRec double CosTheta; double alpha = std::copysign(1.0,record.GetPrimaryHelicity()); // 1 for RH, -1 for LH alpha = (signature.primary_type == LI::dataclasses::Particle::ParticleType::NuF4) ? -1*alpha : alpha; - if(nature==ChiralNature::Majorana) { - CosTheta = random->Uniform(-1,1); + + if(nature == ChiralNature::Majorana) { + CosTheta = random->Uniform(-1,1); } else { - double X = random->Uniform(0,1); - CosTheta = (std::sqrt(1 - 2*alpha*(1 - alpha/2. - 2*X)) - 1)/alpha; + double X = random->Uniform(0,1); + CosTheta = (std::sqrt(1 - 2*alpha*(1 - alpha/2. - 2*X)) - 1)/alpha; } double SinTheta = std::sin(std::acos(CosTheta)); @@ -170,20 +170,25 @@ void NeutrissimoDecay::SampleFinalState(dataclasses::CrossSectionDistributionRec rk::P4 pGamma = pGamma_HNLrest.boost(boost_to_lab); rk::P4 pNu(pHNL.momentum() - pGamma.momentum(),0); // ensures the neutrino has zero mass, avoids rounding errors - secondaries[gamma_index].momentum[0] = pGamma.e(); // pGamma_energy - secondaries[gamma_index].momentum[1] = pGamma.px(); // pGamma_x - secondaries[gamma_index].momentum[2] = pGamma.py(); // pGamma_y - secondaries[gamma_index].momentum[3] = pGamma.pz(); // pGamma_z - secondaries[gamma_index].mass = pGamma.m(); - secondaries[nu_index].momentum[0] = pNu.e(); // pNu_energy - secondaries[nu_index].momentum[1] = pNu.px(); // pNu_x - secondaries[nu_index].momentum[2] = pNu.py(); // pNu_y - secondaries[nu_index].momentum[3] = pNu.pz(); // pNu_z - secondaries[nu_index].mass = pNu.m(); - secondaries[nu_index].helicity = -1*record.primary_helicity; + LI::dataclasses::SecondaryParticleRecord & gamma = record.GetSecondaryParticleRecord(gamma_index); + LI::dataclasses::SecondaryParticleRecord & nu = record.GetSecondaryParticleRecord(nu_index); + + assert(gamma.type == LI::dataclasses::Particle::ParticleType::Gamma); + assert(nu.type == LI::dataclasses::Particle::ParticleType::NuE || + nu.type == LI::dataclasses::Particle::ParticleType::NuMu || + nu.type == LI::dataclasses::Particle::ParticleType::NuTau || + nu.type == LI::dataclasses::Particle::ParticleType::NuEBar || + nu.type == LI::dataclasses::Particle::ParticleType::NuMuBar || + nu.type == LI::dataclasses::Particle::ParticleType::NuTauBar); + + gamma.SetFourMomentum({pGamma.e(), pGamma.px(), pGamma.py(), pGamma.pz()}); + gamma.SetMass(pGamma.m()); + gamma.SetHelicity(std::copysign(1.0, record.primary_helicity)); - record.SetSecondaryParticles(secondaries); + nu.SetFourMomentum({pNu.e(), pNu.px(), pNu.py(), pNu.pz()}); + nu.SetMass(pNu.m()); + nu.SetHelicity(-1*record.primary_helicity); } double NeutrissimoDecay::FinalStateProbability(dataclasses::InteractionRecord const & record) const {