From abf68839d41dd92b138fa66fad83a9abdec45bd8 Mon Sep 17 00:00:00 2001 From: Miaochen Jin Date: Tue, 19 Nov 2024 15:41:19 -0500 Subject: [PATCH] new stuff --- .../dataclasses/private/InteractionRecord.cxx | 20 + projects/dataclasses/private/Particle.cxx | 6 +- .../public/SIREN/dataclasses/Particle.h | 1 + projects/injection/private/Injector.cxx | 16 +- projects/interactions/CMakeLists.txt | 2 + .../private/CharmDISFromSpline.cxx | 2 + .../private/QuarkDISFromSpline.cxx | 367 ++++++++++++++---- .../private/pybindings/QuarkDISFromSpline.h | 90 +++++ .../private/pybindings/interactions.cxx | 3 + .../SIREN/interactions/QuarkDISFromSpline.h | 32 +- python/SIREN_Controller.py | 1 + 11 files changed, 458 insertions(+), 82 deletions(-) create mode 100644 projects/interactions/private/pybindings/QuarkDISFromSpline.h diff --git a/projects/dataclasses/private/InteractionRecord.cxx b/projects/dataclasses/private/InteractionRecord.cxx index 48c52530..4b79b75b 100644 --- a/projects/dataclasses/private/InteractionRecord.cxx +++ b/projects/dataclasses/private/InteractionRecord.cxx @@ -327,13 +327,22 @@ void PrimaryDistributionRecord::FinalizeAvailable(InteractionRecord & record) co } void PrimaryDistributionRecord::Finalize(InteractionRecord & record) const { + //std::cout << "starting finalize" << std::endl; + //std::cout << record.signature << std::endl; + //std::cout << record.primary_momentum[0] << std::endl; + record.signature.primary_type = type; record.primary_id = GetID(); record.interaction_vertex = GetInteractionVertex(); record.primary_initial_position = GetInitialPosition(); + //std::cout << "in primary record finalize" << std::endl; record.primary_mass = GetMass(); record.primary_momentum = GetFourMomentum(); record.primary_helicity = GetHelicity(); + //std::cout << "finished primary record finalize" << std::endl; + + //std::cout << record.signature << std::endl; + //std::cout << record.primary_momentum[0] << std::endl; } ///////////////////////////////////////// @@ -553,11 +562,16 @@ void SecondaryParticleRecord::UpdateMomentum() const { } void SecondaryParticleRecord::Finalize(InteractionRecord & record) const { + ////std::cout << "SecPartRecord::Finalize : starting for " << type << std::endl; + assert(record.signature.secondary_types.at(secondary_index) == type); record.secondary_ids.at(secondary_index) = GetID(); record.secondary_masses.at(secondary_index) = GetMass(); record.secondary_momenta.at(secondary_index) = GetFourMomentum(); record.secondary_helicities.at(secondary_index) = GetHelicity(); + ////std::cout << "SecPartRecord::Finalize : finished for " << type << " with" << std::endl; + ////std::cout << record << std::endl; + } ///////////////////////////////////////// @@ -676,6 +690,11 @@ SecondaryParticleRecord const & CrossSectionDistributionRecord::GetSecondaryPart } void CrossSectionDistributionRecord::Finalize(InteractionRecord & record) const { + //std::cout << "XsecDistRecord::Finalize : starting" << std::endl; + //std::cout << record.signature << std::endl; + //std::cout << signature << std::endl; + //std::cout << primary_momentum[0] << std::endl; + record.target_id = target_id; record.target_mass = target_mass; record.target_helicity = target_helicity; @@ -688,6 +707,7 @@ void CrossSectionDistributionRecord::Finalize(InteractionRecord & record) const record.secondary_helicities.resize(secondary_particles.size()); for(SecondaryParticleRecord const & secondary : secondary_particles) { + //std::cout << "XsecDistRecord::Finalize : going to secondaries: " << secondary << std::endl; secondary.Finalize(record); } } diff --git a/projects/dataclasses/private/Particle.cxx b/projects/dataclasses/private/Particle.cxx index 7952adc1..afc9e559 100644 --- a/projects/dataclasses/private/Particle.cxx +++ b/projects/dataclasses/private/Particle.cxx @@ -93,6 +93,10 @@ bool isHadron(Particle::ParticleType p){ p==Particle::ParticleType::DPlus || p==Particle::ParticleType::DMinus); } - +bool isD(Particle::ParticleType p){ + return(p==Particle::ParticleType::D0 || p==Particle::ParticleType::D0Bar || + p==Particle::ParticleType::DPlus || p==Particle::ParticleType::DMinus); + } + } // namespace utilities } // namespace siren diff --git a/projects/dataclasses/public/SIREN/dataclasses/Particle.h b/projects/dataclasses/public/SIREN/dataclasses/Particle.h index fff79c4a..a272b8c5 100644 --- a/projects/dataclasses/public/SIREN/dataclasses/Particle.h +++ b/projects/dataclasses/public/SIREN/dataclasses/Particle.h @@ -76,6 +76,7 @@ bool isCharged(ParticleType p); bool isNeutrino(ParticleType p); bool isQuark(Particle::ParticleType p); bool isHadron(Particle::ParticleType p); +bool isD(Particle::ParticleType p); } // namespace dataclasses diff --git a/projects/injection/private/Injector.cxx b/projects/injection/private/Injector.cxx index 7ca32979..197da2c8 100644 --- a/projects/injection/private/Injector.cxx +++ b/projects/injection/private/Injector.cxx @@ -159,7 +159,7 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record throw(siren::utilities::InjectionFailure("No particle interaction!")); } - //std::cout << "in sample cross section" << std::endl; + ////std::cout << "in sample cross section" << std::endl; std::set const & possible_targets = interactions->TargetTypes(); @@ -267,6 +267,7 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record } } } + //std::cout << "injector :: sample cross sections: after obtaining signatures" << std::endl; if(total_prob == 0) throw(siren::utilities::InjectionFailure("No valid interactions for this event!")); // Throw a random number @@ -287,12 +288,17 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record record.target_mass = detector_model->GetTargetMass(record.signature.target_type); siren::dataclasses::CrossSectionDistributionRecord xsec_record(record); if(r <= xsec_prob) { - //std::cout << "going into sampel final state" << std::endl; + //std::cout << "injector::sample cross section: going into sampel final state" << std::endl; matching_cross_sections[index]->SampleFinalState(xsec_record, random); + //std::cout << "injector::sample cross section: finished sampling" << std::endl; + } else { matching_decays[index - matching_cross_sections.size()]->SampleFinalState(xsec_record, random); } + ////std::cout << "injector::sample cross section: calling finalizing" << std::endl; xsec_record.Finalize(record); + ////std::cout << "injector::sample cross section: finished finalizing" << std::endl; + } } @@ -341,13 +347,15 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() { // Initial Process while(true) { tries += 1; + ////std::cout << "injector::GenerateEvent : trying to generate with trial number " << tries << std::endl; try { - //std::cout << "generating primary process" << std::endl; + ////std::cout << "injector::GenerateEvent : generating primary process" << std::endl; siren::dataclasses::PrimaryDistributionRecord primary_record(primary_process->GetPrimaryType()); for(auto & distribution : primary_process->GetPrimaryInjectionDistributions()) { distribution->Sample(random, detector_model, primary_process->GetInteractions(), primary_record); } primary_record.Finalize(record); + //std::cout << "injector::GenerateEvent : sampling cross section" << std::endl; SampleCrossSection(record); break; } catch(siren::utilities::InjectionFailure const & e) { @@ -367,6 +375,7 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() { std::shared_ptr parent = tree.add_entry(record); // Secondary Processes + //std::cout << "injector::GenerateEvent : sampling secondary process" << std::endl; std::deque, std::shared_ptr>> secondaries; std::function)> add_secondaries = [&](std::shared_ptr parent) { for(size_t i=0; irecord.signature.secondary_types.size(); ++i) { @@ -399,6 +408,7 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() { } } + //std::cout << "finished sampling secondary process" << std::endl; injected_events += 1; return tree; } diff --git a/projects/interactions/CMakeLists.txt b/projects/interactions/CMakeLists.txt index a8eec569..ddd9936c 100644 --- a/projects/interactions/CMakeLists.txt +++ b/projects/interactions/CMakeLists.txt @@ -21,6 +21,8 @@ LIST (APPEND interactions_SOURCES ${PROJECT_SOURCE_DIR}/projects/interactions/private/CharmHadronization.cxx ${PROJECT_SOURCE_DIR}/projects/interactions/private/CharmMesonDecay.cxx ${PROJECT_SOURCE_DIR}/projects/interactions/private/DMesonELoss.cxx + ${PROJECT_SOURCE_DIR}/projects/interactions/private/QuarkDISFromSpline.cxx + ) add_library(SIREN_interactions OBJECT ${interactions_SOURCES}) set_property(TARGET SIREN_interactions PROPERTY POSITION_INDEPENDENT_CODE ON) diff --git a/projects/interactions/private/CharmDISFromSpline.cxx b/projects/interactions/private/CharmDISFromSpline.cxx index 9844231e..27921cb3 100644 --- a/projects/interactions/private/CharmDISFromSpline.cxx +++ b/projects/interactions/private/CharmDISFromSpline.cxx @@ -45,6 +45,8 @@ namespace { bool kinematicallyAllowed(double x, double y, double E, double M, double m) { if(x > 1) //Eq. 6 right inequality return false; + // this is to get rid of the infinities as a temporary solution + if (x < 1e-6 or y < 1e-6) return false; if(x < ((m * m) / (2 * M * (E - m)))) //Eq. 6 left inequality return false; //denominator of a and b diff --git a/projects/interactions/private/QuarkDISFromSpline.cxx b/projects/interactions/private/QuarkDISFromSpline.cxx index 9844231e..0c1842e0 100644 --- a/projects/interactions/private/QuarkDISFromSpline.cxx +++ b/projects/interactions/private/QuarkDISFromSpline.cxx @@ -1,4 +1,4 @@ -#include "SIREN/interactions/CharmDISFromSpline.h" +#include "SIREN/interactions/QuarkDISFromSpline.h" #include // for map, opera... #include // for set, opera... @@ -27,6 +27,8 @@ #include "SIREN/dataclasses/Particle.h" // for Particle #include "SIREN/utilities/Random.h" // for SIREN_random #include "SIREN/utilities/Constants.h" // for electronMass +#include "SIREN/utilities/Errors.h" // for PythonImplementationError + namespace siren { namespace interactions { @@ -47,6 +49,8 @@ bool kinematicallyAllowed(double x, double y, double E, double M, double m) { return false; if(x < ((m * m) / (2 * M * (E - m)))) //Eq. 6 left inequality return false; + if (x < 1e-6 || y < 1e-6) return false; + //denominator of a and b double d = 2 * (1 + (M * x) / (2 * E)); //the numerator of a (or a*d) @@ -62,47 +66,63 @@ bool kinematicallyAllowed(double x, double y, double E, double M, double m) { } } -CharmDISFromSpline::CharmDISFromSpline() {} +QuarkDISFromSpline::QuarkDISFromSpline() { + // initialize the pdf normalization and cdf table for the hadronization process + normalize_pdf(); + compute_cdf(); +} -CharmDISFromSpline::CharmDISFromSpline(std::vector differential_data, std::vector total_data, int interaction, double target_mass, double minimum_Q2, std::set primary_types, std::set target_types, std::string units) : primary_types_(primary_types), target_types_(target_types), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) { +QuarkDISFromSpline::QuarkDISFromSpline(std::vector differential_data, std::vector total_data, int interaction, double target_mass, double minimum_Q2, std::set primary_types, std::set target_types, std::string units) : primary_types_(primary_types), target_types_(target_types), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) { + normalize_pdf(); + compute_cdf(); LoadFromMemory(differential_data, total_data); InitializeSignatures(); SetUnits(units); } -CharmDISFromSpline::CharmDISFromSpline(std::vector differential_data, std::vector total_data, int interaction, double target_mass, double minimum_Q2, std::vector primary_types, std::vector target_types, std::string units) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) { +QuarkDISFromSpline::QuarkDISFromSpline(std::vector differential_data, std::vector total_data, int interaction, double target_mass, double minimum_Q2, std::vector primary_types, std::vector target_types, std::string units) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) { + normalize_pdf(); + compute_cdf(); LoadFromMemory(differential_data, total_data); InitializeSignatures(); SetUnits(units); } -CharmDISFromSpline::CharmDISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minimum_Q2, std::set primary_types, std::set target_types, std::string units) : primary_types_(primary_types), target_types_(target_types), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) { +QuarkDISFromSpline::QuarkDISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minimum_Q2, std::set primary_types, std::set target_types, std::string units) : primary_types_(primary_types), target_types_(target_types), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) { + normalize_pdf(); + compute_cdf(); LoadFromFile(differential_filename, total_filename); InitializeSignatures(); SetUnits(units); } -CharmDISFromSpline::CharmDISFromSpline(std::string differential_filename, std::string total_filename, std::set primary_types, std::set target_types, std::string units) : primary_types_(primary_types), target_types_(target_types) { +QuarkDISFromSpline::QuarkDISFromSpline(std::string differential_filename, std::string total_filename, std::set primary_types, std::set target_types, std::string units) : primary_types_(primary_types), target_types_(target_types) { + normalize_pdf(); + compute_cdf(); LoadFromFile(differential_filename, total_filename); ReadParamsFromSplineTable(); InitializeSignatures(); SetUnits(units); } -CharmDISFromSpline::CharmDISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minimum_Q2, std::vector primary_types, std::vector target_types, std::string units) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) { +QuarkDISFromSpline::QuarkDISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minimum_Q2, std::vector primary_types, std::vector target_types, std::string units) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) { + normalize_pdf(); + compute_cdf(); LoadFromFile(differential_filename, total_filename); InitializeSignatures(); SetUnits(units); } -CharmDISFromSpline::CharmDISFromSpline(std::string differential_filename, std::string total_filename, std::vector primary_types, std::vector target_types, std::string units) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()) { +QuarkDISFromSpline::QuarkDISFromSpline(std::string differential_filename, std::string total_filename, std::vector primary_types, std::vector target_types, std::string units) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()) { + normalize_pdf(); + compute_cdf(); LoadFromFile(differential_filename, total_filename); ReadParamsFromSplineTable(); InitializeSignatures(); SetUnits(units); } -void CharmDISFromSpline::SetUnits(std::string units) { +void QuarkDISFromSpline::SetUnits(std::string units) { std::transform(units.begin(), units.end(), units.begin(), [](unsigned char c){ return std::tolower(c); }); if(units == "cm") { @@ -114,13 +134,17 @@ void CharmDISFromSpline::SetUnits(std::string units) { } } -void CharmDISFromSpline::SetInteractionType(int interaction) { +void QuarkDISFromSpline::SetInteractionType(int interaction) { interaction_type_ = interaction; } -bool CharmDISFromSpline::equal(CrossSection const & other) const { - const CharmDISFromSpline* x = dynamic_cast(&other); +void QuarkDISFromSpline::SetQuarkType(int q_type) { + quark_type_ = q_type; +} +bool QuarkDISFromSpline::equal(CrossSection const & other) const { + const QuarkDISFromSpline* x = dynamic_cast(&other); + // to do: include more features in the hadronization side to check equivalence if(!x) return false; else @@ -146,7 +170,7 @@ bool CharmDISFromSpline::equal(CrossSection const & other) const { x->total_cross_section_); } -void CharmDISFromSpline::LoadFromFile(std::string dd_crossSectionFile, std::string total_crossSectionFile) { +void QuarkDISFromSpline::LoadFromFile(std::string dd_crossSectionFile, std::string total_crossSectionFile) { differential_cross_section_ = photospline::splinetable<>(dd_crossSectionFile.c_str()); @@ -161,12 +185,12 @@ void CharmDISFromSpline::LoadFromFile(std::string dd_crossSectionFile, std::stri + " dimensions, should have 1, log10(E)"); } -void CharmDISFromSpline::LoadFromMemory(std::vector & differential_data, std::vector & total_data) { +void QuarkDISFromSpline::LoadFromMemory(std::vector & differential_data, std::vector & total_data) { differential_cross_section_.read_fits_mem(differential_data.data(), differential_data.size()); total_cross_section_.read_fits_mem(total_data.data(), total_data.size()); } -double CharmDISFromSpline::GetLeptonMass(siren::dataclasses::ParticleType lepton_type) { +double QuarkDISFromSpline::GetLeptonMass(siren::dataclasses::ParticleType lepton_type) { int32_t lepton_number = std::abs(static_cast(lepton_type)); double lepton_mass; switch(lepton_number) { @@ -192,7 +216,45 @@ double CharmDISFromSpline::GetLeptonMass(siren::dataclasses::ParticleType lepton return lepton_mass; } -void CharmDISFromSpline::ReadParamsFromSplineTable() { +double QuarkDISFromSpline::getHadronMass(siren::dataclasses::ParticleType hadron_type) { + switch(hadron_type){ + case siren::dataclasses::ParticleType::D0: + return( siren::utilities::Constants::D0Mass); + case siren::dataclasses::ParticleType::D0Bar: + return( siren::utilities::Constants::D0Mass); + case siren::dataclasses::ParticleType::DPlus: + return( siren::utilities::Constants::DPlusMass); + case siren::dataclasses::ParticleType::DMinus: + return( siren::utilities::Constants::DPlusMass); + case siren::dataclasses::ParticleType::Charm: + return( siren::utilities::Constants::CharmMass); + case siren::dataclasses::ParticleType::CharmBar: + return( siren::utilities::Constants::CharmMass); + default: + return(0.0); + } +} + + +std::map QuarkDISFromSpline::getIndices(siren::dataclasses::InteractionSignature signature) { + int lepton_id, hadron_id, meson_id; + for (size_t i = 0; i < signature.secondary_types.size(); i++){ + if (isLepton(signature.secondary_types[i])) { + lepton_id = i; + continue; + } else if (isD(signature.secondary_types[i])) { + meson_id = i; + continue; + } else { + hadron_id = i; + continue; + } + } + return {{"lepton", lepton_id}, {"hadron", hadron_id}, {"meson", meson_id}}; +} + + +void QuarkDISFromSpline::ReadParamsFromSplineTable() { // returns true if successfully read target mass bool mass_good = differential_cross_section_.read_key("TARGETMASS", target_mass_); if (mass_good) {std::cout << "read target mass!!" << std::endl;} // for debugging purposes @@ -200,19 +262,24 @@ void CharmDISFromSpline::ReadParamsFromSplineTable() { bool int_good = differential_cross_section_.read_key("INTERACTION", interaction_type_); // returns true if successfully read minimum Q2 bool q2_good = differential_cross_section_.read_key("Q2MIN", minimum_Q2_); + // returns true if successfully read quark type + bool qtype_good = differential_cross_section_.read_key("QUARKTYPE", quark_type_); + if(!int_good) { // assume DIS to preserve compatability with previous versions interaction_type_ = 1; } + if (!qtype_good) { + quark_type_ = 1; // assume quark is produced + } + if(!q2_good) { // assume 1 GeV^2 minimum_Q2_ = 1; } - std::cout << "Q2 good status is " << q2_good << "and is set to " << minimum_Q2_; - if(!mass_good) { if(int_good) { if(interaction_type_ == 1 or interaction_type_ == 2) { @@ -235,23 +302,19 @@ void CharmDISFromSpline::ReadParamsFromSplineTable() { } } } - std::cout << "target mass is " << target_mass_ << std::endl; - } -void CharmDISFromSpline::InitializeSignatures() { +void QuarkDISFromSpline::InitializeSignatures() { signatures_.clear(); for(auto primary_type : primary_types_) { dataclasses::InteractionSignature signature; signature.primary_type = primary_type; - if(not isNeutrino(primary_type)) { throw std::runtime_error("This DIS implementation only supports neutrinos as primaries!"); } - + // first push back the charged lepton product siren::dataclasses::ParticleType charged_lepton_product = siren::dataclasses::ParticleType::unknown; siren::dataclasses::ParticleType neutral_lepton_product = primary_type; - if(primary_type == siren::dataclasses::ParticleType::NuE) { charged_lepton_product = siren::dataclasses::ParticleType::EMinus; } else if(primary_type == siren::dataclasses::ParticleType::NuEBar) { @@ -267,7 +330,6 @@ void CharmDISFromSpline::InitializeSignatures() { } else { throw std::runtime_error("InitializeSignatures: Unkown parent neutrino type!"); } - if(interaction_type_ == 1) { signature.secondary_types.push_back(charged_lepton_product); } else if(interaction_type_ == 2) { @@ -277,20 +339,100 @@ void CharmDISFromSpline::InitializeSignatures() { } else { throw std::runtime_error("InitializeSignatures: Unkown interaction type!"); } + // now push back the hadron product + signature.secondary_types.push_back(siren::dataclasses::ParticleType::Hadrons); + // define the charmed meson types based on the quark type, now considering only D0 and D+ + if (quark_type_ == 1) { + D_types_ = {siren::dataclasses::Particle::ParticleType::D0, + siren::dataclasses::Particle::ParticleType::DPlus}; + } else { + D_types_ = {siren::dataclasses::Particle::ParticleType::D0Bar, + siren::dataclasses::Particle::ParticleType::DMinus}; + } + // push back the meson type + for (auto meson_type : D_types_) { + dataclasses::InteractionSignature full_signature = signature; + full_signature.secondary_types.push_back(meson_type); + // and finally set the target type and push back the entire signature as well as sig by target + for(auto target_type : target_types_) { + full_signature.target_type = target_type; + + signatures_.push_back(full_signature); + + std::pair key(primary_type, target_type); + signatures_by_parent_types_[key].push_back(full_signature); + } + } + } +} - signature.secondary_types.push_back(siren::dataclasses::ParticleType::Charm); - for(auto target_type : target_types_) { - signature.target_type = target_type; +void QuarkDISFromSpline::normalize_pdf() { + if (fragmentation_integral == 0){ + std::function integrand = [&] (double x) -> double { + return (0.8 / x ) / (std::pow(1 - (1 / x) - (0.2 / (1 - x)), 2)); + }; + fragmentation_integral = siren::utilities::rombergIntegrate(integrand, 0.001, 0.999); + } else { + std::cout << "Something is wrong... you already computed the normalization" << std::endl; + return; + } +} - signatures_.push_back(signature); +void QuarkDISFromSpline::compute_cdf() { + // first set the z nodes + std::vector zspline; + for (int i = 0; i < 100; ++i) { + zspline.push_back(0.01 + i * (0.99-0.01) / 100 ); + } - std::pair key(primary_type, target_type); - signatures_by_parent_types_[key].push_back(signature); + // declare the cdf vectors + std::vector cdf_vector; + std::vector cdf_z_nodes; + std::vector pdf_vector; + + cdf_z_nodes.push_back(0); + cdf_vector.push_back(0); + pdf_vector.push_back(0); + + // compute the spline table + for (int i = 0; i < zspline.size(); ++i) { + if (i == 0) { + double cur_z = zspline[i]; + double cur_pdf = sample_pdf(cur_z); + double area = cur_z * cur_pdf * 0.5; + pdf_vector.push_back(cur_pdf); + cdf_vector.push_back(area); + cdf_z_nodes.push_back(cur_z); + continue; } + double cur_z = zspline[i]; + double cur_pdf = sample_pdf(cur_z); + double area = 0.5 * (pdf_vector[i - 1] + cur_pdf) * (zspline[i] - zspline[i - 1]); + pdf_vector.push_back(cur_pdf); + cdf_z_nodes.push_back(cur_z); + cdf_vector.push_back(area + cdf_vector.back()); } + + cdf_z_nodes.push_back(1); + cdf_vector.push_back(1); + pdf_vector.push_back(0); + + + // set the spline table + siren::utilities::TableData1D inverse_cdf_data; + inverse_cdf_data.x = cdf_vector; + inverse_cdf_data.f = cdf_z_nodes; + + inverseCdfTable = siren::utilities::Interpolator1D(inverse_cdf_data); + + return; +} + +double QuarkDISFromSpline::sample_pdf(double x) const { + return (0.8 / x ) / (std::pow(1 - (1 / x) - (0.2 / (1 - x)), 2)) / fragmentation_integral; } -double CharmDISFromSpline::TotalCrossSection(dataclasses::InteractionRecord const & interaction) const { +double QuarkDISFromSpline::TotalCrossSection(dataclasses::InteractionRecord const & interaction) const { siren::dataclasses::ParticleType primary_type = interaction.signature.primary_type; rk::P4 p1(geom3::Vector3(interaction.primary_momentum[1], interaction.primary_momentum[2], interaction.primary_momentum[3]), interaction.primary_mass); double primary_energy; @@ -303,7 +445,7 @@ double CharmDISFromSpline::TotalCrossSection(dataclasses::InteractionRecord cons return TotalCrossSection(primary_type, primary_energy); } -double CharmDISFromSpline::TotalCrossSection(siren::dataclasses::ParticleType primary_type, double primary_energy) const { +double QuarkDISFromSpline::TotalCrossSection(siren::dataclasses::ParticleType primary_type, double primary_energy) const { if(not primary_types_.count(primary_type)) { throw std::runtime_error("Supplied primary not supported by cross section!"); } @@ -328,34 +470,38 @@ double CharmDISFromSpline::TotalCrossSection(siren::dataclasses::ParticleType pr return unit * std::pow(10.0, log_xs); } -double CharmDISFromSpline::DifferentialCrossSection(dataclasses::InteractionRecord const & interaction) const { +double QuarkDISFromSpline::DifferentialCrossSection(dataclasses::InteractionRecord const & interaction) const { rk::P4 p1(geom3::Vector3(interaction.primary_momentum[1], interaction.primary_momentum[2], interaction.primary_momentum[3]), interaction.primary_mass); rk::P4 p2(geom3::Vector3(0, 0, 0), interaction.target_mass); double primary_energy; primary_energy = interaction.primary_momentum[0]; - assert(interaction.signature.secondary_types.size() == 2); - unsigned int lepton_index = (isLepton(interaction.signature.secondary_types[0])) ? 0 : 1; - unsigned int other_index = 1 - lepton_index; + assert(interaction.signature.secondary_types.size() == 3); + std::map secondaries = getIndices(interaction.signature); + unsigned int lepton_index = secondaries["lepton"]; + unsigned int hadron_index = secondaries["hadron"]; + unsigned int meson_index = secondaries["meson"]; std::array const & mom3 = interaction.secondary_momenta[lepton_index]; - std::array const & mom4 = interaction.secondary_momenta[other_index]; - rk::P4 p3(geom3::Vector3(mom3[1], mom3[2], mom3[3]), interaction.secondary_masses[lepton_index]); - rk::P4 p4(geom3::Vector3(mom4[1], mom4[2], mom4[3]), interaction.secondary_masses[other_index]); + std::array const & mom_x = interaction.secondary_momenta[hadron_index]; + std::array const & mom_d = interaction.secondary_momenta[meson_index]; + rk::P4 p3(geom3::Vector3(mom3[1], mom3[2], mom3[3]), interaction.secondary_masses[lepton_index]); + rk::P4 p_x(geom3::Vector3(mom_x[1], mom_x[2], mom_x[3]), interaction.secondary_masses[hadron_index]); + rk::P4 p_d(geom3::Vector3(mom_d[1], mom_d[2], mom_d[3]), interaction.secondary_masses[meson_index]); + rk::P4 p4 = p_x + p_d; // this assume that we are working in a good frame where the hadronization vertex has 4-momentum conserved rk::P4 q = p1 - p3; + // however p4 is not used in computation here so we should be fine... double Q2 = -q.dot(q); double x, y; double lepton_mass = GetLeptonMass(interaction.signature.secondary_types[lepton_index]); - y = 1.0 - p2.dot(p3) / p2.dot(p1); x = Q2 / (2.0 * p2.dot(q)); double log_energy = log10(primary_energy); std::array coordinates{{log_energy, log10(x), log10(y)}}; std::array centers; - if (Q2 < minimum_Q2_ || !kinematicallyAllowed(x, y, primary_energy, target_mass_, lepton_mass) || !differential_cross_section_.searchcenters(coordinates.data(), centers.data())) { // std::cout << "weighting: revert back to saved x and y" << std::endl; @@ -366,11 +512,9 @@ double CharmDISFromSpline::DifferentialCrossSection(dataclasses::InteractionReco Q2 = 2. * E1_lab * E2_lab * x * y; } return DifferentialCrossSection(primary_energy, x, y, lepton_mass, Q2); - - } -double CharmDISFromSpline::DifferentialCrossSection(double energy, double x, double y, double secondary_lepton_mass, double Q2) const { +double QuarkDISFromSpline::DifferentialCrossSection(double energy, double x, double y, double secondary_lepton_mass, double Q2) const { double log_energy = log10(energy); // check preconditions if(log_energy < differential_cross_section_.lower_extent(0) @@ -386,9 +530,6 @@ double CharmDISFromSpline::DifferentialCrossSection(double energy, double x, dou return 0.0; } - // 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 if(std::isnan(Q2)) { Q2 = 2.0 * energy * target_mass_ * x * y; } @@ -413,22 +554,29 @@ double CharmDISFromSpline::DifferentialCrossSection(double energy, double x, dou std::cout << "energy, x, y, Q2 are " << energy << " " << x << " " << y << " " << Q2 << " " << std::endl; std::cout << "spline value read is " << differential_cross_section_.ndsplineeval(coordinates.data(), centers.data(), 0) << std::endl; } - return unit * result; } -double CharmDISFromSpline::InteractionThreshold(dataclasses::InteractionRecord const & interaction) const { +double QuarkDISFromSpline::InteractionThreshold(dataclasses::InteractionRecord const & interaction) const { // Consider implementing DIS thershold at some point return 0; } -void CharmDISFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionRecord & record, std::shared_ptr random) const { +void QuarkDISFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionRecord & record, std::shared_ptr random) const { + // first obtain the indices from secondaries + // std::cout << "in sample final state" << std::endl; + std::map secondary_indices = getIndices(record.signature); + unsigned int lepton_index = secondary_indices["lepton"]; + unsigned int hadron_index = secondary_indices["hadron"]; + unsigned int meson_index = secondary_indices["meson"]; + // Uses Metropolis-Hastings Algorithm! // useful for cases where we don't know the supremum of our distribution, and the distribution is multi-dimensional if (differential_cross_section_.get_ndim() != 3) { throw std::runtime_error("I expected 3 dimensions in the cross section spline, but got " + std::to_string(differential_cross_section_.get_ndim()) +". Maybe your fits file doesn't have the right 'INTERACTION' key?"); } rk::P4 p1(geom3::Vector3(record.primary_momentum[1], record.primary_momentum[2], record.primary_momentum[3]), record.primary_mass); + // std::cout << "quark::sampleFinalState : primary momentum is read to be " << p1 << std::endl; rk::P4 p2(geom3::Vector3(0, 0, 0), record.target_mass); // we assume that: @@ -444,8 +592,7 @@ void CharmDISFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionR p2_lab = p2; primary_energy = p1_lab.e(); - unsigned int lepton_index = (isLepton(record.signature.secondary_types[0])) ? 0 : 1; - unsigned int other_index = 1 - lepton_index; + // correctly assign lepton, hadron and meson index double m = GetLeptonMass(record.signature.secondary_types[lepton_index]); double m1 = record.primary_mass; @@ -566,9 +713,8 @@ void CharmDISFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionR // std::cout << "trial Q is" << trialQ << std::endl; } } - //////////////////////////////////////////////////////////////////////////////////////////// - //////////////////////////////////////////////////////////////////////////////////////////// - //////////////////////////////////////////////////////////////////////////////////////////// + + // scaling down to handle numerical issues double final_x = pow(10., kin_vars[1]); double final_y = pow(10., kin_vars[2]); record.interaction_parameters.clear(); @@ -646,53 +792,132 @@ void CharmDISFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionR rk::P4 p3; rk::P4 p4; - p3 = p3_lab; - p4 = p4_lab; - + p3 = p3_lab; // now we have our lepton momentum set, which should not be modified from here on + p4 = p4_lab; // momentum of the virtual charm + // std::cout << "charm momentum is " << p4 << std::endl; + + // compute the energy and 3-momentum of the virtual charm + // std::cout << "the virtual charm off-shell mass is " << p4.m() << std::endl; + double p3c = std::sqrt(std::pow(p4.px(), 2) + std::pow(p4.py(), 2) + std::pow(p4.pz(), 2)); + double Ec = p4.e(); //energy of primary charm + double mCH = getHadronMass(record.signature.secondary_types[meson_index]); // obtain charmed hadron mass + + // accept-reject sampling for a valid momentum fragmentation + bool frag_accept; + double randValue; + double z; + double ECH; + + // add a maximum number of trials in the while loop + int max_sampling = 500; + int sampling = 0; + + // sample again if this eenrgy is not kinematically allowed + do { + sampling += 1; + if (sampling > max_sampling) { + std::cout << "energy of the charm is " << Ec << " and momentum is " << p3c << std::endl; + std::cout << "desired mass of hadron is " << mCH << std::endl; + // throw(siren::utilities::InjectionFailure("Failed to sample hadronization!")); + break; + } + randValue = random->Uniform(0,1); + z = inverseCdfTable(randValue); + ECH = z * Ec; + if (std::pow(ECH, 2) - std::pow(mCH, 2) <= 0) { + frag_accept = false; + } else { + frag_accept = true; + } + double test_ED = ECH; + double test_EX = (1-z) * Ec; + double test_p3D = std::sqrt(std::pow(test_ED, 2) - std::pow(mCH, 2)); + double test_rD = test_p3D / p3c; + double test_p3X = std::pow((1 - test_rD), 2) * std::pow(p3c, 2); + if (std::pow(test_EX, 2) - std::pow(test_p3X, 2) <= 0) {frag_accept = false;} else {frag_accept = true;} + } while (!frag_accept); + + // set the 3-momentum of the charmed meson and subsequently the 4-momentum + double p3CH = std::sqrt(std::pow(ECH, 2) - std::pow(mCH, 2)); //obtain charmed hadron 3-momentum + double rCH = p3CH/p3c; // ratio of momentum carried away by the charmed hadron, assume collinearity + rk::P4 p4CH(geom3::Vector3(rCH * record.primary_momentum[1], rCH * record.primary_momentum[2], rCH * record.primary_momentum[3]), mCH); + + // the 4 momentum (and the mass) of the resulting hadronic vertex is determined solely via 4-momentum conservation + rk::P4 p4X = p4 - p4CH; + // let's first assume massless hadron final state + // double EX = (1 - z) * Ec; // energy of the hadronic shower + // double p3X = EX; // assume no hadronic mass + // double rX = p3X/p3c; // assume collinear + // rk::P4 p4X(geom3::Vector3(rX * record.primary_momentum[1], rX * record.primary_momentum[2], rX * record.primary_momentum[3]), 0); + + // now we proceed to saving the final state kinematics std::vector & secondaries = record.GetSecondaryParticleRecords(); siren::dataclasses::SecondaryParticleRecord & lepton = secondaries[lepton_index]; - siren::dataclasses::SecondaryParticleRecord & other = secondaries[other_index]; + siren::dataclasses::SecondaryParticleRecord & hadron = secondaries[hadron_index]; + siren::dataclasses::SecondaryParticleRecord & meson = secondaries[meson_index]; + // std::cout << "QuarkDIS::SampleFInalState : the indices are: " << lepton_index << hadron_index<< meson_index << std::endl; lepton.SetFourMomentum({p3.e(), p3.px(), p3.py(), p3.pz()}); + // std::cout << "setting lepton mass with lepton momentum " << p3 << std::endl; 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); + hadron.SetFourMomentum({p4X.e(), p4X.px(), p4X.py(), p4X.pz()}); + // std::cout << "setting hadron mass with hadron momentum " << p4X << std::endl; + hadron.SetMass(p4X.m()); + hadron.SetHelicity(record.target_helicity); + meson.SetFourMomentum({p4CH.e(), p4CH.px(), p4CH.py(), p4CH.pz()}); + // std::cout << "setting meson mass with meson momentum " << p4CH << std::endl; + meson.SetMass(p4CH.m()); + meson.SetHelicity(record.target_helicity); // this needs working on + // std::cout << "finished sampling final state" << std::endl; +} + +double QuarkDISFromSpline::FragmentationFraction(siren::dataclasses::Particle::ParticleType secondary) const { + if (secondary == siren::dataclasses::Particle::ParticleType::D0 || secondary == siren::dataclasses::Particle::ParticleType::D0Bar) { + return 0.6; + } else if (secondary == siren::dataclasses::Particle::ParticleType::DPlus || secondary == siren::dataclasses::Particle::ParticleType::DMinus) { + return 0.23; + } // D_s and Lambda^+ not yet implemented + return 0; } -double CharmDISFromSpline::FinalStateProbability(dataclasses::InteractionRecord const & interaction) const { +double QuarkDISFromSpline::FinalStateProbability(dataclasses::InteractionRecord const & interaction) const { + // first compute the differential and total cross section double dxs = DifferentialCrossSection(interaction); // if (dxs == 0) { // std::cout << "diff xsec gives 0" << std::endl; // } double txs = TotalCrossSection(interaction); + //then compute the fragmentation probability + std::map secondaries = getIndices(interaction.signature); + unsigned int meson_index = secondaries["meson"]; + double fragfrac = FragmentationFraction(interaction.signature.secondary_types[meson_index]); if(dxs == 0) { return 0.0; } else { // if (txs == 0) {std::cout << "wtf??? txs is 0 in final state prob" << txs << std::endl;} // if (std::isinf(dxs)) {std::cout << "dxs is inf in final state prob" << std::endl;} - return dxs / txs; + return dxs / txs * fragfrac; } } -std::vector CharmDISFromSpline::GetPossiblePrimaries() const { +std::vector QuarkDISFromSpline::GetPossiblePrimaries() const { return std::vector(primary_types_.begin(), primary_types_.end()); } -std::vector CharmDISFromSpline::GetPossibleTargetsFromPrimary(siren::dataclasses::ParticleType primary_type) const { +std::vector QuarkDISFromSpline::GetPossibleTargetsFromPrimary(siren::dataclasses::ParticleType primary_type) const { return std::vector(target_types_.begin(), target_types_.end()); } -std::vector CharmDISFromSpline::GetPossibleSignatures() const { +std::vector QuarkDISFromSpline::GetPossibleSignatures() const { return std::vector(signatures_.begin(), signatures_.end()); } -std::vector CharmDISFromSpline::GetPossibleTargets() const { +std::vector QuarkDISFromSpline::GetPossibleTargets() const { return std::vector(target_types_.begin(), target_types_.end()); } -std::vector CharmDISFromSpline::GetPossibleSignaturesFromParents(siren::dataclasses::ParticleType primary_type, siren::dataclasses::ParticleType target_type) const { +std::vector QuarkDISFromSpline::GetPossibleSignaturesFromParents(siren::dataclasses::ParticleType primary_type, siren::dataclasses::ParticleType target_type) const { std::pair key(primary_type, target_type); if(signatures_by_parent_types_.find(key) != signatures_by_parent_types_.end()) { return signatures_by_parent_types_.at(key); @@ -701,7 +926,7 @@ std::vector CharmDISFromSpline::GetPossibleSi } } -std::vector CharmDISFromSpline::DensityVariables() const { +std::vector QuarkDISFromSpline::DensityVariables() const { return std::vector{"Bjorken x", "Bjorken y"}; } diff --git a/projects/interactions/private/pybindings/QuarkDISFromSpline.h b/projects/interactions/private/pybindings/QuarkDISFromSpline.h new file mode 100644 index 00000000..65e60c44 --- /dev/null +++ b/projects/interactions/private/pybindings/QuarkDISFromSpline.h @@ -0,0 +1,90 @@ +#include +#include +#include + +#include +#include +#include + +#include "../../public/SIREN/interactions/CrossSection.h" +#include "../../public/SIREN/interactions/QuarkDISFromSpline.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_QuarkDISFromSpline(pybind11::module_ & m) { + using namespace pybind11; + using namespace siren::interactions; + + class_, CrossSection> quarkdisfromspline(m, "QuarkDISFromSpline"); + + quarkdisfromspline + + .def(init<>()) + .def(init, std::vector, int, double, double, std::set, std::set, std::string>(), + arg("total_xs_data"), + arg("differential_xs_data"), + arg("interaction"), + arg("target_mass"), + arg("minimum_Q2"), + arg("primary_types"), + arg("target_types"), + arg("units") = std::string("cm")) + .def(init, std::vector, int, double, double, std::vector, std::vector, std::string>(), + arg("total_xs_data"), + arg("differential_xs_data"), + arg("interaction"), + arg("target_mass"), + arg("minimum_Q2"), + arg("primary_types"), + arg("target_types"), + arg("units") = std::string("cm")) + .def(init, std::set, std::string>(), + arg("total_xs_filename"), + arg("differential_xs_filename"), + arg("interaction"), + arg("target_mass"), + arg("minimum_Q2"), + arg("primary_types"), + arg("target_types"), + arg("units") = std::string("cm")) + .def(init, std::set, std::string>(), + arg("total_xs_filename"), + arg("differential_xs_filename"), + arg("primary_types"), + arg("target_types"), + arg("units") = std::string("cm")) + .def(init, std::vector, std::string>(), + arg("total_xs_filename"), + arg("differential_xs_filename"), + arg("interaction"), + arg("target_mass"), + arg("minimum_Q2"), + arg("primary_types"), + arg("target_types"), + arg("units") = std::string("cm")) + .def(init, std::vector, std::string>(), + arg("total_xs_filename"), + arg("differential_xs_filename"), + arg("primary_types"), + arg("target_types"), + arg("units") = std::string("cm")) + .def(self == self) + .def("SetInteractionType",&QuarkDISFromSpline::SetInteractionType) + .def("SetQuarkType",&QuarkDISFromSpline::SetQuarkType) + .def("TotalCrossSection",overload_cast(&QuarkDISFromSpline::TotalCrossSection, const_)) + .def("TotalCrossSection",overload_cast(&QuarkDISFromSpline::TotalCrossSection, const_)) + .def("DifferentialCrossSection",overload_cast(&QuarkDISFromSpline::DifferentialCrossSection, const_)) + .def("DifferentialCrossSection",overload_cast(&QuarkDISFromSpline::DifferentialCrossSection, const_)) + .def("InteractionThreshold",&QuarkDISFromSpline::InteractionThreshold) + .def("FragmentationFraction",&QuarkDISFromSpline::FragmentationFraction) + .def("GetPossibleTargets",&QuarkDISFromSpline::GetPossibleTargets) + .def("GetPossibleTargetsFromPrimary",&QuarkDISFromSpline::GetPossibleTargetsFromPrimary) + .def("GetPossiblePrimaries",&QuarkDISFromSpline::GetPossiblePrimaries) + .def("GetPossibleSignatures",&QuarkDISFromSpline::GetPossibleSignatures) + .def("GetPossibleSignaturesFromParents",&QuarkDISFromSpline::GetPossibleSignaturesFromParents) + .def("FinalStateProbability",&QuarkDISFromSpline::FinalStateProbability); +} + diff --git a/projects/interactions/private/pybindings/interactions.cxx b/projects/interactions/private/pybindings/interactions.cxx index 7cfc91cd..9a2c628e 100644 --- a/projects/interactions/private/pybindings/interactions.cxx +++ b/projects/interactions/private/pybindings/interactions.cxx @@ -6,6 +6,7 @@ #include "../../public/SIREN/interactions/InteractionCollection.h" #include "../../public/SIREN/interactions/DISFromSpline.h" #include "../../public/SIREN/interactions/CharmDISFromSpline.h" +#include "../../public/SIREN/interactions/QuarkDISFromSpline.h" #include "../../public/SIREN/interactions/HNLFromSpline.h" #include "../../public/SIREN/interactions/DipoleFromTable.h" #include "../../public/SIREN/interactions/DarkNewsCrossSection.h" @@ -21,6 +22,7 @@ #include "./DarkNewsDecay.h" #include "./DISFromSpline.h" #include "./CharmDISFromSpline.h" +#include "./QuarkDISFromSpline.h" #include "./HNLFromSpline.h" #include "./Decay.h" #include "./NeutrissimoDecay.h" @@ -55,6 +57,7 @@ PYBIND11_MODULE(interactions,m) { register_DarkNewsDecay(m); register_DISFromSpline(m); register_CharmDISFromSpline(m); + register_QuarkDISFromSpline(m); register_HNLFromSpline(m); register_NeutrissimoDecay(m); register_InteractionCollection(m); diff --git a/projects/interactions/public/SIREN/interactions/QuarkDISFromSpline.h b/projects/interactions/public/SIREN/interactions/QuarkDISFromSpline.h index 62dc88a7..f17139bf 100644 --- a/projects/interactions/public/SIREN/interactions/QuarkDISFromSpline.h +++ b/projects/interactions/public/SIREN/interactions/QuarkDISFromSpline.h @@ -27,6 +27,8 @@ #include "SIREN/interactions/CrossSection.h" // for CrossSe... #include "SIREN/dataclasses/InteractionSignature.h" // for Interac... #include "SIREN/dataclasses/Particle.h" // for Particle +#include "SIREN/utilities/Interpolator.h" +#include "SIREN/utilities/Integration.h" namespace siren { namespace dataclasses { class InteractionRecord; } } namespace siren { namespace utilities { class SIREN_random; } } @@ -45,10 +47,18 @@ friend cereal::access; std::set target_types_; std::map> targets_by_primary_types_; std::map, std::vector> signatures_by_parent_types_; - + std::set D_types_; + + // used by the DIS process int interaction_type_; + int quark_type_; double target_mass_; double minimum_Q2_; + + // used by the hadronization process + double fragmentation_integral = 0; // for storing the integrated unnormed pdf + void normalize_pdf(); // for normalizing pdf and integral, to be called at initialization + siren::utilities::Interpolator1D inverseCdfTable; // for storing the CDF-1 table for the hadronization double unit; @@ -62,36 +72,44 @@ friend cereal::access; QuarkDISFromSpline(std::string differential_filename, std::string total_filename, std::vector primary_types, std::vector target_types, std::string units = "cm"); void SetUnits(std::string units); - // this might be integrated later? could also make another initialization method - // problem with current implementation is that EM is not supported b/c at initialization we assume int = 1 - // this sets the isoscalar target mass void SetInteractionType(int interaction); + void SetQuarkType(int q_type); virtual bool equal(CrossSection const & other) const override; + // function definitions needed to compute the DIS vertex double TotalCrossSection(dataclasses::InteractionRecord const &) const override; double TotalCrossSection(siren::dataclasses::ParticleType primary, double energy) const; double DifferentialCrossSection(dataclasses::InteractionRecord const &) const override; double DifferentialCrossSection(double energy, double x, double y, double secondary_lepton_mass, double Q2=std::numeric_limits::quiet_NaN()) const; double InteractionThreshold(dataclasses::InteractionRecord const &) const override; - void SampleFinalState(dataclasses::CrossSectionDistributionRecord &, std::shared_ptr random) const override; + // function definitions needed to compute the hadronization vertex + double FragmentationFraction(siren::dataclasses::Particle::ParticleType secondary) const; + double sample_pdf(double z) const; + void compute_cdf(); + static double getHadronMass(siren::dataclasses::ParticleType hadron_type); + + // used for both processes + void SampleFinalState(dataclasses::CrossSectionDistributionRecord &, std::shared_ptr random) const override; std::vector GetPossibleTargets() const override; std::vector GetPossibleTargetsFromPrimary(siren::dataclasses::ParticleType primary_type) const override; std::vector GetPossiblePrimaries() const override; std::vector GetPossibleSignatures() const override; std::vector GetPossibleSignaturesFromParents(siren::dataclasses::ParticleType primary_type, siren::dataclasses::ParticleType target_type) const override; - virtual double FinalStateProbability(dataclasses::InteractionRecord const & record) const override; + // other utility functions void LoadFromFile(std::string differential_filename, std::string total_filename); void LoadFromMemory(std::vector & differential_data, std::vector & total_data); + // utilities for DIS parametrs double GetMinimumQ2() const {return minimum_Q2_;}; double GetTargetMass() const {return target_mass_;}; int GetInteractionType() const {return interaction_type_;}; - static double GetLeptonMass(siren::dataclasses::ParticleType lepton_type); + static std::map getIndices(siren::dataclasses::InteractionSignature signature); + public: virtual std::vector DensityVariables() const override; diff --git a/python/SIREN_Controller.py b/python/SIREN_Controller.py index 7daf7721..dc725e2e 100644 --- a/python/SIREN_Controller.py +++ b/python/SIREN_Controller.py @@ -489,6 +489,7 @@ def GenerateEvents(self, N=None, fill_tables_at_exit=True): self.global_times.append(t-self.global_start) prev_time = t count += 1 + # print("finished generating one events") if hasattr(self, "DN_processes"): self.DN_processes.SaveCrossSectionTables(fill_tables_at_exit=fill_tables_at_exit) return self.events