Fork/alexwenym

CMakeLists.txt

@@ -68,6 +68,10 @@ include(delabella)
 include(photospline)
 include(googletest)
 
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -I${CFITSIO_INCLUDE_DIR}")                                                                                                         
+set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -I${CFITSIO_INCLUDE_DIR}")                                                                                             
+set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -I${CFITSIO_INCLUDE_DIR}")
+
 # load macros for googletest
 include(testing)
 

projects/crosssections/private/DISFromSpline.cxx

@@ -54,36 +54,54 @@ bool kinematicallyAllowed(double x, double y, double E, double M, double m) {
 
 DISFromSpline::DISFromSpline() {}
 
-DISFromSpline::DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minimum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types) : primary_types_(primary_types), target_types_(target_types), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
+DISFromSpline::DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minimum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types, std::string units) : primary_types_(primary_types), target_types_(target_types), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
     LoadFromMemory(differential_data, total_data);
     InitializeSignatures();
+    SetUnits(units);
 }
 
-DISFromSpline::DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minimum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
+DISFromSpline::DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minimum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types, std::string units) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
     LoadFromMemory(differential_data, total_data);
     InitializeSignatures();
+    SetUnits(units);
 }
 
-DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minimum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types) : primary_types_(primary_types), target_types_(target_types), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
+DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minimum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types, std::string units) : primary_types_(primary_types), target_types_(target_types), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
     LoadFromFile(differential_filename, total_filename);
     InitializeSignatures();
+    SetUnits(units);
 }
 
-DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types) : primary_types_(primary_types), target_types_(target_types) {
+DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types, std::string units) : primary_types_(primary_types), target_types_(target_types) {
     LoadFromFile(differential_filename, total_filename);
     ReadParamsFromSplineTable();
     InitializeSignatures();
+    SetUnits(units);
 }
 
-DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minimum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
+DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minimum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types, std::string units) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
     LoadFromFile(differential_filename, total_filename);
     InitializeSignatures();
+    SetUnits(units);
 }
 
-DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()) {
+DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types, std::string units) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()) {
     LoadFromFile(differential_filename, total_filename);
     ReadParamsFromSplineTable();
     InitializeSignatures();
+    SetUnits(units);
+}
+
+void DISFromSpline::SetUnits(std::string units) {
+    std::transform(units.begin(), units.end(), units.begin(),
+        [](unsigned char c){ return std::tolower(c); });
+    if(units == "cm") {
+        unit = 1.0;
+    } else if(units == "m") {
+        unit = 10000.0;
+    } else {
+        throw std::runtime_error("Cross section units not supported!");
+    }
 }
 
 bool DISFromSpline::equal(CrossSection const & other) const {
@@ -263,7 +281,7 @@ double DISFromSpline::TotalCrossSection(LI::dataclasses::Particle::ParticleType
     total_cross_section_.searchcenters(&log_energy, &center);
     double log_xs = total_cross_section_.ndsplineeval(&log_energy, &center, 0);
 
-    return std::pow(10.0, log_xs);
+    return unit * std::pow(10.0, log_xs);
 }
 
 // No implementation for DIS yet, just use non-target function
@@ -336,7 +354,7 @@ double DISFromSpline::DifferentialCrossSection(double energy, double x, double y
         return 0;
     double result = pow(10., differential_cross_section_.ndsplineeval(coordinates.data(), centers.data(), 0));
     assert(result >= 0);
-    return result;
+    return unit * result;
 }
 
 double DISFromSpline::InteractionThreshold(dataclasses::InteractionRecord const & interaction) const {

projects/crosssections/private/pybindings/crosssections.cxx

@@ -6,7 +6,6 @@
 #include "../../public/LeptonInjector/crosssections/CrossSectionCollection.h"
 #include "../../public/LeptonInjector/crosssections/DISFromSpline.h"
 #include "../../public/LeptonInjector/crosssections/HNLFromSpline.h"
-#include "../../public/LeptonInjector/crosssections/CrossSectionCollection.h"
 #include "../../public/LeptonInjector/crosssections/DipoleFromTable.h"
 #include "../../public/LeptonInjector/crosssections/DarkNewsCrossSection.h"
 #include "../../public/LeptonInjector/crosssections/DarkNewsDecay.h"

projects/crosssections/public/LeptonInjector/crosssections/DISFromSpline.h

@@ -49,15 +49,19 @@ friend cereal::access;
     int interaction_type_;
     double target_mass_;
     double minimum_Q2_;
+
+    double unit;
 
 public:
     DISFromSpline();
-    DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minumum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types);
-    DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minumum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types);
-    DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minumum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types);
-    DISFromSpline(std::string differential_filename, std::string total_filename, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types);
-    DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minumum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types);
-    DISFromSpline(std::string differential_filename, std::string total_filename, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types);
+    DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minumum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types, std::string units = "cm");
+    DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minumum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types, std::string units = "cm");
+    DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minumum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types, std::string units = "cm");
+    DISFromSpline(std::string differential_filename, std::string total_filename, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types, std::string units = "cm");
+    DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minumum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types, std::string units = "cm");
+    DISFromSpline(std::string differential_filename, std::string total_filename, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types, std::string units = "cm");
+
+    void SetUnits(std::string units);
 
     virtual bool equal(CrossSection const & other) const override;
 

projects/detector/private/EarthModel.cxx

@@ -578,6 +578,9 @@ double EarthModel::GetInteractionDensity(Geometry::IntersectionList const & inte
             interaction_density = 0.0;
             for(unsigned int i=0; i<targets.size(); ++i) {
                 interaction_density += density * particle_fractions[i] * total_cross_sections[i];
+                //std::cout << "        EarthModel::GetInteractionDensity: density " << density << std::endl;
+                //std::cout << "        EarthModel::GetInteractionDensity: particle_fractions[i] " << particle_fractions[i] << std::endl;
+                //std::cout << "        EarthModel::GetInteractionDensity: total_cross_sections[i] " << total_cross_sections[i] << std::endl;
             }
             interaction_density *= 100; // cm^-1 --> m^-1
             return true;
@@ -697,6 +700,10 @@ double EarthModel::DistanceForColumnDepthFromPoint(Geometry::IntersectionList co
             EarthSector sector = GetSector(current_intersection->hierarchy);
             double target = column_depth - total_column_depth;
             double distance = sector.density->InverseIntegral(p0+start_point*direction, direction, target, segment_length);
+            //std::cout << "    EarthModel::DistanceForColumnDepthFromPoint: p0+start_point*direction (" << (p0+start_point*direction).GetX() << ", " << (p0+start_point*direction).GetY() << ", " << (p0+start_point*direction).GetZ() << ")" << std::endl;
+            //std::cout << "    EarthModel::DistanceForColumnDepthFromPoint: direction (" << direction.GetX() << ", " << direction.GetY() << ", " << direction.GetZ() << ")" << std::endl;
+            //std::cout << "    EarthModel::DistanceForColumnDepthFromPoint: (target, segment_length) (" << target << ", " << segment_length << ")" <<std::endl;
+
             done = distance >= 0;
             double integral = sector.density->Integral(p0+start_point*direction, direction, segment_length);
             total_column_depth += integral;
@@ -873,9 +880,14 @@ double EarthModel::GetInteractionDepthInCGS(Geometry::IntersectionList const & i
             double segment_length = end_point - start_point;
             EarthSector sector = GetSector(current_intersection->hierarchy);
             double integral = sector.density->Integral(p0+start_point*direction, direction, segment_length);
+
+            //std::cout << "        EarthModel::GetInteractionDepthInCGS: segment_length: " << segment_length << std::endl; 
+            //std::cout << "        EarthModel::GetInteractionDepthInCGS: integral: " << integral << std::endl; 
             std::vector<double> particle_fractions = materials_.GetTargetParticleFraction(sector.material_id, targets.begin(), targets.end());
             for(unsigned int i=0; i<targets.size(); ++i) {
                 interaction_depths[i] += (integral * 100) * particle_fractions[i]; // cm^-3 * m --> cm^-2
+                //std::cout << "        EarthModel::GetInteractionDepthInCGS: particle_fractions[i]: " << particle_fractions[i] << std::endl;
+                //std::cout << "        EarthModel::GetInteractionDepthInCGS: interaction_depths[i]: " << interaction_depths[i] << std::endl; 
             }
         }
         // last_point = end_point;