diff --git a/projects/injection/CMakeLists.txt b/projects/injection/CMakeLists.txt index a1838824..e0956aee 100644 --- a/projects/injection/CMakeLists.txt +++ b/projects/injection/CMakeLists.txt @@ -4,12 +4,7 @@ LIST (APPEND injection_SOURCES ${PROJECT_SOURCE_DIR}/projects/injection/private/Process.cxx ${PROJECT_SOURCE_DIR}/projects/injection/private/Injector.cxx ${PROJECT_SOURCE_DIR}/projects/injection/private/WeightingUtils.cxx - ${PROJECT_SOURCE_DIR}/projects/injection/private/TreeWeighter.cxx ${PROJECT_SOURCE_DIR}/projects/injection/private/Weighter.cxx - #${PROJECT_SOURCE_DIR}/projects/injection/private/ColumnDepthSIREN.cxx - #${PROJECT_SOURCE_DIR}/projects/injection/private/CylinderVolumeSIREN.cxx - #${PROJECT_SOURCE_DIR}/projects/injection/private/DecayRangeSIREN.cxx - #${PROJECT_SOURCE_DIR}/projects/injection/private/RangedSIREN.cxx ) add_library(SIREN_injection OBJECT ${injection_SOURCES}) set_property(TARGET SIREN_injection PROPERTY POSITION_INDEPENDENT_CODE ON) diff --git a/projects/injection/private/TreeWeighter.cxx b/projects/injection/private/TreeWeighter.cxx deleted file mode 100644 index da692f3a..00000000 --- a/projects/injection/private/TreeWeighter.cxx +++ /dev/null @@ -1,166 +0,0 @@ -#include "SIREN/injection/TreeWeighter.h" - -#include // for ite... -#include // for array -#include // for assert -#include // for exp -#include // for ini... -#include // for ope... -#include // for set -#include -#include -#include -#include -#include - // for out... -#include "SIREN/interactions/CrossSection.h" // for Cro... -#include "SIREN/interactions/InteractionCollection.h" // for Cro... -#include "SIREN/dataclasses/InteractionRecord.h" // for Int... -#include "SIREN/dataclasses/InteractionSignature.h" // for Int... -#include "SIREN/detector/DetectorModel.h" // for Ear... -#include "SIREN/detector/Coordinates.h" -#include "SIREN/distributions/Distributions.h" // for Inj... -#include "SIREN/geometry/Geometry.h" // for Geo... -#include "SIREN/injection/Injector.h" // for Inj... -#include "SIREN/injection/Process.h" // for Phy... -#include "SIREN/injection/WeightingUtils.h" // for Cro... -#include "SIREN/math/Vector3D.h" // for Vec... - - -#include "SIREN/injection/Injector.h" - -#include "SIREN/distributions/primary/vertex/VertexPositionDistribution.h" - -#include "SIREN/interactions/CrossSection.h" -#include "SIREN/interactions/InteractionCollection.h" - -#include "SIREN/dataclasses/InteractionSignature.h" - -#include - -namespace siren { -namespace injection { - -using detector::DetectorPosition; -using detector::DetectorDirection; - -//--------------- -// class TreeWeighter -//--------------- - -void TreeWeighter::Initialize() { - int i = 0; - primary_process_weighters.reserve(injectors.size()); - secondary_process_weighter_maps.reserve(injectors.size()); - for(auto const & injector : injectors) { - assert(primary_physical_process->MatchesHead(injector->GetPrimaryProcess())); - primary_process_weighters.push_back(std::make_shared(PrimaryProcessWeighter(primary_physical_process, injector->GetPrimaryProcess(), detector_model))); - std::map> - injector_sec_process_weighter_map; - std::map> - injector_sec_process_map = injector->GetSecondaryProcessMap(); - for(auto const & sec_phys_process : secondary_physical_processes) { - try{ - std::shared_ptr sec_inj_process = injector_sec_process_map.at(sec_phys_process->GetPrimaryType()); - assert(sec_phys_process->MatchesHead(sec_inj_process)); // make sure cross section collection matches - injector_sec_process_weighter_map[sec_phys_process->GetPrimaryType()] = - std::make_shared( - SecondaryProcessWeighter( - sec_phys_process, - sec_inj_process,detector_model - ) - ); - } catch(const std::out_of_range& oor) { - std::cout << "Out of Range error: " << oor.what() << '\n'; - std::cout << "Initialization Incomplete: Particle " << sec_phys_process->GetPrimaryType() << " does not exist in injector\n"; - return; - } - } - if(injector_sec_process_weighter_map.size() != injector_sec_process_map.size()) { - std::cout << "Initialization Incomplete: No one-to-one mapping between injection and physical distributions for injector " << i << "\n"; - return; - } - secondary_process_weighter_maps.push_back(injector_sec_process_weighter_map); - } -} - -double TreeWeighter::EventWeight(siren::dataclasses::InteractionTree const & tree) const { - // The weight is given by - // - // [sum_{injectors i} - // x prod_{tree datum d} - // x (prod_{generation dist j} p_gen^{idj}) - // / (prod_{physical dist j} p_phys^{idj}) ] ^-1 - // - double inv_weight = 0; - for(unsigned int idx = 0; idx < injectors.size(); ++idx) { - double physical_probability = 1.0; - double generation_probability = injectors[idx]->EventsToInject();//GenerationProbability(tree); - for(auto const & datum : tree.tree) { - std::tuple bounds; - if(datum->depth() == 0) { - bounds = injectors[idx]->PrimaryInjectionBounds(datum->record); - physical_probability *= primary_process_weighters[idx]->PhysicalProbability(bounds, datum->record); - generation_probability *= primary_process_weighters[idx]->GenerationProbability(*datum); - } - else { - try { - bounds = injectors[idx]->SecondaryInjectionBounds(datum->record); - double phys_prob = secondary_process_weighter_maps[idx].at(datum->record.signature.primary_type)->PhysicalProbability(bounds, datum->record); - double gen_prob = secondary_process_weighter_maps[idx].at(datum->record.signature.primary_type)->GenerationProbability(*datum); - physical_probability *= phys_prob; - generation_probability *= gen_prob; - } catch(const std::out_of_range& oor) { - std::cout << "Out of Range error: " << oor.what() << '\n'; - return 0; - } - } - } - inv_weight += generation_probability / physical_probability; - } - return 1./inv_weight; -} - -void TreeWeighter::SaveWeighter(std::string const & filename) const { - std::ofstream os(filename+".siren_weighter", std::ios::binary); - ::cereal::BinaryOutputArchive archive(os); - this->save(archive,0); -} - -void TreeWeighter::LoadWeighter(std::string const & filename) { - std::cout << "Weighter loading not yet supported... sorry!\n"; - exit(0); - std::ifstream is(filename+".siren_weighter", std::ios::binary); - ::cereal::BinaryInputArchive archive(is); - //this->load(archive,0); -} - -TreeWeighter::TreeWeighter(std::vector> injectors, std::shared_ptr detector_model, std::shared_ptr primary_physical_process, std::vector> secondary_physical_processes) - : injectors(injectors) - , detector_model(detector_model) - , primary_physical_process(primary_physical_process) - , secondary_physical_processes(secondary_physical_processes) -{ - Initialize(); -} - -TreeWeighter::TreeWeighter(std::vector> injectors, std::shared_ptr detector_model, std::shared_ptr primary_physical_process) - : injectors(injectors) - , detector_model(detector_model) - , primary_physical_process(primary_physical_process) - , secondary_physical_processes(std::vector>()) -{ - Initialize(); -} - -TreeWeighter::TreeWeighter(std::vector> _injectors, std::string filename) { - LoadWeighter(filename); - if(_injectors.size() > 0) { - // overwrite the serialized injectors if the user have provided any - injectors = _injectors; - } - Initialize(); -} - -} // namespace injection -} // namespace siren diff --git a/projects/injection/private/Weighter.cxx b/projects/injection/private/Weighter.cxx index 9088e332..81dbe6df 100644 --- a/projects/injection/private/Weighter.cxx +++ b/projects/injection/private/Weighter.cxx @@ -1,29 +1,42 @@ #include "SIREN/injection/Weighter.h" -#include -#include -#include +#include // for ite... +#include // for array +#include // for assert +#include // for exp +#include // for ini... +#include // for ope... +#include // for set +#include #include #include -#include -#include +#include #include -#include -#include + // for out... +#include "SIREN/interactions/CrossSection.h" // for Cro... +#include "SIREN/interactions/InteractionCollection.h" // for Cro... +#include "SIREN/dataclasses/InteractionRecord.h" // for Int... +#include "SIREN/dataclasses/InteractionSignature.h" // for Int... +#include "SIREN/detector/DetectorModel.h" // for Ear... +#include "SIREN/detector/Coordinates.h" +#include "SIREN/distributions/Distributions.h" // for Inj... +#include "SIREN/geometry/Geometry.h" // for Geo... +#include "SIREN/injection/Injector.h" // for Inj... +#include "SIREN/injection/Process.h" // for Phy... +#include "SIREN/injection/WeightingUtils.h" // for Cro... +#include "SIREN/math/Vector3D.h" // for Vec... + + +#include "SIREN/injection/Injector.h" + +#include "SIREN/distributions/primary/vertex/VertexPositionDistribution.h" #include "SIREN/interactions/CrossSection.h" #include "SIREN/interactions/InteractionCollection.h" -#include "SIREN/dataclasses/InteractionRecord.h" + #include "SIREN/dataclasses/InteractionSignature.h" -#include "SIREN/dataclasses/Particle.h" -#include "SIREN/detector/DetectorModel.h" -#include "SIREN/detector/Coordinates.h" -#include "SIREN/distributions/Distributions.h" -#include "SIREN/distributions/primary/vertex/VertexPositionDistribution.h" -#include "SIREN/geometry/Geometry.h" -#include "SIREN/injection/Injector.h" -#include "SIREN/injection/WeightingUtils.h" -#include "SIREN/math/Vector3D.h" + +#include namespace siren { namespace injection { @@ -31,482 +44,54 @@ namespace injection { using detector::DetectorPosition; using detector::DetectorDirection; -namespace { - template - typename std::iterator_traits::value_type accumulate(InIt begin, InIt end) { - typedef typename std::iterator_traits::value_type real; - real sum = real(0); - real running_error = real(0); - real temp; - real difference; - - for (; begin != end; ++begin) { - difference = *begin; - difference -= running_error; - temp = sum; - temp += difference; - running_error = temp; - running_error -= sum; - running_error -= difference; - sum = std::move(temp); - } - return sum; - } - - template - T accumulate(std::initializer_list list) { - return accumulate(list.begin(), list.end()); - } - - double one_minus_exp_of_negative(double x) { - if(x < 1e-1) { - return std::exp(std::log(x) - x/2.0 + x*x/24.0 - x*x*x*x/2880.0); - } else { - return 1.0 - std::exp(-x); - } - } - double log_one_minus_exp_of_negative(double x) { - if(x < 1e-1) { - return std::log(x) - x/2.0 + x*x/24.0 - x*x*x*x/2880.0; - } else if(x > 3) { - double ex = std::exp(-x); - double ex2 = ex * ex; - double ex3 = ex2 * ex; - double ex4 = ex3 * ex; - double ex5 = ex4 * ex; - double ex6 = ex5 * ex; - return -(ex + ex2 / 2.0 + ex3 / 3.0 + ex4 / 4.0 + ex5 / 5.0 + ex6 / 6.0); - } else { - return std::log(1.0 - std::exp(-x)); - } - } -} - //--------------- -// class LeptonWeighter +// class Weighter //--------------- -double LeptonWeighter::InteractionProbability(std::shared_ptr injector, siren::dataclasses::InteractionRecord const & record) const { - std::tuple bounds = injector->PrimaryInjectionBounds(record); - return InteractionProbability(bounds, record); -} - -double LeptonWeighter::InteractionProbability(std::tuple bounds, siren::dataclasses::InteractionRecord const & record) const { - siren::math::Vector3D interaction_vertex( - record.interaction_vertex[0], - record.interaction_vertex[1], - record.interaction_vertex[2]); - - siren::math::Vector3D primary_direction( - record.primary_momentum[1], - record.primary_momentum[2], - record.primary_momentum[3]); - primary_direction.normalize(); - - siren::geometry::Geometry::IntersectionList intersections = detector_model->GetIntersections(DetectorPosition(interaction_vertex), DetectorDirection(primary_direction)); - std::map>> const & cross_sections_by_target = interactions->GetCrossSectionsByTarget(); - std::vector targets; - targets.reserve(cross_sections_by_target.size()); - std::vector total_cross_sections; - double total_decay_length = interactions->TotalDecayLength(record); - siren::dataclasses::InteractionRecord fake_record = record; - for(auto const & target_xs : cross_sections_by_target) { - targets.push_back(target_xs.first); - fake_record.target_mass = detector_model->GetTargetMass(target_xs.first); - std::vector> const & xs_list = target_xs.second; - double total_xs = 0.0; - for(auto const & xs : xs_list) { - std::vector signatures = xs->GetPossibleSignaturesFromParents(record.signature.primary_type, target_xs.first); - for(auto const & signature : signatures) { - fake_record.signature = signature; - // Add total cross section - total_xs += xs->TotalCrossSection(fake_record); - } - } - total_cross_sections.push_back(total_xs); - } - - double total_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(std::get<0>(bounds)), DetectorPosition(std::get<1>(bounds)), targets, total_cross_sections, total_decay_length); - double interaction_probability; - if(total_interaction_depth < 1e-6) { - interaction_probability = total_interaction_depth; - } else { - interaction_probability = one_minus_exp_of_negative(total_interaction_depth); - } - return interaction_probability; -} - -double LeptonWeighter::UnnormalizedPositionProbability(std::shared_ptr injector, siren::dataclasses::InteractionRecord const & record) const { - std::tuple bounds = injector->PrimaryInjectionBounds(record); - return UnnormalizedPositionProbability(bounds, record); -} - -double LeptonWeighter::UnnormalizedPositionProbability(std::tuple bounds, siren::dataclasses::InteractionRecord const & record) const { - siren::math::Vector3D interaction_vertex( - record.interaction_vertex[0], - record.interaction_vertex[1], - record.interaction_vertex[2]); - - siren::math::Vector3D primary_direction( - record.primary_momentum[1], - record.primary_momentum[2], - record.primary_momentum[3]); - primary_direction.normalize(); - - siren::geometry::Geometry::IntersectionList intersections = detector_model->GetIntersections(DetectorPosition(interaction_vertex), DetectorDirection(primary_direction)); - std::map>> const & cross_sections_by_target = interactions->GetCrossSectionsByTarget(); - - unsigned int n_targets = cross_sections_by_target.size(); - - std::vector targets; targets.reserve(n_targets); - std::vector total_cross_sections; - double total_decay_length = interactions->TotalDecayLength(record); - siren::dataclasses::InteractionRecord fake_record = record; - for(auto const & target_xs : cross_sections_by_target) { - targets.push_back(target_xs.first); - fake_record.target_mass = detector_model->GetTargetMass(target_xs.first); - std::vector> const & xs_list = target_xs.second; - double total_xs = 0.0; - for(auto const & xs : xs_list) { - std::vector signatures = xs->GetPossibleSignaturesFromParents(record.signature.primary_type, target_xs.first); - for(auto const & signature : signatures) { - fake_record.signature = signature; - // Add total cross section - total_xs += xs->TotalCrossSection(fake_record); +void Weighter::Initialize() { + int i = 0; + primary_process_weighters.reserve(injectors.size()); + secondary_process_weighter_maps.reserve(injectors.size()); + for(auto const & injector : injectors) { + assert(primary_physical_process->MatchesHead(injector->GetPrimaryProcess())); + primary_process_weighters.push_back(std::make_shared(PrimaryProcessWeighter(primary_physical_process, injector->GetPrimaryProcess(), detector_model))); + std::map> + injector_sec_process_weighter_map; + std::map> + injector_sec_process_map = injector->GetSecondaryProcessMap(); + for(auto const & sec_phys_process : secondary_physical_processes) { + try{ + std::shared_ptr sec_inj_process = injector_sec_process_map.at(sec_phys_process->GetPrimaryType()); + assert(sec_phys_process->MatchesHead(sec_inj_process)); // make sure cross section collection matches + injector_sec_process_weighter_map[sec_phys_process->GetPrimaryType()] = + std::make_shared( + SecondaryProcessWeighter( + sec_phys_process, + sec_inj_process,detector_model + ) + ); + } catch(const std::out_of_range& oor) { + std::cout << "Out of Range error: " << oor.what() << '\n'; + std::cout << "Initialization Incomplete: Particle " << sec_phys_process->GetPrimaryType() << " does not exist in injector\n"; + return; } } - total_cross_sections.push_back(total_xs); - } - - double total_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(std::get<0>(bounds)), DetectorPosition(std::get<1>(bounds)), targets, total_cross_sections, total_decay_length); - double traversed_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(std::get<0>(bounds)), DetectorPosition(interaction_vertex), targets, total_cross_sections, total_decay_length); - double interaction_density = detector_model->GetInteractionDensity(intersections, DetectorPosition(interaction_vertex), targets, total_cross_sections, total_decay_length); - - double prob_density; - if(total_interaction_depth < 1e-6) { - prob_density = interaction_density; - } else { - prob_density = interaction_density * exp(-traversed_interaction_depth); - } - - return prob_density; -} - -double LeptonWeighter::NormalizedPositionProbability(std::tuple bounds, siren::dataclasses::InteractionRecord const & record) const { - siren::math::Vector3D interaction_vertex( - record.interaction_vertex[0], - record.interaction_vertex[1], - record.interaction_vertex[2]); - - siren::math::Vector3D primary_direction( - record.primary_momentum[1], - record.primary_momentum[2], - record.primary_momentum[3]); - primary_direction.normalize(); - - siren::geometry::Geometry::IntersectionList intersections = detector_model->GetIntersections(DetectorPosition(interaction_vertex), DetectorDirection(primary_direction)); - std::map>> const & cross_sections_by_target = interactions->GetCrossSectionsByTarget(); - - unsigned int n_targets = cross_sections_by_target.size(); - - std::vector targets; targets.reserve(n_targets); - std::vector total_cross_sections; - double total_decay_length = interactions->TotalDecayLength(record); - siren::dataclasses::InteractionRecord fake_record = record; - for(auto const & target_xs : cross_sections_by_target) { - targets.push_back(target_xs.first); - fake_record.target_mass = detector_model->GetTargetMass(target_xs.first); - std::vector> const & xs_list = target_xs.second; - double total_xs = 0.0; - for(auto const & xs : xs_list) { - std::vector signatures = xs->GetPossibleSignaturesFromParents(record.signature.primary_type, target_xs.first); - for(auto const & signature : signatures) { - fake_record.signature = signature; - // Add total cross section - total_xs += xs->TotalCrossSection(fake_record); - } + if(injector_sec_process_weighter_map.size() != injector_sec_process_map.size()) { + std::cout << "Initialization Incomplete: No one-to-one mapping between injection and physical distributions for injector " << i << "\n"; + return; } - total_cross_sections.push_back(total_xs); - } - - double total_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(std::get<0>(bounds)), DetectorPosition(std::get<1>(bounds)), targets, total_cross_sections, total_decay_length); - double traversed_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(std::get<0>(bounds)), DetectorPosition(interaction_vertex), targets, total_cross_sections, total_decay_length); - double interaction_density = detector_model->GetInteractionDensity(intersections, DetectorPosition(interaction_vertex), targets, total_cross_sections, total_decay_length); - - double prob_density; - if(total_interaction_depth < 1e-6) { - prob_density = interaction_density / total_interaction_depth; - } else { - prob_density = interaction_density * exp(-log_one_minus_exp_of_negative(total_interaction_depth) - traversed_interaction_depth); + secondary_process_weighter_maps.push_back(injector_sec_process_weighter_map); } - - return prob_density; } -void LeptonWeighter::Initialize() { - // Clear distributions - unique_distributions.clear(); - common_gen_idxs.clear(); - common_phys_idxs.clear(); - distinct_gen_idxs_by_injector.clear(); - distinct_physical_idxs_by_injector.clear(); - unique_contexts.clear(); - context_idx_by_injector.clear(); - normalization = 1.0; - - // Weights are is given by - // w = (\sum_i (\prod_j p_gen^ij / p_phys^ij) )^-1 - // We first want to determine which pairs of p_gen^ij and p_phys^ij cancel - // Secondly we want to determine which p_gen^j are common across all injectors {i} - // and similarly which p_phys^j are common across all injectors {i} - // The calculation can then be simplified by not computing terms that cancel, - // pulling out common terms, and finding duplicate terms - - // To do this we will track unique terms in each ratio - // Initially we assume all term are unique - - // Initialize the state for physical distributions - // true ==> distribution does not cancel and is not common - std::vector>> physical_init_state; - for(auto physical_dist : physical_distributions) { - physical_init_state.push_back(std::make_pair(true, physical_dist)); - const siren::distributions::PhysicallyNormalizedDistribution* p = dynamic_cast(physical_dist.get()); - if(p) { - if(p->IsNormalizationSet()) { - normalization *= p->GetNormalization(); - } - } - } - std::vector>>> physical_distribution_state(injectors.size(), physical_init_state); - assert(physical_distribution_state.size() == injectors.size()); - // Initialize the state for generation distributions - // true ==> distribution does not cancel and is not common - std::vector>>> generation_distribution_state; - generation_distribution_state.reserve(injectors.size()); - for(auto injector : injectors) { - std::vector> dists = injector->GetPrimaryInjectionDistributions(); - std::vector>> dist_state; - dist_state.reserve(dists.size()); - for(auto dist : dists) { - dist_state.push_back(std::make_pair(true, dist)); - } - generation_distribution_state.push_back(dist_state); - } - assert(generation_distribution_state.size() == injectors.size()); - - // Now we can try to identify term that cancel - for(unsigned int i=0; i>> & phys_dists = physical_distribution_state[i]; - std::vector>> & gen_dists = generation_distribution_state[i]; - // Must check every pair of physical and injection distribution (unless already cancelled) - for(unsigned int phys_idx=0; phys_idx> & phys_dist = phys_dists[phys_idx]; - if(not phys_dist.first) // Skip if already cancelled - continue; - for(unsigned int gen_idx=0; gen_idx> & gen_dist = gen_dists[gen_idx]; - if(not gen_dist.first) { // Skip if already cancelled - continue; - } - // Check if dists are equivalent - // Must consider the DetectorModel and InteractionCollection context in the comparison - std::shared_ptr gen_dist_ptr(gen_dist.second); - bool equivalent_dists = - phys_dist.second->AreEquivalent( // physical dist - detector_model, // physical context - interactions, // physical context - gen_dist_ptr, // generation dist - injectors[i]->GetDetectorModel(), // generation context - injectors[i]->GetInteractions()); // generation context - if(not equivalent_dists) { - continue; - } - phys_dist.first = false; - gen_dist.first = false; - break; // This physical dist is cancelled out so we can skip additional comparisons - } - } - } - - // With cancelled terms marked, we can now collect distributions that are common across all terms - // The one exception to this is vertex position distributions - // Physical vertex position distributions depend on the injection bounds and so cannot be common across terms - - // Physical distributions have the same DetectorModel+CrossSection context so we do not need to compare them - // We just need to check that these distributions have not been cancelled out for any terms - std::vector common_physical_dist_idxs; - for(unsigned int phys_idx=0; phys_idx(physical_distributions[phys_idx].get())) { - user_supplied_position_distribution = true; - continue; - } - // Remove distribution from distinct distributions - for(unsigned int i=0; i common_generation_dist_idxs; - - unsigned int i=0; - std::vector>> & gen_dists_0 = generation_distribution_state[i]; - for(unsigned int gen_idx_0=0; gen_idx_0> & gen_dist_0 = gen_dists_0[gen_idx_0]; - if(not gen_dist_0.first) - continue; - bool is_common = true; - std::vector common_idxs(injectors.size(), 0); - common_idxs[i] = gen_idx_0; - for(unsigned int j=i+1; j>> & gen_dists_1 = generation_distribution_state[j]; - for(unsigned int gen_idx_1=0; gen_idx_1> & gen_dist_1 = gen_dists_1[gen_idx_1]; - if(not gen_dist_1.first) - continue; - bool equivalent_dists = - gen_dist_0.second->AreEquivalent( // gen dist 0 - injectors[i]->GetDetectorModel(), // gen dist 0 context - injectors[i]->GetInteractions(), // gen dist 0 context - (std::shared_ptr)(gen_dist_1.second), // gen dist 1 - injectors[j]->GetDetectorModel(), // gen dist 1 context - injectors[j]->GetInteractions()); // gen dist 1 context - if(not equivalent_dists) - continue; - found_common = true; - common_idxs[j] = gen_idx_1; - break; // We found a gen dist cancelled out so we can skip additional comparisons - } - if(not found_common) { - // No matching distribution in this injector - // Term is not common across injectors - is_common = false; - // We can stop checking other injectors for this term - break; - } - } - if(not is_common) - continue; - // Remove distribution from distinct distribution list - for(unsigned int inj_idx=0; inj_idx dist = generation_distribution_state[0][gen_idx].second; - std::shared_ptr dist_detector = injectors[0]->GetDetectorModel(); - std::shared_ptr dist_interactions = injectors[0]->GetInteractions(); - std::function, std::shared_ptr, std::shared_ptr>)> predicate = [&] (std::tuple, std::shared_ptr, std::shared_ptr> p) -> bool { - return std::get<0>(p)->AreEquivalent(std::get<1>(p), std::get<2>(p), dist, dist_detector, dist_interactions); - }; - auto it = std::find_if(unique_distributions.begin(), unique_distributions.end(), predicate); - if(it != unique_distributions.end()) { - unsigned int index = std::distance(unique_distributions.begin(), it); - common_gen_idxs.push_back(index); - } else { - unique_distributions.push_back(std::make_tuple(dist, injectors[0]->GetDetectorModel(), injectors[0]->GetInteractions())); - common_gen_idxs.push_back(unique_distributions.size()-1); - } - } - - for(unsigned int phys_idx : common_physical_dist_idxs) { - std::shared_ptr dist = physical_distributions[phys_idx]; - std::function, std::shared_ptr, std::shared_ptr>)> predicate = [&] (std::tuple, std::shared_ptr, std::shared_ptr> p) -> bool { - return std::get<0>(p)->AreEquivalent(std::get<1>(p), std::get<2>(p), dist, detector_model, interactions); - }; - auto it = std::find_if(unique_distributions.begin(), unique_distributions.end(), predicate); - if(it != unique_distributions.end()) { - unsigned int index = std::distance(unique_distributions.begin(), it); - common_phys_idxs.push_back(index); - } else { - unique_distributions.push_back(std::make_tuple(dist, detector_model, interactions)); - common_phys_idxs.push_back(unique_distributions.size()-1); - } - } - - for(unsigned int injector_idx=0; injector_idx>> & phys_dists = physical_distribution_state[injector_idx]; - std::vector>> & gen_dists = generation_distribution_state[injector_idx]; - - std::vector gen_idxs; - std::vector phys_idxs; - for(unsigned int gen_idx=0; gen_idx dist = gen_dists[gen_idx].second; - // These are common to all injectors, so we pull information from the first injector - std::shared_ptr dist_detector = injectors[injector_idx]->GetDetectorModel(); - std::shared_ptr dist_interactions = injectors[injector_idx]->GetInteractions(); - std::function, std::shared_ptr, std::shared_ptr>)> predicate = [&] (std::tuple, std::shared_ptr, std::shared_ptr> p) -> bool { - return std::get<0>(p)->AreEquivalent(std::get<1>(p), std::get<2>(p), dist, dist_detector, dist_interactions); - }; - auto it = std::find_if(unique_distributions.begin(), unique_distributions.end(), predicate); - if(it != unique_distributions.end()) { - unsigned int index = std::distance(unique_distributions.begin(), it); - gen_idxs.push_back(index); - } else { - unique_distributions.push_back(std::make_tuple(dist, injectors[injector_idx]->GetDetectorModel(), injectors[injector_idx]->GetInteractions())); - gen_idxs.push_back(unique_distributions.size()-1); - } - } - - for(unsigned int phys_idx=0; phys_idx dist = phys_dists[phys_idx].second; - std::function, std::shared_ptr, std::shared_ptr>)> predicate = [&] (std::tuple, std::shared_ptr, std::shared_ptr> p) -> bool { - return std::get<0>(p)->AreEquivalent(std::get<1>(p), std::get<2>(p), dist, detector_model, interactions); - }; - auto it = std::find_if(unique_distributions.begin(), unique_distributions.end(), predicate); - if(it != unique_distributions.end()) { - unsigned int index = std::distance(unique_distributions.begin(), it); - phys_idxs.push_back(index); - } else { - unique_distributions.push_back(std::make_tuple(dist, detector_model, interactions)); - phys_idxs.push_back(unique_distributions.size()-1); - } - } - distinct_gen_idxs_by_injector.push_back(gen_idxs); - distinct_physical_idxs_by_injector.push_back(phys_idxs); - } - - //TODO - // Find unique contexts - // std::vector, std::shared_ptr>> unique_contexts; - // std::vector context_idx_by_injector; -} - -LeptonWeighter::LeptonWeighter(std::vector> injectors, std::shared_ptr detector_model, std::shared_ptr interactions, std::vector> physical_distributions) - : injectors(injectors) - , detector_model(detector_model) - , interactions(interactions) - , physical_distributions(physical_distributions) -{ - Initialize(); -} - -double LeptonWeighter::EventWeight(siren::dataclasses::InteractionRecord const & record) const { +double Weighter::EventWeight(siren::dataclasses::InteractionTree const & tree) const { // The weight is given by - // w = (\sum_i p_gen^i / p_phys^i)^-1 - + // + // [sum_{injectors i} + // x prod_{tree datum d} + // x (prod_{generation dist j} p_gen^{idj}) + // / (prod_{physical dist j} p_phys^{idj}) ] ^-1 + // // The generation probabilities are different between each injector. // Most of the physical probabilities are common between all injectors. // The physical interaction probability and physical position distribution @@ -519,121 +104,77 @@ double LeptonWeighter::EventWeight(siren::dataclasses::InteractionRecord const & // Thus, the two will cancel out and we are left with only the unnormalized position probability // w = p_physCommon / (\sum_i p_gen^i / p_physPosNonNorm^i) - // The ratio between physical and generation probabilities that differ between injectors - std::vector gen_over_phys; - gen_over_phys.reserve(injectors.size()); + - // From each injector we need the generation probability and the unnormalized position probability (interaction probability * position probability) - for(auto injector : injectors) { - double generation_probability = injector->GenerationProbability(record); - std::tuple bounds = injector->PrimaryInjectionBounds(record); + double inv_weight = 0; + for(unsigned int idx = 0; idx < injectors.size(); ++idx) { double physical_probability = 1.0; - - /* - if(user_supplied_position_distribution) { - // Need pos_prob * int_prob - // pos_prob already supplied - // just need int_prob - physical_probability *= InteractionProbability((std::shared_ptr)injector, record); - } else { - // Need pos_prob * int_prob - // nothing is already supplied - // need pos_prob and int_prob - // pos_prob * int_prob == unnormalized pos_prob - physical_probability *= UnnormalizedPositionProbability((std::shared_ptr)injector, record); + double generation_probability = injectors[idx]->EventsToInject();//GenerationProbability(tree); + for(auto const & datum : tree.tree) { + std::tuple bounds; + if(datum->depth() == 0) { + bounds = injectors[idx]->PrimaryInjectionBounds(datum->record); + physical_probability *= primary_process_weighters[idx]->PhysicalProbability(bounds, datum->record); + generation_probability *= primary_process_weighters[idx]->GenerationProbability(*datum); + } + else { + try { + bounds = injectors[idx]->SecondaryInjectionBounds(datum->record); + double phys_prob = secondary_process_weighter_maps[idx].at(datum->record.signature.primary_type)->PhysicalProbability(bounds, datum->record); + double gen_prob = secondary_process_weighter_maps[idx].at(datum->record.signature.primary_type)->GenerationProbability(*datum); + physical_probability *= phys_prob; + generation_probability *= gen_prob; + } catch(const std::out_of_range& oor) { + std::cout << "Out of Range error: " << oor.what() << '\n'; + return 0; + } + } } - */ - double prob = InteractionProbability(bounds, record); - physical_probability *= prob; - prob = NormalizedPositionProbability(bounds, record); - physical_probability *= prob; - prob = siren::injection::CrossSectionProbability(injector->GetDetectorModel(), injector->GetInteractions(), record); - physical_probability *= prob; - // Number of events is already in GenerationProbability - // double num_events = injector->EventsToInject(); - gen_over_phys.push_back(generation_probability / physical_probability); - } - - // The denominator is the sum over the ratios for each injector - double injection_specific_factors = accumulate(gen_over_phys.begin(), gen_over_phys.end()); - - // One physical probability density is computed for each distribution, independent of the injectors - double common_physical_probability = 1.0; - for(auto physical_distribution : physical_distributions) { - double prob = physical_distribution->GenerationProbability(detector_model, interactions, record); - common_physical_probability *= prob; + inv_weight += generation_probability / physical_probability; } + return 1./inv_weight; +} - double weight = common_physical_probability / injection_specific_factors; - return normalization * weight; +void Weighter::SaveWeighter(std::string const & filename) const { + std::ofstream os(filename+".siren_weighter", std::ios::binary); + ::cereal::BinaryOutputArchive archive(os); + this->save(archive,0); } -double LeptonWeighter::SimplifiedEventWeight(siren::dataclasses::InteractionRecord const & record) const { - std::vector probs; - probs.reserve(unique_distributions.size()); - for(unsigned int i=0; i, - std::shared_ptr, - std::shared_ptr - > const & p = unique_distributions[i]; - probs.push_back(std::get<0>(p)->GenerationProbability(std::get<1>(p), std::get<2>(p), record)); - } +void Weighter::LoadWeighter(std::string const & filename) { + std::cout << "Weighter loading not yet supported... sorry!\n"; + exit(0); + std::ifstream is(filename+".siren_weighter", std::ios::binary); + ::cereal::BinaryInputArchive archive(is); + //this->load(archive,0); +} - double phys_over_gen = 1.0; - for(unsigned int i=0; i> injectors, std::shared_ptr detector_model, std::shared_ptr primary_physical_process, std::vector> secondary_physical_processes) + : injectors(injectors) + , detector_model(detector_model) + , primary_physical_process(primary_physical_process) + , secondary_physical_processes(secondary_physical_processes) +{ + Initialize(); +} - std::vector gen_over_phys; - gen_over_phys.reserve(injectors.size()); - for(unsigned int i=0; iEventsToInject(); - for(unsigned int j=0; jGetDetectorModel(), injectors[i]->GetInteractions(), record); - prob *= cross_section_probability; - for(unsigned int j=0; j)injectors[i], record); - prob /= int_prob; - } else { - // Need pos_prob * int_prob - // nothing is already supplied - // need pos_prob and int_prob - // pos_prob * int_prob == unnormalized pos_prob - double pos_prob = UnnormalizedPositionProbability((std::shared_ptr)injectors[i], record); - prob /= pos_prob; - }*/ - std::tuple bounds = injectors[i]->PrimaryInjectionBounds(record); - double interaction_probability = InteractionProbability(bounds, record); - double normalized_position_probability = NormalizedPositionProbability(bounds, record); - prob /= interaction_probability; - prob /= normalized_position_probability; +Weighter::Weighter(std::vector> injectors, std::shared_ptr detector_model, std::shared_ptr primary_physical_process) + : injectors(injectors) + , detector_model(detector_model) + , primary_physical_process(primary_physical_process) + , secondary_physical_processes(std::vector>()) +{ + Initialize(); +} - // TODO - // Use unique contexts to compute cross section probability - gen_over_phys.push_back(prob); +Weighter::Weighter(std::vector> _injectors, std::string filename) { + LoadWeighter(filename); + if(_injectors.size() > 0) { + // overwrite the serialized injectors if the user have provided any + injectors = _injectors; } - - double gen_over_phys_d = accumulate(gen_over_phys.begin(), gen_over_phys.end()); - double weight = phys_over_gen / gen_over_phys_d; - return normalization * weight; + Initialize(); } } // namespace injection } // namespace siren - diff --git a/projects/injection/private/pybindings/injection.cxx b/projects/injection/private/pybindings/injection.cxx index c2512b6b..044ccd06 100644 --- a/projects/injection/private/pybindings/injection.cxx +++ b/projects/injection/private/pybindings/injection.cxx @@ -3,11 +3,6 @@ #include "../../public/SIREN/injection/Process.h" #include "../../public/SIREN/injection/Injector.h" -//#include "../../public/SIREN/injection/ColumnDepthSIREN.h" -//#include "../../public/SIREN/injection/CylinderVolumeSIREN.h" -//#include "../../public/SIREN/injection/DecayRangeSIREN.h" -//#include "../../public/SIREN/injection/RangedSIREN.h" -#include "../../public/SIREN/injection/TreeWeighter.h" #include "../../public/SIREN/injection/Weighter.h" #include "../../public/SIREN/injection/WeightingUtils.h" @@ -126,17 +121,12 @@ PYBIND11_MODULE(injection,m) { .def("GenerationProbability",&SecondaryProcessWeighter::GenerationProbability) .def("EventWeight",&SecondaryProcessWeighter::EventWeight); - class_>(m, "TreeWeighter") + class_>(m, "Weighter") .def(init>, std::shared_ptr, std::shared_ptr, std::vector>>()) .def(init>, std::shared_ptr, std::shared_ptr>()) .def(init>, std::string>()) - .def("EventWeight",&TreeWeighter::EventWeight) - .def("SaveWeighter",&TreeWeighter::SaveWeighter) - .def("LoadWeighter",&TreeWeighter::LoadWeighter) + .def("EventWeight",&Weighter::EventWeight) + .def("SaveWeighter",&Weighter::SaveWeighter) + .def("LoadWeighter",&Weighter::LoadWeighter) ; - - class_>(m, "LeptonWeighter") - .def(init>, std::shared_ptr, std::shared_ptr, std::vector>>()) - .def("EventWeight",&LeptonWeighter::EventWeight) - .def("SimplifiedEventWeight",&LeptonWeighter::SimplifiedEventWeight); } diff --git a/projects/injection/private/test/CCM_HNL_TEST.cxx b/projects/injection/private/test/CCM_HNL_TEST.cxx index 83b1b830..6704ee50 100644 --- a/projects/injection/private/test/CCM_HNL_TEST.cxx +++ b/projects/injection/private/test/CCM_HNL_TEST.cxx @@ -17,7 +17,7 @@ #include "SIREN/dataclasses/Particle.h" #include "SIREN/injection/Injector.h" #include "SIREN/injection/Process.h" -#include "SIREN/injection/TreeWeighter.h" +#include "SIREN/injection/Weighter.h" #include "SIREN/geometry/Geometry.h" #include "SIREN/geometry/ExtrPoly.h" #include "SIREN/geometry/Sphere.h" @@ -303,8 +303,8 @@ TEST(Injector, Generation) upper_injector->SetStoppingCondition(stopping_condition); lower_injector->SetStoppingCondition(stopping_condition); - std::shared_ptr upper_weighter = std::make_shared(std::vector>{upper_injector}, detector_model, primary_physical_process_upper_injector, secondary_physical_processes); - std::shared_ptr lower_weighter = std::make_shared(std::vector>{lower_injector}, detector_model, primary_physical_process_lower_injector, secondary_physical_processes); + std::shared_ptr upper_weighter = std::make_shared(std::vector>{upper_injector}, detector_model, primary_physical_process_upper_injector, secondary_physical_processes); + std::shared_ptr lower_weighter = std::make_shared(std::vector>{lower_injector}, detector_model, primary_physical_process_lower_injector, secondary_physical_processes); int i = 0; diff --git a/projects/injection/public/SIREN/injection/TreeWeighter.h b/projects/injection/public/SIREN/injection/TreeWeighter.h deleted file mode 100644 index ecd12cb3..00000000 --- a/projects/injection/public/SIREN/injection/TreeWeighter.h +++ /dev/null @@ -1,127 +0,0 @@ -#pragma once -#ifndef SIREN_TreeWeighter_H -#define SIREN_TreeWeighter_H - -#include // for map -#include -#include // for shared_ptr -#include // for vector - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "SIREN/dataclasses/InteractionTree.h" // for InteractionT... -#include "SIREN/dataclasses/Particle.h" // for Particle - -namespace siren { namespace dataclasses { class InteractionRecord; } } -namespace siren { namespace detector { class DetectorModel; } } -namespace siren { namespace distributions { class PrimaryInjectionDistribution; } } -namespace siren { namespace distributions { class SecondaryInjectionDistribution; } } -namespace siren { namespace distributions { class WeightableDistribution; } } -namespace siren { namespace injection { class Injector; } } -namespace siren { namespace injection { class PrimaryInjectionProcess; } } -namespace siren { namespace injection { class SecondaryInjectionProcess; } } -namespace siren { namespace injection { class PhysicalProcess; } } -namespace siren { namespace math { class Vector3D; } } - -namespace siren { -namespace injection { - -// Class handling weight calculation for a single pair of injection and physical processes -template -class ProcessWeighter { -private: - std::shared_ptr phys_process; - std::shared_ptr inj_process; - std::vector> unique_gen_distributions; - std::vector> unique_phys_distributions; - std::shared_ptr detector_model; - std::vector> const & GetInjectionDistributions(); - void Initialize(); - double normalization; -public: - double InteractionProbability(std::tuple const & bounds, siren::dataclasses::InteractionRecord const & record) const; - double NormalizedPositionProbability(std::tuple const & bounds, siren::dataclasses::InteractionRecord const & record) const; - double PhysicalProbability(std::tuple const & bounds, siren::dataclasses::InteractionRecord const & record) const; - double GenerationProbability(siren::dataclasses::InteractionTreeDatum const & datum) const; - double EventWeight(std::tuple const & bounds, siren::dataclasses::InteractionTreeDatum const & datum) const; - ProcessWeighter(std::shared_ptr phys_process, std::shared_ptr inj_process, std::shared_ptr detector_model); - -}; // ProcessWeighter - -typedef ProcessWeighter PrimaryProcessWeighter; -typedef ProcessWeighter SecondaryProcessWeighter; - -// Parent class for calculating event weights -// Assumes there is a unique secondary physical process for each particle type -class TreeWeighter { -private: - // Supplied by constructor - std::vector> injectors; - std::shared_ptr detector_model; - std::shared_ptr primary_physical_process; - std::vector> secondary_physical_processes; - - // Calculated upon initialization - std::vector> primary_process_weighters; - std::vector< - std::map< - siren::dataclasses::ParticleType, - std::shared_ptr - > - > secondary_process_weighter_maps; - - void Initialize(); -public: - double EventWeight(siren::dataclasses::InteractionTree const & tree) const; - TreeWeighter(std::vector> injectors, std::shared_ptr detector_model, std::shared_ptr primary_physical_process, std::vector> secondary_physical_processes); - TreeWeighter(std::vector> injectors, std::shared_ptr detector_model, std::shared_ptr primary_physical_process); - TreeWeighter(std::vector> injectors, std::string filename); - void SaveWeighter(std::string const & filename) const; - void LoadWeighter(std::string const & filename); - - template - void save(Archive & archive, std::uint32_t const version) const { - if(version == 0) { - archive(::cereal::make_nvp("Injectors", injectors)); - archive(::cereal::make_nvp("DetectorModel", detector_model)); - archive(::cereal::make_nvp("PrimaryPhysicalProcess", primary_physical_process)); - archive(::cereal::make_nvp("SecondaryPhysicalProcesses", secondary_physical_processes)); - } else { - throw std::runtime_error("TreeWeighter only supports version <= 0!"); - } - } - - template - void load(Archive & archive, std::uint32_t const version) const { - if(version == 0) { - archive(::cereal::make_nvp("Injectors", injectors)); - archive(::cereal::make_nvp("DetectorModel", detector_model)); - archive(::cereal::make_nvp("PrimaryPhysicalProcess", primary_physical_process)); - archive(::cereal::make_nvp("SecondaryPhysicalProcesses", secondary_physical_processes)); - } else { - throw std::runtime_error("TreeWeighter only supports version <= 0!"); - } - } - -}; // TreeWeighter - - -} //namespace injection -} //namespace siren - -#include "TreeWeighter.tcc" - -CEREAL_CLASS_VERSION(siren::injection::TreeWeighter, 0); - - -#endif // SIREN_TreeWeighter_H diff --git a/projects/injection/public/SIREN/injection/Weighter.h b/projects/injection/public/SIREN/injection/Weighter.h index d0658af3..93e845cb 100644 --- a/projects/injection/public/SIREN/injection/Weighter.h +++ b/projects/injection/public/SIREN/injection/Weighter.h @@ -2,69 +2,126 @@ #ifndef SIREN_Weighter_H #define SIREN_Weighter_H -#include // for tuple -#include // for shared_ptr -#include // for vector +#include // for map +#include +#include // for shared_ptr +#include // for vector #include +#include #include #include #include +#include +#include +#include #include #include #include -namespace siren { namespace interactions { class InteractionCollection; } } +#include "SIREN/dataclasses/InteractionTree.h" // for InteractionT... +#include "SIREN/dataclasses/Particle.h" // for Particle + namespace siren { namespace dataclasses { class InteractionRecord; } } namespace siren { namespace detector { class DetectorModel; } } +namespace siren { namespace distributions { class PrimaryInjectionDistribution; } } +namespace siren { namespace distributions { class SecondaryInjectionDistribution; } } namespace siren { namespace distributions { class WeightableDistribution; } } namespace siren { namespace injection { class Injector; } } +namespace siren { namespace injection { class PrimaryInjectionProcess; } } +namespace siren { namespace injection { class SecondaryInjectionProcess; } } +namespace siren { namespace injection { class PhysicalProcess; } } namespace siren { namespace math { class Vector3D; } } namespace siren { namespace injection { -class LeptonWeighter { +// Class handling weight calculation for a single pair of injection and physical processes +template +class ProcessWeighter { private: - std::vector> injectors; + std::shared_ptr phys_process; + std::shared_ptr inj_process; + std::vector> unique_gen_distributions; + std::vector> unique_phys_distributions; std::shared_ptr detector_model; - //TODO Think about whether we want to pass a CrossSection collection, or a vector of cross sections - //TODO Think about what to do with multiple neutrino primary types. Do we want to support mutiple types across one InteractionCollection, across one Injector, across one LeptonWeighter? - std::shared_ptr interactions; - std::vector> physical_distributions; - - std::vector, std::shared_ptr, std::shared_ptr>> unique_distributions; - std::vector common_gen_idxs; - std::vector common_phys_idxs; - std::vector> distinct_gen_idxs_by_injector; - std::vector> distinct_physical_idxs_by_injector; - std::vector, std::shared_ptr>> unique_contexts; - std::vector context_idx_by_injector; + std::vector> const & GetInjectionDistributions(); + void Initialize(); double normalization; +public: + double InteractionProbability(std::tuple const & bounds, siren::dataclasses::InteractionRecord const & record) const; + double NormalizedPositionProbability(std::tuple const & bounds, siren::dataclasses::InteractionRecord const & record) const; + double PhysicalProbability(std::tuple const & bounds, siren::dataclasses::InteractionRecord const & record) const; + double GenerationProbability(siren::dataclasses::InteractionTreeDatum const & datum) const; + double EventWeight(std::tuple const & bounds, siren::dataclasses::InteractionTreeDatum const & datum) const; + ProcessWeighter(std::shared_ptr phys_process, std::shared_ptr inj_process, std::shared_ptr detector_model); + +}; // ProcessWeighter + +typedef ProcessWeighter PrimaryProcessWeighter; +typedef ProcessWeighter SecondaryProcessWeighter; + +// Parent class for calculating event weights +// Assumes there is a unique secondary physical process for each particle type +class Weighter { +private: + // Supplied by constructor + std::vector> injectors; + std::shared_ptr detector_model; + std::shared_ptr primary_physical_process; + std::vector> secondary_physical_processes; - bool user_supplied_position_distribution = false; + // Calculated upon initialization + std::vector> primary_process_weighters; + std::vector< + std::map< + siren::dataclasses::ParticleType, + std::shared_ptr + > + > secondary_process_weighter_maps; void Initialize(); public: - //TODO Think about the relationship between interaction probability and the positional distribution. Check that the math works out - //TODO Add versions of these functions that take precomputed intersections - double InteractionProbability(std::shared_ptr injector, siren::dataclasses::InteractionRecord const & record) const; - double InteractionProbability(std::tuple bounds, siren::dataclasses::InteractionRecord const & record) const; - double UnnormalizedPositionProbability(std::shared_ptr injector, siren::dataclasses::InteractionRecord const & record) const; - double UnnormalizedPositionProbability(std::tuple bounds, siren::dataclasses::InteractionRecord const & record) const; - double NormalizedPositionProbability(std::tuple bounds, siren::dataclasses::InteractionRecord const & record) const; - //TODO Add a function to check that we have the right match up of variables between generator and physical distribution - //TODO Figure out a way to check that physical and generation probabilities match, and ignore those when weighting - LeptonWeighter(std::vector> injectors, std::shared_ptr detector_model, std::shared_ptr interactions, std::vector> physical_distributions); - double EventWeight(siren::dataclasses::InteractionRecord const & record) const; - double SimplifiedEventWeight(siren::dataclasses::InteractionRecord const & record) const; -}; + double EventWeight(siren::dataclasses::InteractionTree const & tree) const; + Weighter(std::vector> injectors, std::shared_ptr detector_model, std::shared_ptr primary_physical_process, std::vector> secondary_physical_processes); + Weighter(std::vector> injectors, std::shared_ptr detector_model, std::shared_ptr primary_physical_process); + Weighter(std::vector> injectors, std::string filename); + void SaveWeighter(std::string const & filename) const; + void LoadWeighter(std::string const & filename); + + template + void save(Archive & archive, std::uint32_t const version) const { + if(version == 0) { + archive(::cereal::make_nvp("Injectors", injectors)); + archive(::cereal::make_nvp("DetectorModel", detector_model)); + archive(::cereal::make_nvp("PrimaryPhysicalProcess", primary_physical_process)); + archive(::cereal::make_nvp("SecondaryPhysicalProcesses", secondary_physical_processes)); + } else { + throw std::runtime_error("Weighter only supports version <= 0!"); + } + } + + template + void load(Archive & archive, std::uint32_t const version) const { + if(version == 0) { + archive(::cereal::make_nvp("Injectors", injectors)); + archive(::cereal::make_nvp("DetectorModel", detector_model)); + archive(::cereal::make_nvp("PrimaryPhysicalProcess", primary_physical_process)); + archive(::cereal::make_nvp("SecondaryPhysicalProcesses", secondary_physical_processes)); + } else { + throw std::runtime_error("Weighter only supports version <= 0!"); + } + } + +}; // Weighter + } //namespace injection } //namespace siren -CEREAL_CLASS_VERSION(siren::injection::LeptonWeighter, 0); +#include "Weighter.tcc" +CEREAL_CLASS_VERSION(siren::injection::Weighter, 0); -#endif // SIREN_Weighter_H +#endif // SIREN_Weighter_H diff --git a/projects/injection/public/SIREN/injection/TreeWeighter.tcc b/projects/injection/public/SIREN/injection/Weighter.tcc similarity index 98% rename from projects/injection/public/SIREN/injection/TreeWeighter.tcc rename to projects/injection/public/SIREN/injection/Weighter.tcc index ae6fa48d..fa85d0f5 100644 --- a/projects/injection/public/SIREN/injection/TreeWeighter.tcc +++ b/projects/injection/public/SIREN/injection/Weighter.tcc @@ -1,7 +1,7 @@ #pragma once -#ifndef SIREN_TreeWeighter_TCC -#define SIREN_TreeWeighter_TCC -#include "SIREN/injection/TreeWeighter.h" +#ifndef SIREN_Weighter_TCC +#define SIREN_Weighter_TCC +#include "SIREN/injection/Weighter.h" #include // for ite... #include // for array @@ -262,4 +262,4 @@ std::vector