new stuff

Harvard-Neutrino · Nov 19, 2024 · abf6883 · abf6883
1 parent 8ea969d
commit abf6883
Show file tree

Hide file tree

Showing 11 changed files with 458 additions and 82 deletions.
diff --git a/projects/dataclasses/private/InteractionRecord.cxx b/projects/dataclasses/private/InteractionRecord.cxx
@@ -327,13 +327,22 @@ void PrimaryDistributionRecord::FinalizeAvailable(InteractionRecord & record) co
 }
 
 void PrimaryDistributionRecord::Finalize(InteractionRecord & record) const {
+    //std::cout << "starting finalize" << std::endl;
+    //std::cout << record.signature << std::endl;
+    //std::cout << record.primary_momentum[0] << std::endl;
+
     record.signature.primary_type = type;
     record.primary_id = GetID();
     record.interaction_vertex = GetInteractionVertex();
     record.primary_initial_position = GetInitialPosition();
+    //std::cout << "in primary record finalize" << std::endl;
     record.primary_mass = GetMass();
     record.primary_momentum = GetFourMomentum();
     record.primary_helicity = GetHelicity();
+    //std::cout << "finished primary record finalize" << std::endl;
+
+    //std::cout << record.signature << std::endl;
+    //std::cout << record.primary_momentum[0] << std::endl;
 }
 
 /////////////////////////////////////////
@@ -553,11 +562,16 @@ void SecondaryParticleRecord::UpdateMomentum() const {
 }
 
 void SecondaryParticleRecord::Finalize(InteractionRecord & record) const {
+    ////std::cout << "SecPartRecord::Finalize : starting for " << type << std::endl;
+
     assert(record.signature.secondary_types.at(secondary_index) == type);
     record.secondary_ids.at(secondary_index) = GetID();
     record.secondary_masses.at(secondary_index) = GetMass();
     record.secondary_momenta.at(secondary_index) = GetFourMomentum();
     record.secondary_helicities.at(secondary_index) = GetHelicity();
+    ////std::cout << "SecPartRecord::Finalize : finished for " << type << " with" << std::endl;
+    ////std::cout << record << std::endl;
+
 }
 
 /////////////////////////////////////////
@@ -676,6 +690,11 @@ SecondaryParticleRecord const & CrossSectionDistributionRecord::GetSecondaryPart
 }
 
 void CrossSectionDistributionRecord::Finalize(InteractionRecord & record) const {
+    //std::cout << "XsecDistRecord::Finalize : starting" << std::endl;
+    //std::cout << record.signature << std::endl;
+    //std::cout << signature << std::endl;
+    //std::cout << primary_momentum[0] << std::endl;
+
     record.target_id = target_id;
     record.target_mass = target_mass;
     record.target_helicity = target_helicity;
@@ -688,6 +707,7 @@ void CrossSectionDistributionRecord::Finalize(InteractionRecord & record) const
     record.secondary_helicities.resize(secondary_particles.size());
 
     for(SecondaryParticleRecord const & secondary : secondary_particles) {
+        //std::cout << "XsecDistRecord::Finalize : going to secondaries: " << secondary << std::endl;
         secondary.Finalize(record);
     }
 }

diff --git a/projects/dataclasses/private/Particle.cxx b/projects/dataclasses/private/Particle.cxx
@@ -93,6 +93,10 @@ bool isHadron(Particle::ParticleType p){
 				p==Particle::ParticleType::DPlus || p==Particle::ParticleType::DMinus);
 	}
 
-
+bool isD(Particle::ParticleType p){
+		return(p==Particle::ParticleType::D0 || p==Particle::ParticleType::D0Bar ||
+				p==Particle::ParticleType::DPlus || p==Particle::ParticleType::DMinus);
+	}
+
 } // namespace utilities
 } // namespace siren
diff --git a/projects/dataclasses/public/SIREN/dataclasses/Particle.h b/projects/dataclasses/public/SIREN/dataclasses/Particle.h
@@ -76,6 +76,7 @@ bool isCharged(ParticleType p);
 bool isNeutrino(ParticleType p);
 bool isQuark(Particle::ParticleType p);
 bool isHadron(Particle::ParticleType p);
+bool isD(Particle::ParticleType p);
 
 
 } // namespace dataclasses

diff --git a/projects/injection/private/Injector.cxx b/projects/injection/private/Injector.cxx
@@ -159,7 +159,7 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record
         throw(siren::utilities::InjectionFailure("No particle interaction!"));
     }
 
-    //std::cout << "in sample cross section" << std::endl;
+    ////std::cout << "in sample cross section" << std::endl;
 
     std::set<siren::dataclasses::ParticleType> const & possible_targets = interactions->TargetTypes();
 
@@ -267,6 +267,7 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record
                 }
             }
         }
+        //std::cout << "injector :: sample cross sections: after obtaining signatures" << std::endl;
         if(total_prob == 0)
             throw(siren::utilities::InjectionFailure("No valid interactions for this event!"));
         // Throw a random number
@@ -287,12 +288,17 @@ 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;
+            //std::cout << "injector::sample cross section: going into sampel final state" << std::endl;
             matching_cross_sections[index]->SampleFinalState(xsec_record, random);
+            //std::cout << "injector::sample cross section: finished sampling" << std::endl;
+
         } else {
             matching_decays[index - matching_cross_sections.size()]->SampleFinalState(xsec_record, random);
         }
+        ////std::cout << "injector::sample cross section: calling finalizing" << std::endl;
         xsec_record.Finalize(record);
+        ////std::cout << "injector::sample cross section: finished finalizing" << std::endl;
+
     }
 }
 
@@ -341,13 +347,15 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() {
     // Initial Process
     while(true) {
         tries += 1;
+        ////std::cout << "injector::GenerateEvent : trying to generate with trial number " << tries << std::endl;
         try {
-            //std::cout << "generating primary process" << std::endl;
+            ////std::cout << "injector::GenerateEvent : 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);
             }
             primary_record.Finalize(record);
+            //std::cout << "injector::GenerateEvent : sampling cross section" << std::endl;
             SampleCrossSection(record);
             break;
         } catch(siren::utilities::InjectionFailure const & e) {
@@ -367,6 +375,7 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() {
     std::shared_ptr<siren::dataclasses::InteractionTreeDatum> parent = tree.add_entry(record);
 
     // Secondary Processes
+    //std::cout << "injector::GenerateEvent : sampling secondary process" << 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) {
@@ -399,6 +408,7 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() {
 
         }
     }
+    //std::cout << "finished sampling secondary process" << std::endl;
     injected_events += 1;
     return tree;
 }

diff --git a/projects/interactions/CMakeLists.txt b/projects/interactions/CMakeLists.txt
@@ -21,6 +21,8 @@ LIST (APPEND interactions_SOURCES
     ${PROJECT_SOURCE_DIR}/projects/interactions/private/CharmHadronization.cxx
     ${PROJECT_SOURCE_DIR}/projects/interactions/private/CharmMesonDecay.cxx
     ${PROJECT_SOURCE_DIR}/projects/interactions/private/DMesonELoss.cxx
+    ${PROJECT_SOURCE_DIR}/projects/interactions/private/QuarkDISFromSpline.cxx
+
 )
 add_library(SIREN_interactions OBJECT ${interactions_SOURCES})
 set_property(TARGET SIREN_interactions PROPERTY POSITION_INDEPENDENT_CODE ON)

diff --git a/projects/interactions/private/CharmDISFromSpline.cxx b/projects/interactions/private/CharmDISFromSpline.cxx
@@ -45,6 +45,8 @@ namespace {
 bool kinematicallyAllowed(double x, double y, double E, double M, double m) {
     if(x > 1) //Eq. 6 right inequality
         return false;
+    // this is to get rid of the infinities as a temporary solution
+    if (x < 1e-6 or y < 1e-6) return false;
     if(x < ((m * m) / (2 * M * (E - m)))) //Eq. 6 left inequality
         return false;
     //denominator of a and b