Scalar dimuon

.github/workflows/build_wheels.yml

@@ -3,6 +3,7 @@ name: Build
 on: [pull_request]
 
 jobs:
+
   build_wheels_macos:
     name: Build wheels on ${{ matrix.os }}
     runs-on: ${{ matrix.os }}

projects/dataclasses/public/SIREN/dataclasses/ParticleTypes.def

@@ -196,3 +196,4 @@ X(STauPlus, -2000009131)
 X(STauMinus, -2000009132)
 X(SMPPlus, -2000009500)
 X(SMPMinus, -2000009501)
+X(Scalar, 2000000001)

projects/interactions/CMakeLists.txt

@@ -16,6 +16,7 @@ LIST (APPEND interactions_SOURCES
     ${PROJECT_SOURCE_DIR}/projects/interactions/private/DummyCrossSection.cxx
     ${PROJECT_SOURCE_DIR}/projects/interactions/private/pyDarkNewsCrossSection.cxx
     ${PROJECT_SOURCE_DIR}/projects/interactions/private/pyDarkNewsDecay.cxx
+    ${PROJECT_SOURCE_DIR}/projects/interactions/private/ScalarDecay.cxx
 )
 add_library(SIREN_interactions OBJECT ${interactions_SOURCES})
 set_property(TARGET SIREN_interactions PROPERTY POSITION_INDEPENDENT_CODE ON)

projects/interactions/private/Decay.cxx

@@ -22,9 +22,10 @@ bool Decay::operator==(Decay const & other) const {
 
 double Decay::TotalDecayLength(dataclasses::InteractionRecord const & interaction) const {
     double tau = 1./TotalDecayWidth(interaction); // in inverse GeV
-    std::array<double, 4> const & p4 = interaction.primary_momentum;
+    std::array<double, 4> const & p4 = interaction.primary_momentum; // only energy is gauranteed to be set here!
     double const & mass = interaction.primary_mass;
-    rk::P4 p1(geom3::Vector3(p4[1], p4[2], p4[3]), mass);
+    double momentum = sqrt(p4[0]*p4[0] - mass*mass); // three-momentum
+    rk::P4 p1(geom3::Vector3(0, 0, momentum), mass); // just assume along z direction for now
     return p1.beta() * p1.gamma() * tau * siren::utilities::Constants::hbarc;
 }
 

projects/interactions/private/ScalarDecay.cxx

@@ -0,0 +1,147 @@
+#include "SIREN/interactions/ScalarDecay.h"
+
+#include <array>                                              // for array
+#include <cmath>                                              // for copysign
+#include <tuple>                                              // for tie
+#include <string>                                             // for basic_s...
+
+#include <rk/geom3.hh>                                        // for Vector3
+#include <rk/rk.hh>                                           // for P4, Boost
+
+#include "SIREN/interactions/Decay.h"               // for Decay
+#include "SIREN/dataclasses/InteractionRecord.h"     // for Interac...
+#include "SIREN/dataclasses/InteractionSignature.h"  // for Interac...
+#include "SIREN/dataclasses/Particle.h"              // for Particle
+#include "SIREN/math/Vector3D.h"                     // for Vector3D
+#include "SIREN/utilities/Constants.h"               // for GeV, pi
+#include "SIREN/utilities/Random.h"                  // for SIREN_random
+
+namespace siren {
+namespace interactions {
+
+bool ScalarDecay::equal(Decay const & other) const {
+    const ScalarDecay* x = dynamic_cast<const ScalarDecay*>(&other);
+
+    if(!x)
+        return false;
+    else
+        return
+            std::tie(
+                    primary_types,
+                    mass,
+                    coupling)
+            ==
+            std::tie(
+                    x->primary_types,
+                    x->mass,
+                    x->coupling);
+}
+
+double ScalarDecay::TotalDecayWidth(dataclasses::InteractionRecord const & record) const {
+    return TotalDecayWidth(record.signature.primary_type);
+}
+
+double ScalarDecay::TotalDecayWidth(siren::dataclasses::ParticleType primary) const {
+    return coupling*coupling * mass / (16*siren::utilities::Constants::pi) * siren::utilities::Constants::GeV;
+}
+
+double ScalarDecay::TotalDecayWidthForFinalState(dataclasses::InteractionRecord const & record) const {
+    siren::dataclasses::InteractionSignature const & signature = record.signature;
+    if(signature.secondary_types.size()!=2) return 0;
+    // make sure this is a dimuon decay
+    if(!(signature.secondary_types[0]==siren::dataclasses::ParticleType::MuPlus ||
+         signature.secondary_types[0]==siren::dataclasses::ParticleType::MuMinus)) return 0;
+    if(!(signature.secondary_types[1]==siren::dataclasses::ParticleType::MuPlus ||
+         signature.secondary_types[1]==siren::dataclasses::ParticleType::MuMinus)) return 0;
+    // if so, return the total decay width
+    // this could be updated to just the dimuon decay width in the future
+    // right now, requires user to apply branching ratio on the back end
+    return TotalDecayWidth(signature.primary_type);
+}
+
+std::vector<std::string> ScalarDecay::DensityVariables() const {
+    return std::vector<std::string>{"CosTheta"};
+}
+
+
+std::vector<dataclasses::InteractionSignature> ScalarDecay::GetPossibleSignatures() const {
+    std::vector<dataclasses::InteractionSignature> signatures;
+    for(auto primary : primary_types) {
+        std::vector<dataclasses::InteractionSignature> new_signatures = GetPossibleSignaturesFromParent(primary);
+        signatures.insert(signatures.end(),new_signatures.begin(),new_signatures.end());
+    }
+    return signatures;
+}
+
+std::vector<dataclasses::InteractionSignature> ScalarDecay::GetPossibleSignaturesFromParent(siren::dataclasses::ParticleType primary) const {
+    std::vector<dataclasses::InteractionSignature> signatures;
+    dataclasses::InteractionSignature signature;
+    signature.primary_type = primary;
+    signature.target_type = siren::dataclasses::ParticleType::Decay;
+    signature.secondary_types.resize(2);
+    if(primary==siren::dataclasses::ParticleType::Scalar) {
+        signature.secondary_types[0] = siren::dataclasses::ParticleType::MuPlus;
+        signature.secondary_types[1] = siren::dataclasses::ParticleType::MuMinus;
+        signatures.push_back(signature);
+    }
+    return signatures;
+}
+
+double ScalarDecay::DifferentialDecayWidth(dataclasses::InteractionRecord const & record) const {
+    double DecayWidth = TotalDecayWidthForFinalState(record);
+    return DecayWidth/2.; // for costheta between -1,1
+}
+
+void ScalarDecay::SampleFinalState(dataclasses::CrossSectionDistributionRecord & record, std::shared_ptr<siren::utilities::SIREN_random> random) const {
+
+
+    double CosTheta = random->Uniform(-1,1);
+    double SinTheta = std::sin(std::acos(CosTheta));
+
+    rk::P4 pScalar(geom3::Vector3(record.primary_momentum[1], record.primary_momentum[2], record.primary_momentum[3]), record.primary_mass);
+    rk::Boost boost_to_lab = pScalar.labBoost();
+
+    geom3::UnitVector3 x_dir = geom3::UnitVector3::xAxis();
+    geom3::Vector3 pScalar_mom = pScalar.momentum();
+    geom3::UnitVector3 pScalar_dir = pScalar_mom.direction();
+    geom3::Rotation3 x_to_pScalar_rot = geom3::rotationBetween(x_dir, pScalar_dir);
+
+    double phi = random->Uniform(0, 2.0 * M_PI);
+    geom3::Rotation3 rand_rot(pScalar_dir, phi);
+
+    double muonMomentum = sqrt(std::pow(mass/2.,2) - std::pow(siren::utilities::Constants::muonMass,2));
+    rk::P4 pMuPlus_ScalarRest(muonMomentum*geom3::Vector3(CosTheta,SinTheta,0),siren::utilities::Constants::muonMass);
+    pMuPlus_ScalarRest.rotate(x_to_pScalar_rot);
+    pMuPlus_ScalarRest.rotate(rand_rot);
+
+    rk::P4 pMuPlus = pMuPlus_ScalarRest.boost(boost_to_lab);
+    rk::P4 pMuMinus(pScalar.momentum() - pMuPlus.momentum(),siren::utilities::Constants::muonMass);
+
+
+    siren::dataclasses::SecondaryParticleRecord & MuPlus = record.GetSecondaryParticleRecord(0);
+    siren::dataclasses::SecondaryParticleRecord & MuMinus = record.GetSecondaryParticleRecord(1);
+
+    assert(MuPlus.type == siren::dataclasses::ParticleType::MuPlus);
+    assert(MuMinus.type == siren::dataclasses::ParticleType::MuMinus);
+
+
+    MuPlus.SetFourMomentum({pMuPlus.e(), pMuPlus.px(), pMuPlus.py(), pMuPlus.pz()});
+    MuPlus.SetMass(pMuPlus.m());
+    MuPlus.SetHelicity(0.5); // right-handed
+
+    MuMinus.SetFourMomentum({pMuMinus.e(), pMuMinus.px(), pMuMinus.py(), pMuMinus.pz()});
+    MuPlus.SetMass(pMuMinus.m());
+    MuPlus.SetHelicity(-0.5); // left-handed
+}
+
+double ScalarDecay::FinalStateProbability(dataclasses::InteractionRecord const & record) const {
+  double dd = DifferentialDecayWidth(record);
+  double td = TotalDecayWidthForFinalState(record);
+  if (dd == 0) return 0.;
+  else if (td == 0) return 0.;
+  else return dd/td;
+}
+
+} // namespace interactions
+} // namespace siren
+

projects/interactions/private/pybindings/ScalarDecay.h

@@ -0,0 +1,35 @@
+#include <set>
+#include <memory>
+#include <string>
+
+#include <pybind11/pybind11.h>
+#include <pybind11/operators.h>
+#include <pybind11/stl.h>
+
+#include "../../public/SIREN/interactions/CrossSection.h"
+#include "../../public/SIREN/interactions/ScalarDecay.h"
+#include "../../../dataclasses/public/SIREN/dataclasses/Particle.h"
+#include "../../../dataclasses/public/SIREN/dataclasses/InteractionRecord.h"
+#include "../../../dataclasses/public/SIREN/dataclasses/InteractionSignature.h"
+#include "../../../geometry/public/SIREN/geometry/Geometry.h"
+#include "../../../utilities/public/SIREN/utilities/Random.h"
+
+void register_ScalarDecay(pybind11::module_ & m) {
+    using namespace pybind11;
+    using namespace siren::interactions;
+
+    class_<ScalarDecay, std::shared_ptr<ScalarDecay>, Decay> scalardecay(m, "ScalarDecay");
+
+    scalardecay
+        .def(init<double, double>())
+        .def(self == self)
+        .def("GetMass",&ScalarDecay::GetMass)
+        .def("TotalDecayWidth",overload_cast<siren::dataclasses::InteractionRecord const &>(&ScalarDecay::TotalDecayWidth, const_))
+        .def("TotalDecayWidth",overload_cast<siren::dataclasses::ParticleType>(&ScalarDecay::TotalDecayWidth, const_))
+        .def("TotalDecayWidthForFinalState",&ScalarDecay::TotalDecayWidthForFinalState)
+        .def("DifferentialDecayWidth",&ScalarDecay::DifferentialDecayWidth)
+        .def("GetPossibleSignatures",&ScalarDecay::GetPossibleSignatures)
+        .def("GetPossibleSignaturesFromParent",&ScalarDecay::GetPossibleSignaturesFromParent)
+        .def("FinalStateProbability",&ScalarDecay::FinalStateProbability)
+        ;
+}

projects/interactions/private/pybindings/interactions.cxx

@@ -9,6 +9,7 @@
 #include "../../public/SIREN/interactions/DipoleFromTable.h"
 #include "../../public/SIREN/interactions/DarkNewsCrossSection.h"
 #include "../../public/SIREN/interactions/DarkNewsDecay.h"
+#include "../../public/SIREN/interactions/ScalarDecay.h"
 
 #include "./CrossSection.h"
 #include "./DipoleFromTable.h"
@@ -20,6 +21,7 @@
 #include "./NeutrissimoDecay.h"
 #include "./InteractionCollection.h"
 #include "./DummyCrossSection.h"
+#include "./ScalarDecay.h"
 
 #include <pybind11/pybind11.h>
 #include <pybind11/operators.h>
@@ -42,4 +44,5 @@ PYBIND11_MODULE(interactions,m) {
     register_NeutrissimoDecay(m);
     register_InteractionCollection(m);
     register_DummyCrossSection(m);
+    register_ScalarDecay(m);
 }

projects/interactions/public/SIREN/interactions/ScalarDecay.h

@@ -0,0 +1,91 @@
+#pragma once
+#ifndef SIREN_ScalarDecay_H
+#define SIREN_ScalarDecay_H
+
+#include <set>                                    // for set
+#include <memory>                                 // for shared_ptr
+#include <string>                                 // for string
+#include <vector>                                 // for vector
+#include <cstdint>                                // for uint32_t
+#include <stdexcept>                              // for runtime_error
+
+#include <cereal/cereal.hpp>
+#include <cereal/archives/json.hpp>
+#include <cereal/archives/binary.hpp>
+#include <cereal/types/map.hpp>
+#include <cereal/types/set.hpp>
+#include <cereal/types/polymorphic.hpp>
+#include <cereal/types/base_class.hpp>
+#include <cereal/types/utility.hpp>
+
+#include "SIREN/interactions/Decay.h"   // for Decay
+#include "SIREN/dataclasses/Particle.h"  // for Particle
+
+namespace siren { namespace dataclasses { class InteractionRecord; } }
+namespace siren { namespace dataclasses { class CrossSectionDistributionRecord; } }
+namespace siren { namespace dataclasses { struct InteractionSignature; } }
+namespace siren { namespace utilities { class SIREN_random; } }
+
+namespace siren {
+namespace interactions {
+
+class ScalarDecay : public Decay {
+friend cereal::access;
+protected:
+ScalarDecay() {};
+private:
+    double mass;
+    double coupling; // coupling
+    const std::set<siren::dataclasses::ParticleType> primary_types = {siren::dataclasses::ParticleType::Scalar};
+public:
+    ScalarDecay(double mass, double coupling) : mass(mass), coupling(coupling) {};
+    virtual bool equal(Decay const & other) const override;
+    double GetMass() const {return mass;};
+    virtual double TotalDecayWidth(dataclasses::InteractionRecord const &) const override;
+    virtual double TotalDecayWidth(siren::dataclasses::ParticleType primary) const override;
+    virtual double TotalDecayWidthForFinalState(dataclasses::InteractionRecord const &) const override;
+    virtual double DifferentialDecayWidth(dataclasses::InteractionRecord const &) const override;
+    virtual void SampleFinalState(dataclasses::CrossSectionDistributionRecord &, std::shared_ptr<siren::utilities::SIREN_random>) const override;
+    virtual std::vector<siren::dataclasses::InteractionSignature> GetPossibleSignatures() const override;
+    virtual std::vector<siren::dataclasses::InteractionSignature> GetPossibleSignaturesFromParent(siren::dataclasses::ParticleType primary) const override;
+    virtual double FinalStateProbability(dataclasses::InteractionRecord const & record) const override;
+public:
+    virtual std::vector<std::string> DensityVariables() const override;
+    template<typename Archive>
+    void save(Archive & archive, std::uint32_t const version) const {
+        if(version == 0) {
+            archive(::cereal::make_nvp("PrimaryTypes", primary_types));
+            archive(::cereal::make_nvp("Mass", mass));
+            archive(::cereal::make_nvp("Coupling", coupling));
+            archive(::cereal::make_nvp("Decay", cereal::virtual_base_class<Decay>(this)));
+        } else {
+            throw std::runtime_error("ScalarDecay only supports version <= 0!");
+        }
+    }
+    template<typename Archive>
+    void load_and_construct(Archive & archive, cereal::construct<ScalarDecay> & construct, std::uint32_t version) {
+        if(version == 0) {
+            std::set<siren::dataclasses::ParticleType> _primary_types;
+            double _mass;
+            double _coupling;
+
+            archive(::cereal::make_nvp("PrimaryTypes", _primary_types));
+            archive(::cereal::make_nvp("Mass", _mass));
+            archive(::cereal::make_nvp("Coupling", _coupling));
+            construct(_mass, _coupling, _primary_types);
+            archive(::cereal::make_nvp("Decay", cereal::virtual_base_class<Decay>(construct.ptr())));
+        } else {
+            throw std::runtime_error("ScalarDecay only supports version <= 0!");
+        }
+    }
+
+};
+
+} // namespace interactions
+} // namespace siren
+
+CEREAL_CLASS_VERSION(siren::interactions::ScalarDecay, 0);
+CEREAL_REGISTER_TYPE(siren::interactions::ScalarDecay);
+CEREAL_REGISTER_POLYMORPHIC_RELATION(siren::interactions::Decay, siren::interactions::ScalarDecay);
+
+#endif // SIREN_ScalarDecay_H

projects/utilities/public/SIREN/utilities/Integration.h

@@ -106,7 +106,7 @@ struct TrapezoidalIntegrator{
 template<typename FuncType>
 double rombergIntegrate(const FuncType& func, double a, double b, double tol=1e-6){
   const unsigned int order=5;
-  const unsigned int maxIter=20;
+  const unsigned int maxIter=100;
   if(tol<0)
      throw(std::runtime_error("Integration tolerance must be positive"));
 

resources/Fluxes/SCALAR_ATMO/SCALAR_ATMO-v1.0/FluxCalculator.py

@@ -0,0 +1,24 @@
+import os
+s_per_year = 60*60*24*365.25
+def MakeFluxFile(tag, abs_flux_dir):
+    '''
+    Accepts the following tags:
+        phi
+    '''
+    if tag not in ["phi"]:
+        print(" tag %s is not valid"%(tag))
+        exit(0)
+    input_flux_file = os.path.join(abs_flux_dir,
+                                   "%s.dat"%tag)
+    output_flux_file = os.path.join(abs_flux_dir,
+                                    "%s.txt"%tag)
+    with open(input_flux_file,"r") as fin:
+        all_lines = fin.readlines()
+        headers = all_lines[0].strip().split()
+        data = [line.strip().split() for line in all_lines[1:]]
+        with open(output_flux_file,"w") as fout:
+            for row in data:
+                energy,flux = float(row[0]),float(row[1])
+                flux = flux*s_per_year*1e4*4*3.14159/(energy**3) # put flux in units of nu/m^2/GeV/yr
+                print(energy,flux,file=fout)
+    return output_flux_file