Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Dev/d meson #74

Open
wants to merge 13 commits into
base: main
Choose a base branch
from
Open

Dev/d meson #74

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
8b5eb26
finished implementing D meson DB, but lower energies decay still is a…
Jun 12, 2024
f1fb7cf
add examples
Jun 12, 2024
3ce2b8f
preparing for PR
MiaochenJin Aug 22, 2024
a6d66aa
Update DetectorModel.cxx for PR
nickkamp1 Aug 22, 2024
4414cbb
Update Injector.cxx to fixes failed_tries bug
nickkamp1 Aug 22, 2024
e1da327
Merge branch 'Harvard-Neutrino:main' into dev/DMeson
MiaochenJin Aug 22, 2024
7e21f18
updated file for PR comments, preparing for PR again
MiaochenJin Aug 22, 2024
9f5b7d9
Delete resources/Examples/DMesonExample/paper.mplstyle
MiaochenJin Aug 22, 2024
a9c3ac1
Delete resources/Examples/DMesonExample/parse_output.py
MiaochenJin Aug 22, 2024
6a54bd0
Merge branch 'Harvard-Neutrino:main' into dev/DMeson
MiaochenJin Oct 9, 2024
8ea969d
update handling of numerical instabilities
MiaochenJin Nov 13, 2024
abf6883
new stuff
MiaochenJin Nov 19, 2024
4dfb9ce
quarkDIS completed
MiaochenJin Dec 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions projects/dataclasses/private/Particle.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,5 +83,16 @@ bool isCharged(ParticleType p){
p==ParticleType::Hadrons);
}

bool isQuark(Particle::ParticleType p){
return(p==Particle::ParticleType::Charm || p==Particle::ParticleType::CharmBar);
}

bool isHadron(Particle::ParticleType p){
return(p==Particle::ParticleType::Hadrons ||
p==Particle::ParticleType::D0 || p==Particle::ParticleType::D0Bar ||
p==Particle::ParticleType::DPlus || p==Particle::ParticleType::DMinus);
}


} // namespace utilities
} // namespace siren
3 changes: 3 additions & 0 deletions projects/dataclasses/public/SIREN/dataclasses/Particle.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,9 @@ class Particle {
bool isLepton(ParticleType p);
bool isCharged(ParticleType p);
bool isNeutrino(ParticleType p);
bool isQuark(Particle::ParticleType p);
bool isHadron(Particle::ParticleType p);


} // namespace dataclasses
} // namespace siren
Expand Down
6 changes: 6 additions & 0 deletions projects/dataclasses/public/SIREN/dataclasses/ParticleTypes.def
View file
Open in desktop
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@nickkamp1 be sure to update this when you merge your HNL branch

Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,12 @@ X(TauPlus, -15)
X(TauMinus, 15)
X(NuTau, 16)
X(NuTauBar, -16)
X(Charm, 4)
X(CharmBar, -4)
X(K0, 311)
X(K0Bar, -311)



/* Nuclei */
X(HNucleus, 1000010010)
Expand Down
6 changes: 5 additions & 1 deletion projects/detector/private/DetectorModel.cxx
View file
Open in desktop
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Was there a specific issue that the distance.magnitude() <= 1e-5 checks fixed?

I'm wondering if there is another underlying issue that this is covering up, or if we really need this fix.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For the direction.magnitude() <= 1e-5 checks, you should put the 1e-5 constant into the projects/detector/public/SIREN/detector/DetectorModel.h header

class DetectorModel {
public:
    constexpr static double distance_threshold = 1e-5;

and then use this constant to perform the checks in the .cxx 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) {
if(direction.magnitude() == 0 || direction.magnitude() <= 1e-6) {
direction = intersections.direction;
} else {
direction.normalize();
Expand Down Expand Up @@ -978,6 +978,7 @@ double DetectorModel::GetInteractionDepthInCGS(Geometry::IntersectionList const
std::vector<siren::dataclasses::ParticleType> const & targets,
std::vector<double> const & total_cross_sections,
double const & total_decay_length) const {

if(p0 == p1) {
return 0.0;
}
Expand All @@ -986,6 +987,9 @@ double DetectorModel::GetInteractionDepthInCGS(Geometry::IntersectionList const
if(distance == 0.0) {
return 0.0;
}
if(direction.magnitude() <= 1e-6) {
direction = intersections.direction;
}
direction.normalize();

double dot = intersections.direction * direction;
Expand Down
17 changes: 17 additions & 0 deletions projects/distributions/private/secondary/vertex/SecondaryBoundedVertexDistribution.cxx
View file
Open in desktop
MiaochenJin marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Instead of adding special login in this class checks for a hadronization interaction to change the class behavior, you should implement a new secondary vertex distribution that has the appropriate behavior for hadronization.

Original file line number Diff line number Diff line change
Expand Up @@ -53,10 +53,17 @@ 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;

siren::math::Vector3D endcap_0 = pos;
// skip computation of endpoint if interaction is hadronization
if (interactions->HasHadronizations()) {
record.SetLength(0);
return;
}
siren::math::Vector3D endcap_1 = endcap_0 + max_length * dir;

siren::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), max_length);
Expand Down Expand Up @@ -116,11 +123,21 @@ 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);

siren::math::Vector3D endcap_0 = record.primary_initial_position;
// hadrnoization treated differently, but also check for misconducting
if (interactions->HasHadronizations()) {
if (vertex == endcap_0) {
return 1.0;
} else {
return 0.0;
}
}
siren::math::Vector3D endcap_1 = endcap_0 + max_length * dir;

siren::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), max_length);
Expand Down
21 changes: 21 additions & 0 deletions projects/distributions/private/secondary/vertex/SecondaryPhysicalVertexDistribution.cxx
View file
Open in desktop
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same comment as for the bounded vertex distribution:

Instead of adding special login in this class checks for a hadronization interaction to change the class behavior, you should implement a new secondary vertex distribution that has the appropriate behavior for hadronization.

Original file line number Diff line number Diff line change
Expand Up @@ -51,10 +51,21 @@ 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 @@ -75,6 +86,7 @@ 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 @@ -96,11 +108,20 @@ 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);

siren::math::Vector3D endcap_0 = record.primary_initial_position;
if (interactions->HasHadronizations()) {
if (vertex == endcap_0) {
return 1.0;
} else {
return 0.0;
}
}

siren::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), std::numeric_limits<double>::infinity());
path.ClipToOuterBounds();
Expand Down
1 change: 1 addition & 0 deletions projects/distributions/private/secondary/vertex/SecondaryVertexPositionDistribution.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ 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
Loading