Skip to content

Commit

Permalink
updated file for PR comments, preparing for PR again
Browse files Browse the repository at this point in the history
  • Loading branch information
@MiaochenJin
MiaochenJin committed Aug 22, 2024
1 parent e1da327 commit 7e21f18
Show file tree
Hide file tree
Showing 13 changed files with 69 additions and 590 deletions.
1 change: 1 addition & 0 deletions projects/dataclasses/public/SIREN/dataclasses/ParticleTypes.def
View file
Original file line number Diff line number Diff line change
Expand Up @@ -194,6 +194,7 @@ X(NuclInt, -2000001004)
X(MuPair, -2000001005)
X(Hadrons, -2000001006)
X(Decay, -2000001007)
X(Hadronization, 99)
X(ContinuousEnergyLoss, -2000001111)
X(FiberLaser, -2000002100)
X(N2Laser, -2000002101)
Expand Down
5 changes: 3 additions & 2 deletions projects/detector/private/DetectorModel.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -683,7 +683,7 @@ double DetectorModel::GetInteractionDensity(Geometry::IntersectionList const & i
std::vector<double> const & total_cross_sections,
double const & total_decay_length) const {
Vector3D direction = p0 - intersections.position;
if(direction.magnitude() == 0 || direction.magnitude() <= 1e-6) {
if(direction.magnitude() == 0 || direction.magnitude() <= 1e-5) {
direction = intersections.direction;
} else {
direction.normalize();
Expand Down Expand Up @@ -987,8 +987,9 @@ double DetectorModel::GetInteractionDepthInCGS(Geometry::IntersectionList const
if(distance == 0.0) {
return 0.0;
}
if(direction.magnitude() <= 1e-6) {
if(direction.magnitude() <= 1e-5) {
direction = intersections.direction;
// return 1e-05; // I have to do this to ensure that it works
}
direction.normalize();

Expand Down
4 changes: 0 additions & 4 deletions projects/distributions/private/secondary/vertex/SecondaryBoundedVertexDistribution.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,8 +53,6 @@ double log_one_minus_exp_of_negative(double x) {


void SecondaryBoundedVertexDistribution::SampleVertex(std::shared_ptr<siren::utilities::SIREN_random> rand, std::shared_ptr<siren::detector::DetectorModel const> detector_model, std::shared_ptr<siren::interactions::InteractionCollection const> interactions, siren::dataclasses::SecondaryDistributionRecord & record) const {
// std::cout << "in sample bounded vertex" << std::endl;

siren::math::Vector3D pos = record.initial_position;
siren::math::Vector3D dir = record.direction;

Expand Down Expand Up @@ -123,8 +121,6 @@ void SecondaryBoundedVertexDistribution::SampleVertex(std::shared_ptr<siren::uti
}

double SecondaryBoundedVertexDistribution::GenerationProbability(std::shared_ptr<siren::detector::DetectorModel const> detector_model, std::shared_ptr<siren::interactions::InteractionCollection const> interactions, siren::dataclasses::InteractionRecord const & record) const {
// std::cout << "in sample bounded vertex gen prob" << std::endl;

siren::math::Vector3D dir(record.primary_momentum[1], record.primary_momentum[2], record.primary_momentum[3]);
dir.normalize();
siren::math::Vector3D vertex(record.interaction_vertex);
Expand Down
10 changes: 0 additions & 10 deletions projects/distributions/private/secondary/vertex/SecondaryPhysicalVertexDistribution.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,22 +51,15 @@ double log_one_minus_exp_of_negative(double x) {


void SecondaryPhysicalVertexDistribution::SampleVertex(std::shared_ptr<siren::utilities::SIREN_random> rand, std::shared_ptr<siren::detector::DetectorModel const> detector_model, std::shared_ptr<siren::interactions::InteractionCollection const> interactions, siren::dataclasses::SecondaryDistributionRecord & record) const {
// // std::cout << "in sample physical vertex" << std::endl;
siren::math::Vector3D pos = record.initial_position;
siren::math::Vector3D dir = record.direction;
// // std::cout << "in sample physical vertex-1" << std::endl;

siren::math::Vector3D endcap_0 = pos;
// treat hadronizations differntely
if (interactions->HasHadronizations()) {
// std::cout << "in sample physical vertex-hadron" << std::endl;

record.SetLength(0);
return;
}
// // std::cout << "in sample physical vertex-shouldnm't be here" << std::endl;


siren::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), std::numeric_limits<double>::infinity());
path.ClipToOuterBounds();

Expand All @@ -86,7 +79,6 @@ void SecondaryPhysicalVertexDistribution::SampleVertex(std::shared_ptr<siren::ut

double total_interaction_depth = path.GetInteractionDepthInBounds(targets, total_cross_sections, total_decay_length);
if(total_interaction_depth == 0) {
// std::cout << "error is here" << std::endl;
throw(siren::utilities::InjectionFailure("No available interactions along path!"));
}

Expand All @@ -108,8 +100,6 @@ void SecondaryPhysicalVertexDistribution::SampleVertex(std::shared_ptr<siren::ut
}

double SecondaryPhysicalVertexDistribution::GenerationProbability(std::shared_ptr<siren::detector::DetectorModel const> detector_model, std::shared_ptr<siren::interactions::InteractionCollection const> interactions, siren::dataclasses::InteractionRecord const & record) const {
// std::cout << "in sample physical vertex gen prob" << std::endl;

siren::math::Vector3D dir(record.primary_momentum[1], record.primary_momentum[2], record.primary_momentum[3]);
dir.normalize();
siren::math::Vector3D vertex(record.interaction_vertex);
Expand Down
1 change: 0 additions & 1 deletion projects/distributions/private/secondary/vertex/SecondaryVertexPositionDistribution.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@ namespace distributions {
// class SecondaryVertexPositionDistribution : InjectionDistribution
//---------------
void SecondaryVertexPositionDistribution::Sample(std::shared_ptr<siren::utilities::SIREN_random> rand, std::shared_ptr<siren::detector::DetectorModel const> detector_model, std::shared_ptr<siren::interactions::InteractionCollection const> interactions, siren::dataclasses::SecondaryDistributionRecord & record) const {
// // // // // std::cout << "sampling vertex" << std::endl;
SampleVertex(rand, detector_model, interactions, record);
}

Expand Down
108 changes: 4 additions & 104 deletions projects/injection/private/Injector.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,7 +188,6 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record
double fake_prob;
// if contains hadronization, then perform only hadronization
if (interactions->HasHadronizations()) {
// std::cout << "saw hadronization" << std::endl;
double total_frag_prob = 0;
std::vector<double> frag_probs;
for(auto const & hadronization : interactions->GetHadronizations() ) {
Expand All @@ -209,9 +208,6 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record
}
}

// std::cout << "Hadronization finished signatures" << std::endl;


// now choose the specific charmed hadron to fragment into
double r = random->Uniform(0, total_frag_prob);
unsigned int index = 0;
Expand All @@ -224,34 +220,22 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record

matching_hadronizations[index]->SampleFinalState(xsec_record, random);
xsec_record.Finalize(record);

// std::cout << "hadronization done!" << std::endl;

} else {
if (interactions->HasCrossSections()) {
// std::cout << "saw xsec" << std::endl;
for(auto const target : available_targets) {
if(possible_targets.find(target) != possible_targets.end()) {
// std::cout << "saw xsec: in first for loop" << std::endl;
// Get target density
double target_density = detector_model->GetParticleDensity(intersections, DetectorPosition(interaction_vertex), target);
// Loop over cross sections that have this target
std::vector<std::shared_ptr<siren::interactions::CrossSection>> const & target_cross_sections = interactions->GetCrossSectionsForTarget(target);
for(auto const & cross_section : target_cross_sections) {
// std::cout << "saw xsec: in second for loop" << std::endl;
// 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) {
// std::cout << "saw xsec: in third for loop" << std::endl;

fake_record.signature = signature;
fake_record.target_mass = detector_model->GetTargetMass(target);
// Add total cross section times density to the total prob
// std::cout << "about to sample total xsec" << std::endl;

fake_prob = target_density * cross_section->TotalCrossSection(fake_record);
// std::cout << "finished sampling total xsec" << std::endl;

total_prob += fake_prob;
xsec_prob += fake_prob;
// Add total prob to probs
Expand All @@ -266,7 +250,6 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record
}
}
if (interactions->HasDecays()) {
// std::cout << "saw decay" << std::endl;
for(auto const & decay : interactions->GetDecays() ) {
for(auto const & signature : decay->GetPossibleSignaturesFromParent(record.signature.primary_type)) {
fake_record.signature = signature;
Expand All @@ -282,7 +265,6 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record
}
}
}
// std::cout << "continuing...." << std::endl;
if(total_prob == 0)
throw(siren::utilities::InjectionFailure("No valid interactions for this event!"));
// Throw a random number
Expand All @@ -293,7 +275,6 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record
record.signature.target_type = matching_targets[index];
record.signature = matching_signatures[index];
double selected_prob = 0.0;
// std::cout << "finished initializing stuff" << std::endl;
for(unsigned int i=0; i<probs.size(); ++i) {
if(matching_signatures[index] == matching_signatures[i]) {
selected_prob += (i > 0 ? probs[i] - probs[i - 1] : probs[i]);
Expand All @@ -303,14 +284,9 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record
throw(siren::utilities::InjectionFailure("No valid interactions for this event!"));
record.target_mass = detector_model->GetTargetMass(record.signature.target_type);
siren::dataclasses::CrossSectionDistributionRecord xsec_record(record);
// std::cout << "finished sampling from primary process" << std::endl;
if(r <= xsec_prob) {
// std::cout << "about to sample primary process final state" << std::endl;

matching_cross_sections[index]->SampleFinalState(xsec_record, random);
} else {
// std::cout << "about to sample primary process final state" << std::endl;

matching_decays[index - matching_cross_sections.size()]->SampleFinalState(xsec_record, random);
}
xsec_record.Finalize(record);
Expand All @@ -321,8 +297,6 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record
//
// Modifies an InteractionRecord with the new event
siren::dataclasses::InteractionRecord Injector::SampleSecondaryProcess(siren::dataclasses::SecondaryDistributionRecord & secondary_record) const {
// std::cout << "sampling secondary" << std::endl;
// std::cout << "secondary record type is " << secondary_record.type << " " << secondary_record.id << std::endl;
std::shared_ptr<siren::injection::SecondaryInjectionProcess> secondary_process = secondary_process_map.at(secondary_record.type);
std::shared_ptr<siren::interactions::InteractionCollection> secondary_interactions = secondary_process->GetInteractions();
std::vector<std::shared_ptr<siren::distributions::SecondaryInjectionDistribution>> secondary_distributions = secondary_process->GetSecondaryInjectionDistributions();
Expand All @@ -332,32 +306,15 @@ siren::dataclasses::InteractionRecord Injector::SampleSecondaryProcess(siren::da
size_t failed_tries = 0;
while(true) {
tries += 1;
// // std::cout << "gotcha" << std::endl;
// for(auto & distribution : secondary_distributions) {
// // std::cout << "sample distribution" << std::endl;
// distribution->Sample(random, detector_model, secondary_process->GetInteractions(), secondary_record);
// }
// siren::dataclasses::InteractionRecord record;
// secondary_record.Finalize(record);
// // // std::cout << "sample distribution" << std::endl;

// SampleCrossSection(record, secondary_interactions);
// return record;

try {
for(auto & distribution : secondary_distributions) {
// std::cout << "sample distribution" << std::endl;
distribution->Sample(random, detector_model, secondary_process->GetInteractions(), secondary_record);
}
siren::dataclasses::InteractionRecord record;
secondary_record.Finalize(record);
// // std::cout << "sample distribution" << std::endl;

SampleCrossSection(record, secondary_interactions);
return record;
} catch(siren::utilities::InjectionFailure const & e) {
// std::cout << "caught error" << std::endl;

failed_tries += 1;
if(failed_tries > max_tries) {
throw(siren::utilities::InjectionFailure("Failed to generate secondary process!"));
Expand All @@ -379,7 +336,6 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() {
size_t tries = 0;
size_t failed_tries = 0;
// Initial Process
// std::cout << "Sampling primary interactions" << std::endl;
while(true) {
tries += 1;
try {
Expand All @@ -388,7 +344,6 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() {
distribution->Sample(random, detector_model, primary_process->GetInteractions(), primary_record);
}
primary_record.Finalize(record);
// std::cout << "primary record fixed" << std:: endl;
SampleCrossSection(record);
break;
} catch(siren::utilities::InjectionFailure const & e) {
Expand All @@ -408,12 +363,9 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() {
std::shared_ptr<siren::dataclasses::InteractionTreeDatum> parent = tree.add_entry(record);

// Secondary Processes
// std::cout << "Sampling primary interactions 2" << std::endl;

std::deque<std::tuple<std::shared_ptr<siren::dataclasses::InteractionTreeDatum>, std::shared_ptr<siren::dataclasses::SecondaryDistributionRecord>>> secondaries;
std::function<void(std::shared_ptr<siren::dataclasses::InteractionTreeDatum>)> add_secondaries = [&](std::shared_ptr<siren::dataclasses::InteractionTreeDatum> parent) {
for(size_t i=0; i<parent->record.signature.secondary_types.size(); ++i) {
// // std::cout << "for loop 1" << std::endl;
siren::dataclasses::ParticleType const & type = parent->record.signature.secondary_types[i];
std::map<siren::dataclasses::ParticleType, std::shared_ptr<siren::injection::SecondaryInjectionProcess>>::iterator it = secondary_process_map.find(type);
if(it == secondary_process_map.end()) {
Expand All @@ -430,70 +382,18 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() {
};

add_secondaries(parent);
// std::cout << "num secondaries: " << secondaries.size() << std::endl;
while(secondaries.size() > 0) {
// // std::cout << "while loop 1" << std::endl;

for(int i = secondaries.size() - 1; i >= 0; --i) {
// // std::cout << "for loop 2" << std::endl;

std::shared_ptr<siren::dataclasses::InteractionTreeDatum> parent = std::get<0>(secondaries[i]);
std::shared_ptr<siren::dataclasses::SecondaryDistributionRecord> secondary_dist = std::get<1>(secondaries[i]);
// // std::cout << "for loop 2-1" << std::endl;

secondaries.erase(secondaries.begin() + i);

// std::cout << "Primary Type: " << secondary_dist->record.signature.primary_type << std::endl;
// std::cout << "Secondary Types: ";
// for (const auto& type : secondary_dist->record.signature.secondary_types) {
// std::cout << type << " ";
// }
// std::cout << std::endl;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// this try-catch clock is to debug the secondary process LE no available interaction error
try{
siren::dataclasses::InteractionRecord secondary_record = SampleSecondaryProcess(*secondary_dist);
std::shared_ptr<siren::dataclasses::InteractionTreeDatum> secondary_datum = tree.add_entry(secondary_record, parent);
add_secondaries(secondary_datum);
} catch (const std::exception& e) {
std::cerr << "Error occurred: " << e.what() << std::endl;

// Print the primary type and secondary types for debugging
std::cerr << "Primary Type: " << secondary_dist->record.signature.primary_type << std::endl;
std::cerr << "Secondary Types: ";
for (const auto& type : secondary_dist->record.signature.secondary_types) {
std::cerr << type << " ";
}
std::cerr << std::endl;

// Print the primary momentum
std::cerr << "Primary Momentum: ";
for (double component : secondary_dist->record.primary_momentum) {
std::cerr << component << " ";
}
std::cerr << std::endl;

// Print the secondary IDs
std::cerr << "Secondary IDs: ";
for (const auto& id : secondary_dist->record.secondary_ids) {
std::cerr << id << " ";
}
std::cerr << std::endl;
throw;
} catch (...) {
std::cerr << "Unknown exception caught!" << std::endl;
throw;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


siren::dataclasses::InteractionRecord secondary_record = SampleSecondaryProcess(*secondary_dist);
std::shared_ptr<siren::dataclasses::InteractionTreeDatum> secondary_datum = tree.add_entry(secondary_record, parent);
add_secondaries(secondary_datum);


}
// // std::cout << "while loop 1-2" << std::endl;

}
injected_events += 1;
return tree;
Expand Down
2 changes: 1 addition & 1 deletion projects/injection/private/Weighter.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,7 @@ double Weighter::EventWeight(siren::dataclasses::InteractionTree const & tree) c
double generation_probability = injectors[idx]->EventsToInject();//GenerationProbability(tree);
for(auto const & datum : tree.tree) {
// skip weighting if record contains hadronization
if (isQuark(datum->record.signature.primary_type)) {
if (datum->record.signature.target_type == siren::dataclasses::Particle::ParticleType::Hadronization) {
continue;
}
std::tuple<siren::math::Vector3D, siren::math::Vector3D> bounds;
Expand Down
Loading

0 comments on commit 7e21f18

Please sign in to comment.