diff --git a/configure b/configure index a165ef43..b08e6eea 100755 --- a/configure +++ b/configure @@ -709,6 +709,7 @@ OBJECTS = $(patsubst src/%.cpp,build/%.o,$(SOURCES)) EXAMPLES := examples/Single_energy/single_energy \ examples/Multiple_energy/multiple_energy \ examples/Atm_default/atm_default \ + examples/Emitting_Atmosphere/emitting_atmosphere \ examples/Bodies/bodies \ examples/Xsections/xsections \ examples/NSI/nsi \ @@ -799,6 +800,10 @@ examples/Atm_default/atm_default : $(DYN_PRODUCT) examples/Atm_default/main.cpp @echo Compiling atmospheric example @$(CXX) $(EXAMPLES_FLAGS) examples/Atm_default/main.cpp -lnuSQuIDS $(LDFLAGS) -o $@ +examples/Emitting_Atmosphere/emitting_atmosphere : $(DYN_PRODUCT) examples/Emitting_Atmosphere/main.cpp + @echo Compiling emitting atmosphere example + @$(CXX) $(EXAMPLES_FLAGS) examples/Emitting_Atmosphere/main.cpp -lnuSQuIDS $(LDFLAGS) -o $@ + build/exBody.o : examples/Bodies/exBody.h examples/Bodies/exBody.cpp @$(CXX) $(EXAMPLES_FLAGS) -c examples/Bodies/exBody.cpp -o $@ diff --git a/examples/Emitting_Atmosphere/main.cpp b/examples/Emitting_Atmosphere/main.cpp new file mode 100644 index 00000000..7a66ac05 --- /dev/null +++ b/examples/Emitting_Atmosphere/main.cpp @@ -0,0 +1,189 @@ + /****************************************************************************** + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + * Authors: * + * Carlos Arguelles (University of Wisconsin Madison) * + * carguelles@icecube.wisc.edu * + * Jordi Salvado (University of Wisconsin Madison) * + * jsalvado@icecube.wisc.edu * + * Christopher Weaver (University of Wisconsin Madison) * + * chris.weaver@icecube.wisc.edu * + ******************************************************************************/ + +#include +#include +#include +#include + +/* + * The problem of solving the propagation of the atmospheric neutrinos is an + * energy and zenith dependent problem, for this we include the class nuSQUIDSAtm + * that allows to solve a set of nuSUIDS energy dependent objects to take in to account the + * zenith dependence. The problem of flux insertion into the atmosphere due to cosmic ray sources + * during evolution requires a class derived from EarthAtm to be called into nuSQUIDSAtm + */ + + +using namespace nusquids; + + +//Function that sets the initial flux if you want to add flux in at the top of the atmosphere +//for this example we set it to 0 +double flux_function(double enu, double cz){ + return 0; +} + + +int main() +{ + //Units and constants class + squids::Const units; + //Number of neutrinos (3) standard 4 1-sterile + //Number of neutrinos, 4 is the case with one sterile neutrino + std::cout << "Enter number of neutrino flavors to consider:\n"; + std::cout << "(3) Three Active Neutrinos, " << "(4) 3+1 Three Active and One Sterile Neutrino" << std::endl; + unsigned int numneu; + std::cin >>numneu; + if( not(numneu==3 || numneu==4)){ + throw std::runtime_error("Only (3) or (4) are valid options"); + } + std::cout << "Consider Earth absorption and interactions? " << std::endl; + std::string use_int; + std::cin >> use_int; + bool interactions = false; + if(use_int=="yes" || use_int=="y") + interactions = true; + + //Minimum and maximum values for the energy and cosine zenith, notice that the energy needs to have the + //units, if they are omitted the input is in eV. + double Emin=1e-2*units.GeV; + double Emax=1.e6*units.GeV; + double czmin=-1; + double czmax=1; + + //Declaration of the atmospheric object + //Must include the nuSQUIDS template or nuSQUIDSNSI user defined class derived from nuSQUIDS in order to + //use the EmittingEarthAtm class or another user defined class derived from EarthAtm + std::cout << "Begin: constructing nuSQuIDS-Atm object" << std::endl; + nuSQUIDSAtm nus_atm(linspace(czmin,czmax,40),logspace(Emin,Emax,100),numneu,both,interactions); + std::cout << "End: constructing nuSQuIDS-Atm object" << std::endl; + + //Sets height of atmosphere in km (defaults to 22 km when not set) + nus_atm.Set_AtmHeight(60); + + //Functions Necessary for injection + //Passes the energy range into the emitting bodies for interpolation + nus_atm.Set_AtmEmissionEnergies(); + //Allows the nusquids objects to accept neutrino sources + nus_atm.Set_NeutrinoSources(true); + + std::cout << "Begin: setting mixing angles." << std::endl; + // set mixing angles, mass differences and cp phases + nus_atm.Set_MixingAngle(0,1,0.563942); + nus_atm.Set_MixingAngle(0,2,0.154085); + nus_atm.Set_MixingAngle(1,2,0.785398); + + nus_atm.Set_SquareMassDifference(1,7.65e-05); + nus_atm.Set_SquareMassDifference(2,0.00247); + + nus_atm.Set_CPPhase(0,2,0); + if(numneu > 3){ + nus_atm.Set_SquareMassDifference(3,-1.); + nus_atm.Set_MixingAngle(1,3,0.160875); + } + std::cout << "End: setting mixing angles." << std::endl; + + //Setup integration precision + nus_atm.Set_rel_error(1.0e-6); + nus_atm.Set_abs_error(1.0e-6); + nus_atm.Set_GSL_step(gsl_odeiv2_step_rk4); + + //Array that contains the values of the energies and cosine of the zenith, is the same length for every zenith + auto e_range = nus_atm.GetERange(); + auto cz_range = nus_atm.GetCosthRange(); + + std::cout << "Begin: setting initial state." << std::endl; + + //Construct the initial state at the top of the atmosphere, we set it to 0 here + marray inistate{nus_atm.GetNumCos(),nus_atm.GetNumE(),2,numneu}; + std::fill(inistate.begin(),inistate.end(),0); + for ( int ci = 0 ; ci < nus_atm.GetNumCos(); ci++){ + for ( int ei = 0 ; ei < nus_atm.GetNumE(); ei++){ + for ( int rho = 0; rho < 2; rho ++ ){ + for (int flv = 0; flv < numneu; flv++){ + inistate[ci][ei][rho][flv] = (flv == 0) ? flux_function(e_range[ei], cz_range[ci]) : 0;//set the flux for the muon flavor: 1 + inistate[ci][ei][rho][flv] = (flv == 1) ? flux_function(e_range[ei], cz_range[ci]) : 0;//set the flux for the muon flavor: 1 + } + } + } + } + + //Set the initial state in the atmSQuIDS object + nus_atm.Set_initial_state(inistate,flavor); + std::cout << "End: setting initial state." << std::endl; + + + //Set to true the monitoring prgress bar and the vacuum oscillations + nus_atm.Set_ProgressBar(true); + nus_atm.Set_IncludeOscillations(true); + + //Here we do the evolution of all the states + std::cout << "Begin: Evolution" << std::endl; + nus_atm.EvolveState(); + std::cout << "End: Evolution" << std::endl; + + //We can save the current state in HDF5 format for future use. + //nus_atm.WriteStateHDF5("./atmospheric_emission_example_numneu_"+std::to_string(numneu)+".hdf5"); + + //This file will contain the final flux, since initially we set it to 1 for the muon, + //this can be read as the muon ration F_final/F_initial in cos(zenith) and energy. + std::ofstream file("fluxes_flavor.txt"); + + + //Set the resolution and the ranges for the ouput, remember that an interpolation in energy is used in + //in the interaction picture, the vacuum oscillations are solve analytically with arbitrary Energy precision. + //For the zenith a linear interpolation is used. + int Nen=700; + int Ncz=100; + double lEmin=log10(Emin); + double lEmax=log10(Emax);; + + //Writing to the file! + file << "# log10(E) cos(zenith) E flux_i . . . ." << std::endl; + for(double cz=czmin;cz> plt; + if(plt=="yes" || plt=="y") + return system("./plot.plt"); + + + return 0; +} diff --git a/examples/Emitting_Atmosphere/plot.plt b/examples/Emitting_Atmosphere/plot.plt new file mode 100644 index 00000000..a4ccc5de --- /dev/null +++ b/examples/Emitting_Atmosphere/plot.plt @@ -0,0 +1,23 @@ +#!/usr/bin/env gnuplot +if ( GPVAL_VERSION >= 4.4 && strstrt(GPVAL_TERMINALS, 'wxt') > 0 ) set terminal wxt persist +if ( GPVAL_VERSION >= 4.4 && strstrt(GPVAL_TERMINALS, 'qt') > 0 ) set terminal qt persist +if ( GPVAL_VERSION >= 4.4 && strstrt(GPVAL_TERMINALS, 'wxt') == 0 && strstrt(GPVAL_TERMINALS, 'qt') == 0 ) print "wxt and qt terminals not available, proceeding with default" +if ( GPVAL_VERSION < 4.4 ) print "gnuplot is too old to check for available terminals" ; print "attempting to use wxt terminal and hoping for the best" ; set terminal wxt persist + +#set cbrange [0:1] +set logscale cb +set xlabel "log_{10}(E/GeV)" +set ylabel "Cos(zenith)" +set pm3d map +splot "fluxes_flavor.txt" u 1:2:5 + +set xrange [2.:6.] +set terminal png size 800,600 enhanced font 'Verdana,10' +set output "plot.png" +#set terminal postscript enhance eps color +#set output "plot.eps" +#set terminal svg size 350,262 fname 'Verdana' fsize 10 +#set output "plot.svg" +replot + + diff --git a/include/nuSQuIDS/body.h b/include/nuSQuIDS/body.h index 5f6b3db6..9f575731 100644 --- a/include/nuSQuIDS/body.h +++ b/include/nuSQuIDS/body.h @@ -669,6 +669,7 @@ class EarthAtm: public Body{ double x_ye_min; /// \brief Electron fraction at maximum radius. double x_ye_max; + public: /// \brief Default constructor using supplied PREM. EarthAtm(); @@ -699,6 +700,7 @@ class EarthAtm: public Body{ /// \brief EarthAtm trajectory class Track: public Body::Track{ friend class EarthAtm; + friend class EmittingEarthAtm; private: /// \brief Cosine of the zenith angle. double cosphi; @@ -762,6 +764,90 @@ class EarthAtm: public Body{ Track MakeTrackWithCosine(double cosphi); }; + + +/// \class EmittingEarthAtm +/// \brief A model of the production of neutrino cosmic ray flux production within the atmosphere. +class EmittingEarthAtm: public EarthAtm{ + protected: + + /// \brief Min/Max coszen. + double czen_min; + double czen_max; + /// \brief Min/Max height. + double height_min; + double height_max; + /// \brief Min/Max energy. + double energy_min; + double energy_max; + + /// \brief vector containing interpolators for the differential flux + std::vector> Interpolators; + + /// \brief indices for the electron and muon flavors + /// defaulted to the 0 and 1 states (first and second0 + /// WARNING!!!: careful when changing these as other parts of nuSQuIDS are hard coded to + /// default to 0, 1, & 2 for nu_e, nu_mu, and nu_tau respectively + /// Do this only if you have defined your own HI matrix by overriding squids::SU_vector HI + /// and if you have explicitly initialized all non-used mixing angles to 0 + int nu_e_index = 0; + int nu_mu_index = 1; + + public: + + /// \brief Default constructor using supplied MCEQ neutrino flux. + EmittingEarthAtm(); + + /// \brief Constructor from a user supplied production model. + /// @param prodmodel Path to the production model data file. + /// \details The input file should have 3+2n columns for n=the number of produced flavors. + /// The first one is cosine(zenith) in descending order from or between 1 and -1 + /// representing the range of coszen values corresponding to tracks relevant to your propagation + /// The second column must contain the height within the atmosphere in cm + /// in descending order from or between 6000000 (60 km) and 0 + /// The third column must correspond to the energy in GeV + /// in descending order over any relevant range + /// The remaining columns must correspond to the differential produced neutrino fluxes for each flavor + /// The order is assumed nu_e, nu_e_bar, nu_mu, nu_mu_bar with additional flavors added + /// with the corresponding antineutrino differential produced flux in the following column + /// the produced flux corresponds to the coszen, height, and energy in the row + EmittingEarthAtm(std::string prodmodel); + + /// \brief Deconstructor + ~EmittingEarthAtm(); + + /// \brief Function that sets the height of the atmosphere + /// \details Used by nuSQUIDSATM to set the height for every + /// instance of the body + void SetAtmosphereHeight(double height){ + atm_height = height; + earth_with_atm_radius = radius + atm_height; + } + /// \brief set index of produced flavors + /// \details Set the flavor index of the electron and muon neutrinos + /// or to different flavor indices as user input data requires + /// WARNING!!!: careful when using this as other parts of nuSQuIDS are hard coded to + /// default to 0, 1, & 2 for nu_e, nu_mu, and nu_tau respectively + /// Do this only if you have defined your own HI matrix by overriding squids::SU_vector HI + /// and if you have explicitly initialized all non-used mixing angles to 0 + void Set_ProducedFlavors(int nu_e_index_, int nu_mu_index_){ + nu_e_index = nu_e_index_; + nu_mu_index = nu_mu_index_; + } + + //Important function that overrides the default flux injection with ATM data + void injected_neutrino_flux(marray& flux, const GenericTrack& track, const nuSQUIDS& nusquids) override; + + /// \brief Returns the body identifier. + static unsigned int GetId() {return 8;} + /// \brief Returns the name of the body. + static std::string GetName() {return "EmittingEarthAtm";} +}; + + + + + // type defining typedef Body::Track Track; diff --git a/include/nuSQuIDS/nuSQuIDS.h b/include/nuSQuIDS/nuSQuIDS.h index a2ef4095..74d08731 100644 --- a/include/nuSQuIDS/nuSQuIDS.h +++ b/include/nuSQuIDS/nuSQuIDS.h @@ -1063,13 +1063,16 @@ class nuSQUIDS: public squids::SQuIDS { /** * The following class provides functionalities * for atmospheric neutrino experiments - * where a collection of trayectories is explored. + * where a collection of trajectories is explored. */ -template::value>::type > +template class nuSQUIDSAtm { + static_assert(std::is_base_of::value, "BaseNusType must be derived from nuSQUIDS"); + static_assert(std::is_base_of::value, "BaseBodyType must be derived from EarthAtm"); + // Class implementation here public: - using BaseSQUIDS = BaseType; + using BaseSQUIDS = BaseNusType; private: /// \brief Random number generator gsl_rng * r_gsl; @@ -1099,7 +1102,10 @@ class nuSQUIDSAtm { std::vector nusq_array; /// \brief Contains the Earth in atmospheric configuration. - std::shared_ptr earth_atm; + std::shared_ptr earth_atm; + + + /// \brief Contains the trajectories for each nuSQUIDS object, i.e. zenith. std::vector> track_array; /// \brief Contains the neutrino cross section object @@ -1131,18 +1137,20 @@ class nuSQUIDSAtm { const gsl_rng_type * T_gsl = gsl_rng_default; r_gsl = gsl_rng_alloc (T_gsl); - earth_atm = std::make_shared(); + earth_atm = std::make_shared(); + for(double costh : costh_array) track_array.push_back(std::make_shared(earth_atm->MakeTrackWithCosine(costh))); - + + for(unsigned int i = 0; i < costh_array.extent(0); i++){ nusq_array.emplace_back(args...); nusq_array.back().Set_Body(earth_atm); nusq_array.back().Set_Track(track_array[i]); } - ncs=nusq_array.front().GetNeutrinoCrossSections(); - + enu_array = nusq_array.front().GetERange(); + ncs=nusq_array.front().GetNeutrinoCrossSections(); log_enu_array.resize(std::vector{enu_array.size()}); //std::transform(enu_array.begin(), enu_array.end(), log_enu_array.begin(), // [](int enu) { return log(enu); }); @@ -1210,7 +1218,7 @@ class nuSQUIDSAtm { /************************************************************************************ * PUBLIC MEMBERS TO EVALUATE/SET/GET STUFF *************************************************************************************/ - + /// \brief Sets the initial state in the multiple energy mode /// when only considering neutrinos or antineutrinos /// @param ini_state Initial neutrino state. @@ -1272,6 +1280,16 @@ class nuSQUIDSAtm { } iinistate = true; } + + /// \brief Sets the index of the produced flavors. defaults to 0,1 for nuE and nuMu + void Set_AtmProducedFlavors(int nu_e_index_, int nu_mu_index_){ + earth_atm->Set_ProducedFlavors(nu_e_index_, nu_mu_index_); + } + + /// \brief Sets the height of the atmosphere for every body + void Set_AtmHeight(double height){ + earth_atm->SetAtmosphereHeight(height); + } /// \brief Evolves the system. void EvolveState(){ diff --git a/src/body.cpp b/src/body.cpp index 96d6a0b7..a873a474 100644 --- a/src/body.cpp +++ b/src/body.cpp @@ -29,6 +29,7 @@ #include #include #include +#include #include #include @@ -40,6 +41,7 @@ #include #include +#include // Macros #define SQR(x) ((x)*(x)) // x^2 @@ -116,6 +118,14 @@ unsigned int readUIntAttribute(hid_t object, std::string name){ return target; } +void unique_sort(std::vector& vec, double& min_val, double& max_val) { + std::sort(vec.begin(), vec.end()); + auto last = std::unique(vec.begin(), vec.end()); + vec.erase(last, vec.end()); + min_val = vec.front(); + max_val = vec.back(); +} + } // close unnamed namespace @@ -902,6 +912,158 @@ void EarthAtm::SetAtmosphereHeight(double height){ earth_with_atm_radius = radius + atm_height; } + + +/* +---------------------------------------------------------------------- + EMITTINGEARTHATM CLASS DEFINITIONS +---------------------------------------------------------------------- +*/ + +EmittingEarthAtm::EmittingEarthAtm():EmittingEarthAtm(getResourcePath()+"/atmos_prod/PROD_MODEL_MCEQ.dat") +{} + +EmittingEarthAtm::EmittingEarthAtm(std::string filepath):EarthAtm() +{ + + marray atm_prod_model = quickread(filepath); + arraysize = atm_prod_model.extent(0); + + std::vector diff_enu_prod(arraysize); + std::vector diff_munu_prod(arraysize); + + std::vector atm_coszen(arraysize); + std::vector atm_heights(arraysize); + std::vector atm_energy(arraysize); + + for (unsigned int i=0; i < arraysize;i++){ + atm_coszen[i] = atm_prod_model[i][0]; + atm_heights[i] = atm_prod_model[i][1]; + atm_energy[i] = atm_prod_model[i][2]; + diff_enu_prod[i] = atm_prod_model[i][3]; + diff_munu_prod[i] = atm_prod_model[i][4]; + } + + + unique_sort(atm_coszen, czen_min, czen_max); + long unsigned int czen_number = atm_coszen.size(); + marray coszens; + coszens.resize(0, czen_number); + for (unsigned int i=0; i < czen_number;i++){ + coszens[i] = atm_coszen[i]; + } + + unique_sort(atm_heights, height_min, height_max); + long unsigned int height_number = atm_heights.size(); + marray heights; + heights.resize(0, height_number); + for (unsigned int i=0; i < height_number; i++){ + heights[i] = atm_heights[i]; + } + + unique_sort(atm_energy, energy_min, energy_max); + long unsigned int energy_number = atm_energy.size(); + marray energies; + energies.resize(0, energy_number); + for (unsigned int i=0; i < energy_number; i++){ + energies[i] = atm_energy[i]; + } + + int num_prod_flav = (atm_prod_model.extent(1) - 3) / 2; + + std::vector>> flux_matrices(num_prod_flav, std::vector>(2)); + + marray nu_e_flux({czen_number, height_number, energy_number}); + marray nu_e_bar_flux({czen_number, height_number, energy_number}); + marray nu_mu_flux({czen_number, height_number, energy_number}); + marray nu_mu_bar_flux({czen_number, height_number, energy_number}); + + // Initialize each marray in flux_matrices + for (int flav = 0; flav < num_prod_flav; ++flav) { + for (int rho = 0; rho < 2; ++rho) { + flux_matrices[flav][rho] = marray({czen_number, height_number, energy_number}); + } + } + + + for (size_t i = 0; i < atm_prod_model.extent(0); ++i) { + double cz_ = atm_prod_model[i][0]; + double h_ = atm_prod_model[i][1]; + double E_ = atm_prod_model[i][2]; + // Find indices in coszens, heights, and energies + size_t czIndex = std::find(coszens.begin(), coszens.end(), cz_) - coszens.begin(); + size_t hIndex = std::find(heights.begin(), heights.end(), h_) - heights.begin(); + size_t EIndex = std::find(energies.begin(), energies.end(), E_) - energies.begin(); + + for (int flav = 0; flav < num_prod_flav; ++flav) { + for (int rho = 0; rho < 2; ++rho){ + flux_matrices[flav][rho][czIndex][hIndex][EIndex] = atm_prod_model[i][3+2*flav+rho]; + } + } + } + + Interpolators.resize(num_prod_flav); + for (int flav = 0; flav < num_prod_flav; ++flav) { + for (int rho = 0; rho < 2; ++rho) { + // Create a TriCubicInterpolator and add it to the current inner vector + Interpolators[flav].push_back(nusquids::TriCubicInterpolator(flux_matrices[flav][rho], energies, heights, coszens)); + } + } + + +} + + +// Overrides the neutrino flux injection +void EmittingEarthAtm::injected_neutrino_flux(marray& flux, const GenericTrack& track, const nuSQUIDS& nusquids) { + + + auto ESpace = nusquids.GetERange(); + // height calculation + const EarthAtm::Track& track_earthatm = static_cast(track); + double czNode = track_earthatm.cosphi; + double xkm = track.GetX() / param.km; + double sinsqphi = 1 - SQR(czNode); + double dL = sqrt(SQR(earth_with_atm_radius) - SQR(radius) * sinsqphi) + radius * czNode; + double L = sqrt(SQR(radius + atm_height) - SQR(radius) * sinsqphi) - radius * czNode; + double r2 = SQR(earth_with_atm_radius) + SQR(xkm) - (L / param.km + dL) * xkm; + double r = (r2 > 0 ? sqrt(r2) : 0); + double rel_r = r / earth_with_atm_radius; + double Height = earth_with_atm_radius * (rel_r - radius / earth_with_atm_radius)*param.km/param.cm; + + // Check if czNode and Height are within the required ranges + bool isZenithAndHeightInRange = (czNode >= czen_min && czNode <= czen_max) && + (Height >= height_min && Height <= height_max); + + int num_E = ESpace.extent(0); + for (unsigned int ei = 0; ei < num_E; ei++) { + + double Energy = ESpace[ei] / (param.GeV); + bool isInRange = isZenithAndHeightInRange && (Energy >= energy_min && Energy <= energy_max); + + //initializes all flavors' flux to 0 + for (int flav = 0; flav < nusquids.GetNumNeu(); ++flav) { + for (int rho = 0; rho < 2; ++rho){ + flux[ei][rho][flav] = 0; + } + } + + if (isInRange) { + for (int rho = 0; rho < 2; ++rho){ + flux[ei][rho][nu_e_index] = Interpolators[0][rho](Energy, Height, czNode); //energy passed in GeV Height passed in cm + if(Interpolators.size() > 1){ + flux[ei][rho][nu_mu_index] = Interpolators[1][rho](Energy, Height, czNode); + } + } + } + + } + +} + + +EmittingEarthAtm::~EmittingEarthAtm(){} + // body registration stuff std::map(hid_t)>>* body_registry=NULL; @@ -943,5 +1105,6 @@ NUSQUIDS_REGISTER_BODY(ConstantDensity); NUSQUIDS_REGISTER_BODY(VariableDensity); NUSQUIDS_REGISTER_BODY(Earth); NUSQUIDS_REGISTER_BODY(EarthAtm); +NUSQUIDS_REGISTER_BODY(EmittingEarthAtm); NUSQUIDS_REGISTER_BODY(Sun); NUSQUIDS_REGISTER_BODY(SunASnu);