Skip to content

Commit

Permalink
update handling of numerical instabilities
Browse files Browse the repository at this point in the history
  • Loading branch information
@MiaochenJin
MiaochenJin committed Nov 13, 2024
1 parent 6a54bd0 commit 8ea969d
Show file tree
Hide file tree
Showing 15 changed files with 1,107 additions and 21 deletions.
1 change: 1 addition & 0 deletions projects/dataclasses/private/pybindings/dataclasses.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,7 @@ PYBIND11_MODULE(dataclasses,m) {
.def(init<>())
.def("__str__", [](InteractionRecord const & r) { std::stringstream ss; ss << r; return ss.str(); })
.def_readwrite("signature",&InteractionRecord::signature)
.def_readwrite("primary_initial_position",&InteractionRecord::primary_initial_position)
.def_readwrite("primary_mass",&InteractionRecord::primary_mass)
.def_readwrite("primary_momentum",&InteractionRecord::primary_momentum)
.def_readwrite("primary_helicity",&InteractionRecord::primary_helicity)
Expand Down
5 changes: 5 additions & 0 deletions projects/distributions/private/secondary/vertex/SecondaryBoundedVertexDistribution.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,6 +188,11 @@ double SecondaryBoundedVertexDistribution::GenerationProbability(std::shared_ptr
prob_density = interaction_density * exp(-log_one_minus_exp_of_negative(total_interaction_depth) - traversed_interaction_depth);
}

if (prob_density == 0) {
std::cout << "observed prob density 0 in physical vertex under process " << record.signature.primary_type << " to "
<< record.signature.secondary_types[0] << " with total depth " << total_interaction_depth << std::endl;
}

return prob_density;
}

Expand Down
5 changes: 5 additions & 0 deletions projects/distributions/private/secondary/vertex/SecondaryPhysicalVertexDistribution.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,6 +148,11 @@ double SecondaryPhysicalVertexDistribution::GenerationProbability(std::shared_pt
prob_density = interaction_density * exp(-log_one_minus_exp_of_negative(total_interaction_depth) - traversed_interaction_depth);
}

if (prob_density == 0) {
std::cout << "observed prob density 0 in physical vertex under process " << record.signature.primary_type << " to "
<< record.signature.secondary_types[0] << " with total depth " << total_interaction_depth << std::endl;
}

return prob_density;
}

Expand Down
4 changes: 4 additions & 0 deletions projects/injection/private/Injector.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,6 +159,8 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record
throw(siren::utilities::InjectionFailure("No particle interaction!"));
}

//std::cout << "in sample cross section" << std::endl;

std::set<siren::dataclasses::ParticleType> const & possible_targets = interactions->TargetTypes();

siren::math::Vector3D interaction_vertex(
Expand Down Expand Up @@ -285,6 +287,7 @@ 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;
matching_cross_sections[index]->SampleFinalState(xsec_record, random);
} else {
matching_decays[index - matching_cross_sections.size()]->SampleFinalState(xsec_record, random);
Expand Down Expand Up @@ -339,6 +342,7 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() {
while(true) {
tries += 1;
try {
//std::cout << "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);
Expand Down
36 changes: 36 additions & 0 deletions projects/injection/private/Weighter.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,9 +110,11 @@ double Weighter::EventWeight(siren::dataclasses::InteractionTree const & tree) c

double inv_weight = 0;
for(unsigned int idx = 0; idx < injectors.size(); ++idx) {
// std::cout << "New Event" << std::endl;
double physical_probability = 1.0;
double generation_probability = injectors[idx]->EventsToInject();//GenerationProbability(tree);
for(auto const & datum : tree.tree) {
// std::cout << "new depth " << datum->depth() << std::endl;
// skip weighting if record contains hadronization
if (datum->record.signature.target_type == siren::dataclasses::Particle::ParticleType::Hadronization) {
continue;
Expand All @@ -122,6 +124,17 @@ double Weighter::EventWeight(siren::dataclasses::InteractionTree const & tree) c
bounds = injectors[idx]->PrimaryInjectionBounds(datum->record);
physical_probability *= primary_process_weighters[idx]->PhysicalProbability(bounds, datum->record);
generation_probability *= primary_process_weighters[idx]->GenerationProbability(*datum);
// for debugging purposes: nan weights are frequnetly detected
if (physical_probability == 0) {
std::cout << "zero physics depth 0: " << datum->record.signature.primary_type << std::endl;
} else if (std::isinf(physical_probability)) {
std::cout << "inf physics depth 0: " << datum->record.signature.primary_type << std::endl;
}
if (generation_probability == 0) {
std::cout << "zero gen depth 0: " << datum->record.signature.primary_type << std::endl;
} else if (std::isinf(generation_probability)) {
std::cout << "inf gen depth 0: " << datum->record.signature.primary_type << std::endl;
}
}
else {
try {
Expand All @@ -130,13 +143,36 @@ double Weighter::EventWeight(siren::dataclasses::InteractionTree const & tree) c
double gen_prob = secondary_process_weighter_maps[idx].at(datum->record.signature.primary_type)->GenerationProbability(*datum);
physical_probability *= phys_prob;
generation_probability *= gen_prob;
// if (phys_prob == 0) {
// std::cout << "zero physics: " << datum->record.signature.primary_type << std::endl;
// } else if (std::isinf(phys_prob)) {
// std::cout << "inf physics: " << datum->record.signature.primary_type << std::endl;
// }
// if (gen_prob == 0) {
// std::cout << "zero gen: " << datum->record.signature.primary_type << std::endl;
// } else if (std::isinf(gen_prob)) {
// std::cout << "inf gen: " << datum->record.signature.primary_type << std::endl;
// }
} catch(const std::out_of_range& oor) {
std::cout << "Out of Range error: " << oor.what() << '\n';
return 0;
}
}
}
inv_weight += generation_probability / physical_probability;

// if (physical_probability == 0) {
// std::cout << "Event has 0 physical probability, leading to: " << inv_weight << " " << 1./inv_weight << std::endl;
// } else if (physical_probability != physical_probability) {
// std::cout << "Event has inf physical probability, leading to: " << inv_weight << " " << 1./inv_weight << std::endl;
// }
// if (generation_probability == 0) {
// std::cout << "Event has 0 generation probability, leading to: " << inv_weight << " " << 1./inv_weight << std::endl;
// } else if (generation_probability != generation_probability) {
// std::cout << "Event has inf generation probability, leading to: " << inv_weight << " " << 1./inv_weight << std::endl;
// }
// std::cout << "gen and physics prob is " << generation_probability << " " << physical_probability << std::endl;
// std::cout << "inverse weight and final weight is " << inv_weight << " " << 1./inv_weight << std::endl;
}
return 1./inv_weight;
}
Expand Down
19 changes: 17 additions & 2 deletions projects/injection/private/WeightingUtils.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,15 +69,30 @@ double CrossSectionProbability(std::shared_ptr<siren::detector::DetectorModel co
// Loop over cross section signatures with the same target
std::vector<siren::dataclasses::InteractionSignature> signatures = cross_section->GetPossibleSignaturesFromParents(record.signature.primary_type, target);
for(auto const & signature : signatures) {
// check here for 0 generation probability
fake_record.signature = signature;
fake_record.target_mass = detector_model->GetTargetMass(target);
// Add total cross section times density to the total prob
double target_prob = target_density * cross_section->TotalCrossSection(fake_record);
double total_xs = cross_section->TotalCrossSection(fake_record);
double target_prob = target_density * total_xs;
// if (total_xs == 0) {
// std::cout << "total cross section give 0 for process of " << record.signature.primary_type << std::endl;
// std::cout << "for signature " << fake_record.signature << std::endl;
// } else if (std::isinf(total_xs)) {
// std::cout << "total cross section give inf for process of " << record.signature.primary_type << std::endl;
// std::cout << "for signature " << fake_record.signature << std::endl;
// }
total_prob += target_prob;
// Add up total cross section times density times final state prob for matching signatures
if(signature == record.signature) {
// selected_prob += target_prob;
selected_final_state += target_prob * cross_section->FinalStateProbability(record);
double final_prob = cross_section->FinalStateProbability(record);
// if (final_prob == 0) {
// std::cout << "final state prob give 0 for process of " << record.signature.primary_type << std::endl;
// } else if (std::isinf(final_prob)) {
// std::cout << "final state prob give inf for process of " << record.signature.primary_type << std::endl;
// }
selected_final_state += target_prob * final_prob;
}
}
}
Expand Down
23 changes: 22 additions & 1 deletion projects/injection/public/SIREN/injection/Weighter.tcc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -204,16 +204,28 @@ double ProcessWeighter<ProcessType>::PhysicalProbability(std::tuple<siren::math:

double physical_probability = 1.0;
double prob = InteractionProbability(bounds, record);
// if (prob == 0) {
// std::cout << "Interaction probability is 0" << std::endl;
// }
physical_probability *= prob;

prob = NormalizedPositionProbability(bounds, record);
// if (prob == 0) {
// std::cout << "Position probability is 0" << std::endl;
// }
physical_probability *= prob;

prob = siren::injection::CrossSectionProbability(detector_model, phys_process->GetInteractions(), record);
// if (prob == 0) {
// std::cout << "XSec probability is 0" << std::endl;
// }
physical_probability *= prob;

for(auto physical_dist : unique_phys_distributions) {
physical_probability *= physical_dist->GenerationProbability(detector_model, phys_process->GetInteractions(), record);
// if (physical_dist->GenerationProbability(detector_model, phys_process->GetInteractions(), record) == 0) {
// std::cout << "physical dist Generation probablity is 0" << std::endl;
// }
}

return normalization * physical_probability;
Expand All @@ -222,10 +234,19 @@ double ProcessWeighter<ProcessType>::PhysicalProbability(std::tuple<siren::math:
template<typename ProcessType>
double ProcessWeighter<ProcessType>::GenerationProbability(siren::dataclasses::InteractionTreeDatum const & datum ) const {
double gen_probability = siren::injection::CrossSectionProbability(detector_model, inj_process->GetInteractions(), datum.record);

// if (gen_probability == 0) {
// std::cout << "Gen Cross section probability is 0" << std::endl;
// }
for(auto gen_dist : unique_gen_distributions) {
// if (gen_dist->GenerationProbability(detector_model, inj_process->GetInteractions(), datum.record) == 0) {
// std::cout << "generation dist gen probability is 0" << std::endl;
// }
gen_probability *= gen_dist->GenerationProbability(detector_model, inj_process->GetInteractions(), datum.record);
}

// if (gen_probability == 0) {
// std::cout << "tcc file gen prob is 0" << std::endl;
// }
return gen_probability;
}

Expand Down
Loading

0 comments on commit 8ea969d

Please sign in to comment.