Charm DIS injection

projects/dataclasses/private/Particle.cxx

@@ -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

projects/dataclasses/public/SIREN/dataclasses/Particle.h

@@ -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

projects/dataclasses/public/SIREN/dataclasses/ParticleTypes.def

@@ -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)

projects/detector/private/DetectorModel.cxx

@@ -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();
@@ -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;
     }
@@ -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;

projects/distributions/private/secondary/vertex/SecondaryBoundedVertexDistribution.cxx

@@ -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);
@@ -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);

projects/distributions/private/secondary/vertex/SecondaryPhysicalVertexDistribution.cxx

@@ -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();
@@ -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!"));
     }
 
@@ -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();

projects/distributions/private/secondary/vertex/SecondaryVertexPositionDistribution.cxx

@@ -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);
 }