From 1fdfc32725c7a337dd7461fe7b43cb276428eaf7 Mon Sep 17 00:00:00 2001 From: ConnorSponsler <33155897+ConnorSponsler@users.noreply.github.com> Date: Tue, 7 Nov 2023 16:58:34 -0500 Subject: [PATCH 01/12] Add files via upload --- src/body.cpp | 141 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 141 insertions(+) diff --git a/src/body.cpp b/src/body.cpp index 96d6a0b7..608fe656 100644 --- a/src/body.cpp +++ b/src/body.cpp @@ -29,6 +29,7 @@ #include #include #include +#include #include #include @@ -902,6 +903,145 @@ void EarthAtm::SetAtmosphereHeight(double height){ earth_with_atm_radius = radius + atm_height; } + + +/* +---------------------------------------------------------------------- + EMITTINGEARTHATM CLASS DEFINITIONS +---------------------------------------------------------------------- +*/ + +// constructor without passed energy, used by nuSQUIDSAtm +EmittingEarthAtm::EmittingEarthAtm():EmittingEarthAtm(getResourcePath()+"/atmos_prod/PROD_MODEL_MCEQ.dat") +{} + +// constructor with passed energy space +EmittingEarthAtm::EmittingEarthAtm(marray ESpace):EmittingEarthAtm(getResourcePath()+"/atmos_prod/PROD_MODEL_MCEQ.dat") +{ Set_EmissionEnergies(ESpace); } + + +EmittingEarthAtm::EmittingEarthAtm(std::string filepath):EarthAtm() +{ + + marray atm_prod_model = quickread(getResourcePath()+"atmos_prod/PROD_MODEL_MCEQ.dat"); + arraysize = atm_prod_model.extent(0); + + diff_enu_prod.resize(arraysize); + diff_munu_prod.resize(arraysize); + atm_coszen.resize(arraysize); + atm_heights.resize(arraysize); + atm_energy.resize(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, czen_number); + coszens.resize(0, atm_coszen.size()); + for (unsigned int i=0; i < czen_number;i++){ + coszens[i] = atm_coszen[i]; + } + + unique_sort(atm_heights, height_min, height_max, height_number); + heights.resize(0,atm_heights.size()); + for (unsigned int i=0; i < height_number; i++){ + heights[i] = atm_heights[i]; + } + + unique_sort(atm_energy, energy_min, energy_max, energy_number); + energies.resize(0,atm_energy.size()); + for (unsigned int i=0; i < energy_number; i++){ + energies[i] = atm_energy[i]; + } + + marray nuEmatrix{czen_number, height_number, energy_number}; + marray nuMumatrix{czen_number, height_number, energy_number}; + + + std::ifstream InputFile(filepath); + double cz_; + double h_; + double E_; + double nuEflux; + double nuMuflux; + + while (InputFile >> cz_ >> h_ >> E_ >> nuEflux >> nuMuflux) { + auto it1 = std::find(coszens.begin(), coszens.end(), cz_); + auto it2 = std::find(heights.begin(), heights.end(), h_); + auto it3 = std::find(energies.begin(), energies.end(), E_); + if (it1 != coszens.end() && it2 != heights.end() && it3 != energies.end()) { + size_t czIndex = std::distance(coszens.begin(), it1); + size_t hIndex = std::distance(heights.begin(), it2); + size_t EIndex = std::distance(energies.begin(), it3); + nuEmatrix[czIndex][hIndex][EIndex] = nuEflux; + nuMumatrix[czIndex][hIndex][EIndex] = nuMuflux; + } + } + InputFile.close(); + + Interpolators.push_back( nusquids::TriCubicInterpolator(nuEmatrix, energies, heights, coszens) ); + Interpolators.push_back( nusquids::TriCubicInterpolator(nuMumatrix, energies, heights, coszens) ); + +} + + +// Overrides the neutrino flux injection +void EmittingEarthAtm::injected_neutrino_flux(marray& flux, const GenericTrack& track, const nuSQUIDS& nusquids) { + + if(Get_EnergyInit() == false){ + throw std::runtime_error("Energy Space not Passed to EmittingEarthAtm instance"); + } + // 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); // here atm_height is a member of EarthAtm in km + double H = Height*100000; + + + for (unsigned int ei=0; ei < ESpace.extent(0); ei++) { + double E = ESpace[ei]/(param.GeV); + if(czNode >= czen_min && czNode<=czen_max && H>=height_min && H<=height_max && E>=energy_min && E<=energy_max){ + flux[ei][0][E_nu_index] = Interpolators[0](E, H, czNode); //(height*100000 passed in cm, Energy passed in GeV) + flux[ei][0][Mu_nu_index] = Interpolators[1](E, H, czNode); //(height*100000 passed in cm, Energy passed in GeV) + } else { + flux[ei][0][E_nu_index] = 0; //(height*100000 passed in cm, Energy passed in GeV) + flux[ei][0][Mu_nu_index] = 0; //(height*100000 passed in cm, Energy passed in GeV) + } + } + } + +EmittingEarthAtm::~EmittingEarthAtm(){} + + +void EmittingEarthAtm::unique_sort(std::vector& vec, double& min_val, double& max_val, long unsigned int& number_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(); + number_val = vec.size(); +} + +void Get_EnergyEmit(){ + +} + + + + + // body registration stuff std::map(hid_t)>>* body_registry=NULL; @@ -943,5 +1083,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); From 448ef9bb322930858c061b297e05a52b82ca3481 Mon Sep 17 00:00:00 2001 From: ConnorSponsler <33155897+ConnorSponsler@users.noreply.github.com> Date: Tue, 7 Nov 2023 16:59:40 -0500 Subject: [PATCH 02/12] Add files via upload --- include/nuSQuIDS/body.h | 121 ++++++++++++++++++++++++++++++++++++ include/nuSQuIDS/nuSQuIDS.h | 42 ++++++++++--- 2 files changed, 154 insertions(+), 9 deletions(-) diff --git a/include/nuSQuIDS/body.h b/include/nuSQuIDS/body.h index 5f6b3db6..1848595a 100644 --- a/include/nuSQuIDS/body.h +++ b/include/nuSQuIDS/body.h @@ -669,6 +669,9 @@ class EarthAtm: public Body{ double x_ye_min; /// \brief Electron fraction at maximum radius. double x_ye_max; + + marray ESpace; + public: /// \brief Default constructor using supplied PREM. EarthAtm(); @@ -699,6 +702,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 +766,123 @@ 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 Atm differential nu electron flux array + std::vector diff_enu_prod; + marray enu_flux; + /// \brief Atm differential nu muon flux array + std::vector diff_munu_prod; + marray munu_flux; + /// \brief Atm coszen array + std::vector atm_coszen; + marray coszens; + /// \brief Atm heights array + std::vector atm_heights; + marray heights; + /// \brief Atm energies array + std::vector atm_energy; + marray energies; + + /// \brief Minimum coszen. + double czen_min; + /// \brief Maximum coszen. + double czen_max; + /// \brief number of coszens + long unsigned int czen_number; + /// \brief Minimum height. + double height_min; + /// \brief Maximum height. + double height_max; + /// \brief number of heights + long unsigned int height_number; + /// \brief Minimum energy. + double energy_min; + /// \brief Maximum energy. + double energy_max; + /// \brief number of energies + long unsigned int energy_number; + /// \brief vector containing interpolators for the differential flux + std::vector Interpolators; + /// \brief indices for the electron and muon flavors + /// defaulted to the first and second flavor states + int E_nu_index = 0; + int Mu_nu_index = 1; + + bool init_emit_energies = false; + + public: + + /// \brief Default constructor using supplied MCEQ neutrino flux. + EmittingEarthAtm(); + + /// \brief constructor with user passed energy space + EmittingEarthAtm(marray ESpace); + + /// \brief Constructor from a user supplied production model. + /// @param prodmodel Path to the production model data file. + /// \details The input file should have five columns. + /// The first one must run 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 3rd column must correspond to the energy in GeV + /// in descending order over any relevant range + /// The 4th and 5th columns must correspond to the differential E_nu and Mu_nu + /// flux corresponding to the coszen, height, and energy in the row + EmittingEarthAtm(std::string prodmodel); + + /// \brief Deconstructor + ~EmittingEarthAtm(); + + /// \brief Function that passes the user's list of energies to the class + /// \details Used by nuSQUIDSATM implementation to provide the + /// class with the energy list required for retrieval of interpolation + void Set_EmissionEnergies(const marray& ESpace_){ + ESpace = ESpace_; + init_emit_energies = true; + } + /// \brief Returns true if neutrino flux emission from bodies is considered + bool Get_EnergyInit(){ + return init_emit_energies; + } + /// \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 + void Set_ProducedFlavors(int E_nu_index_, int Mu_nu_index_){ + E_nu_index = E_nu_index_; + Mu_nu_index = Mu_nu_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; + + //Helpful function for data parsing + void unique_sort(std::vector& vec, double& min_val, double& max_val, long unsigned int& number_val); + + + /// \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..04303548 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,22 @@ class nuSQUIDSAtm { } iinistate = true; } + + /// \brief Passes the log scale energy array into the user constructed body + void Set_AtmEmissionEnergies(){ + earth_atm->Set_EmissionEnergies(enu_array); + } + + + /// \brief Sets the index of the produced flavors. defaults to 0,1 for nuE and nuMu + void Set_AtmProducedFlavors(int E_nu_index_, int Mu_nu_index_){ + earth_atm->Set_ProducedFlavors(E_nu_index_, Mu_nu_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(){ From 58aeed1019aee7a0b80c49e0266dda2dcf25aad0 Mon Sep 17 00:00:00 2001 From: ConnorSponsler <33155897+ConnorSponsler@users.noreply.github.com> Date: Tue, 7 Nov 2023 17:00:50 -0500 Subject: [PATCH 03/12] Add files via upload --- examples/Emitting_Atmosphere/main.cpp | 189 ++++++++++++++++++++++++++ examples/Emitting_Atmosphere/plot.plt | 23 ++++ 2 files changed, 212 insertions(+) create mode 100644 examples/Emitting_Atmosphere/main.cpp create mode 100644 examples/Emitting_Atmosphere/plot.plt 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 + + From 9523342a61085ff2ebcc8fbe057b5f1b5576b4a7 Mon Sep 17 00:00:00 2001 From: ConnorSponsler <33155897+ConnorSponsler@users.noreply.github.com> Date: Tue, 7 Nov 2023 17:03:51 -0500 Subject: [PATCH 04/12] Update configure with new example --- configure | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/configure b/configure index a165ef43..edbae8cb 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 $@ From 375655dbf277324032fb58f9e539af257f43bb9c Mon Sep 17 00:00:00 2001 From: ConnorSponsler <33155897+ConnorSponsler@users.noreply.github.com> Date: Fri, 17 Nov 2023 20:30:48 -0500 Subject: [PATCH 05/12] Fix allignment of new example --- configure | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/configure b/configure index edbae8cb..b08e6eea 100755 --- a/configure +++ b/configure @@ -709,7 +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/Emitting_Atmosphere/emitting_atmosphere \ examples/Bodies/bodies \ examples/Xsections/xsections \ examples/NSI/nsi \ From dc543406443a2a647f6c5aba8f41c197b3e400d8 Mon Sep 17 00:00:00 2001 From: ConnorSponsler <33155897+ConnorSponsler@users.noreply.github.com> Date: Sun, 19 Nov 2023 16:11:59 -0500 Subject: [PATCH 06/12] Add files via upload adjust EmittingEarthAtm implementation --- src/body.cpp | 194 ++++++++++++++++++++++++++++++--------------------- 1 file changed, 113 insertions(+), 81 deletions(-) diff --git a/src/body.cpp b/src/body.cpp index 608fe656..70b78f31 100644 --- a/src/body.cpp +++ b/src/body.cpp @@ -41,6 +41,7 @@ #include #include +#include // Macros #define SQR(x) ((x)*(x)) // x^2 @@ -117,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 @@ -911,26 +920,21 @@ void EarthAtm::SetAtmosphereHeight(double height){ ---------------------------------------------------------------------- */ -// constructor without passed energy, used by nuSQUIDSAtm EmittingEarthAtm::EmittingEarthAtm():EmittingEarthAtm(getResourcePath()+"/atmos_prod/PROD_MODEL_MCEQ.dat") {} -// constructor with passed energy space -EmittingEarthAtm::EmittingEarthAtm(marray ESpace):EmittingEarthAtm(getResourcePath()+"/atmos_prod/PROD_MODEL_MCEQ.dat") -{ Set_EmissionEnergies(ESpace); } - - EmittingEarthAtm::EmittingEarthAtm(std::string filepath):EarthAtm() { - marray atm_prod_model = quickread(getResourcePath()+"atmos_prod/PROD_MODEL_MCEQ.dat"); + marray atm_prod_model = quickread(filepath); arraysize = atm_prod_model.extent(0); - diff_enu_prod.resize(arraysize); - diff_munu_prod.resize(arraysize); - atm_coszen.resize(arraysize); - atm_heights.resize(arraysize); - atm_energy.resize(arraysize); + 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]; @@ -941,106 +945,134 @@ EmittingEarthAtm::EmittingEarthAtm(std::string filepath):EarthAtm() } - unique_sort(atm_coszen, czen_min, czen_max, czen_number); - coszens.resize(0, atm_coszen.size()); + unique_sort(atm_coszen, czen_min, czen_max); + 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, height_number); - heights.resize(0,atm_heights.size()); + unique_sort(atm_heights, height_min, height_max); + 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, energy_number); - energies.resize(0,atm_energy.size()); + unique_sort(atm_energy, energy_min, energy_max); + 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]; } - - marray nuEmatrix{czen_number, height_number, energy_number}; - marray nuMumatrix{czen_number, height_number, energy_number}; - - std::ifstream InputFile(filepath); - double cz_; - double h_; - double E_; - double nuEflux; - double nuMuflux; - - while (InputFile >> cz_ >> h_ >> E_ >> nuEflux >> nuMuflux) { - auto it1 = std::find(coszens.begin(), coszens.end(), cz_); - auto it2 = std::find(heights.begin(), heights.end(), h_); - auto it3 = std::find(energies.begin(), energies.end(), E_); - if (it1 != coszens.end() && it2 != heights.end() && it3 != energies.end()) { - size_t czIndex = std::distance(coszens.begin(), it1); - size_t hIndex = std::distance(heights.begin(), it2); - size_t EIndex = std::distance(energies.begin(), it3); - nuEmatrix[czIndex][hIndex][EIndex] = nuEflux; - nuMumatrix[czIndex][hIndex][EIndex] = nuMuflux; + 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({static_cast(czen_number), + static_cast(height_number), + static_cast(energy_number)}); + marray nu_e_bar_flux({static_cast(czen_number), + static_cast(height_number), + static_cast(energy_number)}); + marray nu_mu_flux({static_cast(czen_number), + static_cast(height_number), + static_cast(energy_number)}); + marray nu_mu_bar_flux({static_cast(czen_number), + static_cast(height_number), + static_cast(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({static_cast(czen_number), + static_cast(height_number), + static_cast(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]; } } - InputFile.close(); + } - Interpolators.push_back( nusquids::TriCubicInterpolator(nuEmatrix, energies, heights, coszens) ); - Interpolators.push_back( nusquids::TriCubicInterpolator(nuMumatrix, energies, heights, coszens) ); + 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) { + - if(Get_EnergyInit() == false){ - throw std::runtime_error("Energy Space not Passed to EmittingEarthAtm instance"); + 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; } - // 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); // here atm_height is a member of EarthAtm in km - double H = Height*100000; - - - for (unsigned int ei=0; ei < ESpace.extent(0); ei++) { - double E = ESpace[ei]/(param.GeV); - if(czNode >= czen_min && czNode<=czen_max && H>=height_min && H<=height_max && E>=energy_min && E<=energy_max){ - flux[ei][0][E_nu_index] = Interpolators[0](E, H, czNode); //(height*100000 passed in cm, Energy passed in GeV) - flux[ei][0][Mu_nu_index] = Interpolators[1](E, H, czNode); //(height*100000 passed in cm, Energy passed in GeV) - } else { - flux[ei][0][E_nu_index] = 0; //(height*100000 passed in cm, Energy passed in GeV) - flux[ei][0][Mu_nu_index] = 0; //(height*100000 passed in cm, Energy passed in GeV) + } + + 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(){} - - -void EmittingEarthAtm::unique_sort(std::vector& vec, double& min_val, double& max_val, long unsigned int& number_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(); - number_val = vec.size(); -} + } -void Get_EnergyEmit(){ - } - - +EmittingEarthAtm::~EmittingEarthAtm(){} // body registration stuff From 27a6e9ceef798fa95823b2344df6b0c691c14669 Mon Sep 17 00:00:00 2001 From: ConnorSponsler <33155897+ConnorSponsler@users.noreply.github.com> Date: Sun, 19 Nov 2023 16:20:26 -0500 Subject: [PATCH 07/12] Add files via upload adjustments to EmittingEarthAtm implementation --- body.h | 873 +++++++++++++++++++++++ nuSQuIDS.h | 1951 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 2824 insertions(+) create mode 100644 body.h create mode 100644 nuSQuIDS.h diff --git a/body.h b/body.h new file mode 100644 index 00000000..c35d7aeb --- /dev/null +++ b/body.h @@ -0,0 +1,873 @@ + /****************************************************************************** + * 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 * + ******************************************************************************/ + + +#ifndef NUSQUIDS_BODY_H +#define NUSQUIDS_BODY_H + +#if __cplusplus < 201103L +#error C++11 compiler required. Update your compiler and use the flag -std=c++11 +#endif + +#include +#include +#include +#include +#include +#include +#include "nuSQuIDS/marray.h" +#include "nuSQuIDS/tools.h" + +namespace nusquids{ + +// Forward declaration of nuSQuIDS class +class nuSQUIDS; + +/// \class Body +/// \brief Abstract body class. +/// \details This abstract class serves as a prototype +/// of neutrino propagation environments. When implementing +/// the subclasses one must define the density() and ye() +class Body{ + protected: + /// \brief Body parameters. + std::vector BodyParams; + /// \brief bool that signal if the body is constant density or not + bool is_constant_density = false; + public: + /// \brief Body Constructor; + Body(){} + virtual ~Body(){} + + /// \brief Serialization function + virtual void Serialize(hid_t group) const=0; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + + /// \class Track + /// \brief Trajectory subclass. Specifies the trajectory + /// to take within the body. + class Track{ + friend class Body; + protected: + /// \brief Current position. + double x; + /// \brief Initial position. + double xini; + /// \brief Final position. + double xend; + public: + /// \brief Default constructor. + Track(double x,double xini, double xend): x(x), xini(xini),xend(xend) {} + /// \brief Default constructor. + Track(double xini, double xend): Track(xini,xini,xend) {} + /// \brief Serialization function + virtual void Serialize(hid_t group) const=0; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + /// \brief Get track object name + static std::string GetName() {return "BodyTrack";}; + /// \brief Sets the current position along the trajectory. + void SetX(double y){ x = y; } + /// \brief Returns the current position along the trajectory. + double GetX() const { return x; } + /// \brief Returns the initial position measured along the trajectory. + double GetInitialX() const { return xini; } + /// \brief Returns the final position measured along the trajectory. + double GetFinalX() const { return xend; } + /// \brief Reverses the track initial and final position + void ReverseTrack() { + double xtmp = xend; + xend = xini; + xini = xtmp; + } + /// \brief Returns parameters that define the trajectory. + std::vector GetTrackParams() const { + std::vector TrackParams{xini,xend}; + FillDerivedParams(TrackParams); + return TrackParams; + } + /// Should be implemented by derived classes to append their + /// additional parameters to TrackParams + virtual void FillDerivedParams(std::vector& TrackParams) const{} + }; + /// \brief Return the density at a given trajectory object. + virtual double density(const Track&) const {return 0.0;} + /// \brief Return the electron fraction at a given trajectory object. + virtual double ye(const Track&) const {return 1.0;} + /// \brief Returns parameters that define the body. + const std::vector& GetBodyParams() const { return BodyParams;} + /// \brief Returns the body identifier. + static unsigned int GetId() {return 0;} + /// \brief Returns the name of the body. + static std::string GetName() {return "Body";} + /// \brief Return true if the body is a constant density. + virtual bool IsConstantDensity() const {return is_constant_density;} + /// \brief Return true if the body is a constant density. + virtual void SetIsConstantDensity(bool icd) {is_constant_density = icd;} + /// \brief Sets the injected neutrino flux from the Body at a given Track location. + /// \details This function is called by nuSQuIDS during the calculation to obtain the injected flux. + /// @param flux Is a three dimensional array with dimensions: energy, rho, and flavor, which has previously been allocated by nuSQuIDS. + /// @param Track Trajectory object. + /// @param nuSQuIDS nuSQuIDS object that queries the flux. + virtual void injected_neutrino_flux(marray& flux, const Track&, const nuSQUIDS&) {} +}; + +// type defining +typedef Body::Track GenericTrack; + +/// \class Vacuum +/// \brief Vacuum body. +/// \details A zero density body. +class Vacuum: public Body { + public: + /// \brief Constructor. + Vacuum():Body(){} + + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + + /// \class Track + /// \brief Vacuum trajectory + class Track: public Body::Track { + public : + /// \brief Rectilinear trajectory from xini to xend. + /// @param x current position [eV^-1]. + /// @param xini Initial position [eV^-1]. + /// @param xend Final position [eV^-1]. + Track(double x,double xini,double xend):Body::Track(x,xini,xend){} + /// \brief Rectilinear trajectory from xini to xend. + /// @param xini Initial position [eV^-1]. + /// @param xend Final position [eV^-1]. + Track(double xini,double xend):Track(xini,xini,xend){} + /// \brief Construct a trajectory of distance xend. + /// @param xend Final position in eV^-1. + /// \details In this case initial position is assumed 0. + Track(double xend):Track(0.0,xend){} + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + /// \brief Get track object name + static std::string GetName() {return "VacuumTrack";} + }; + /// \brief Returns the body identifier. + static unsigned int GetId() {return 1;} + /// \brief Returns the name of the body. + static std::string GetName() {return "Vacuum";} + /// \brief Returns the density in g/cm^3 + double density(const GenericTrack&) const; + /// \brief Returns the electron fraction + double ye(const GenericTrack&) const; + /// \brief Returns true as this body has constant density by definition + bool IsConstantDensity() const; +}; + +/// \class ConstantDensity +/// \brief Constant density body. +/// \details A body with constant density +class ConstantDensity: public Body{ + private: + /// \brief Constant density in g/cm^3 + const double constant_density; + /// \brief Constant electron fraction + const double constant_ye; + public: + /// \brief Constructor + /// @param density Density in g/cm^3. + /// @param ye electron fraction. + ConstantDensity(double density,double ye); + + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + + /// \class Track + /// \brief Constant density trajectory + class Track: public Body::Track{ + public : + /// \brief Rectilinear trajectory from xini to xend. + /// @param x current position [eV^-1]. + /// @param xini Initial position [eV^-1]. + /// @param xend Final position [eV^-1]. + Track(double x,double xini,double xend):Body::Track(x,xini,xend){} + /// \brief Construct a trajectory between an initial position and final position. + /// @param xini Initial position in eV^-1. + /// @param xend Final position in eV^-1. + Track(double xini,double xend):Track(xini,xini,xend){} + /// \brief Construct a trajectory of distance xend. + /// @param xend Final position in eV^-1. + /// \details In this case initial position is assumed 0. + Track(double xend):Track(0.0,xend){} + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + /// \brief Get track object name + static std::string GetName() {return "ConstantDensityTrack";} + }; + + /// \brief Returns the body identifier. + static unsigned int GetId() {return 2;} + /// \brief Returns the name of the body. + static std::string GetName() {return "ConstantDensity";} + /// \brief Returns the density in g/cm^3 + double density(const GenericTrack&) const; + /// \brief Returns the electron fraction + double ye(const GenericTrack&) const; + /// \brief Returns true as this body has constant density by definition + bool IsConstantDensity() const; +}; + +/// \class VariableDensity +/// \brief Variable user defined density body. +/// \details The user provides arrays which contain the body +/// information. +class VariableDensity: public Body{ + private: + /// \brief Array of length \c arraysize containing spline position nodes. + std::vector x_arr; + /// \brief Array of length \c arraysize containing the density at each position nodes. + std::vector density_arr; + /// \brief Array of length \c arraysize containing the electron fraction at each position nodes. + std::vector ye_arr; + /// \brief Size of array used to create spline. + unsigned int arraysize; + + /// \brief Minimum value of \c x_arr. + double x_max; + /// \brief Maximum value of \c x_arr. + double x_min; + + /// \brief Density spline + AkimaSpline inter_density; + /// \brief Electron fraction spline + AkimaSpline inter_ye; + public: + /// \brief Constructor. + /// @param x Vector containing position nodes in cm. + /// @param density Density, in g/cm^3, at each of the nodes. + /// @param ye Electron fraction at each of the nodes. + /// \pre All input vectors must be of equal size. + VariableDensity(std::vector x,std::vector density,std::vector ye); + /// \brief Destructor. + ~VariableDensity(); + + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + + /// \class Track + /// \brief Variable density trajectory + class Track: public Body::Track{ + public : + /// \brief Rectilinear trajectory from xini to xend. + /// @param x current position [eV^-1]. + /// @param xini Initial position [eV^-1]. + /// @param xend Final position [eV^-1]. + Track(double x,double xini,double xend):Body::Track(x,xini,xend){} + /// \brief Construct a trajectory between an initial position and final position. + /// @param xini Initial position in eV^-1. + /// @param xend Final position in eV^-1. + Track(double xini,double xend):Track(xini,xini,xend){} + /// \brief Construct a trajectory of distance xend. + /// @param xend Final position in eV^-1. + /// \details In this case initial position is assumed 0. + Track(double xend):Track(0.0,xend){} + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + /// \brief Get track object name + static std::string GetName() {return "VariableDensityTrack";} + + }; + /// \brief Returns the body identifier. + static unsigned int GetId() {return 3;} + /// \brief Returns the name of the body. + static std::string GetName() {return "VariableDensity";} + /// \brief Returns the density in g/cm^3 + double density(const GenericTrack&) const; + /// \brief Returns the electron fraction + double ye(const GenericTrack&) const; +}; + +/// \class Earth +/// \brief Earth model based on PREM. +class Earth: public Body{ + private: + /// \brief Radius of the Earth. + double radius; + /// \brief Earth radius position array + std::vector earth_radius; + /// \brief Earth density array + std::vector earth_density; + /// \brief Earth electron fraction array + std::vector earth_ye; + /// \brief Data arrays size + unsigned int arraysize; + + /// \brief Density spline + AkimaSpline inter_density; + /// \brief Electron fraction spline + AkimaSpline inter_ye; + + /// \brief Minimum radius. + double x_radius_min; + /// \brief Maximum radius. + double x_radius_max; + /// \brief Density at minimum radius. + double x_rho_min; + /// \brief Density at maximum radius. + double x_rho_max; + /// \brief Electron fraction at minimum radius. + double x_ye_min; + /// \brief Electron fraction at maximum radius. + double x_ye_max; + public: + /// \brief Default constructor using supplied PREM. + Earth(); + /// \brief Constructor from a user supplied Earth model. + /// @param earthmodel Path to the Earth model file. + /// \details The input file should have three columns. + /// The first one must run from zero to one representing + /// the center and surface of the Earth respectively. The + /// second column must contain the Earth density in g/cm^3 at + /// a given position, while the third column must contain + /// the electron fraction. + Earth(std::string earthmodel); + /// \brief Constructor in which the user provides, as vectors, the + /// Earth properties. + /// @param x Vector containing position nodes in cm. + /// @param rho Density, in g/cm^3, at each of the nodes. + /// @param ye Electron fraction at each of the nodes. + /// \pre All input vectors must be of equal size. + Earth(std::vector x,std::vector rho,std::vector ye); + /// \brief Destructor. + ~Earth(); + + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + + /// \class Track + /// \brief Earth trajectory + class Track: public Body::Track{ + private: + /// \brief Baseline of the neutrino experiment in natural units. + const double baseline; + public : + /// \brief Rectilinear trajectory from xini to xend. + /// @param x current position [eV^-1]. + /// @param xini Initial position [eV^-1]. + /// @param xend Final position [eV^-1]. + Track(double x,double xini,double xend,double baseline):Body::Track(x,xini,xend),baseline(baseline){} + /// \brief Construct a trajectory between an initial position and final position. + /// @param xini Initial position in eV^-1. + /// @param xend Final position in eV^-1. + /// @param baseline Baseline of experiment in eV^-1. + Track(double xini,double xend,double baseline):Track(xini,xini,xend,baseline){} + /// \brief Construct a trajectory where the neutrino travels a distance given by baseline. + /// @param baseline Traverse distance in eV^-1. + /// \details In this case \c xini = 0, \c xend = \c baseline. + Track(double baseline):Track(0.,baseline,baseline){} + /// \brief Returns the neutrino baseline in natural units. + double GetBaseline() const {return baseline;} + virtual void FillDerivedParams(std::vector& TrackParams) const; + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + /// \brief Get track object name + static std::string GetName() {return "EarthTrack";} + }; + /// \brief Returns the body identifier. + static unsigned int GetId() {return 4;} + + /// \brief Returns the name of the body. + static std::string GetName() {return "Earth";} + + /// \brief Returns the density in g/cm^3 + double density(const GenericTrack&) const; + /// \brief Returns the electron fraction + double ye(const GenericTrack&) const; + + /// \brief Returns the radius of the Earth in natural units. + double GetRadius() const {return radius;} +}; + +/// \class Sun +/// \brief A model of the Sun. +class Sun: public Body{ + private: + /// \brief Array that contains all the Solar model parameters. + marray sun_model; + /// \brief Array that contains all the Solar electron density parameters. + marray sun_model_nele; + /// \brief Array of length \c arraysize containing spline position nodes. + std::vector sun_radius; + /// \brief Array of length \c arraysize containing the density at position nodes. + std::vector sun_density; + /// \brief Array of length \c arraysize containing the hydrogen fraction at position nodes. + std::vector sun_xh; + + /// \brief Size of \c sun_radius array. + unsigned int arraysize; + + /// \brief Radius of the Sun. + double radius; + + /// \brief Density spline + AkimaSpline inter_density; + /// \brief Hydrogen fraction spline + AkimaSpline inter_xh; + + /// \brief Returns the density in g/cm^3 at a given radius fraction x + /// @param x Radius fraction: 0:center, 1:surface. + double rdensity(double x) const; + /// \brief Returns the electron fraction at a given radius fraction x + /// @param x Radius fraction: 0:center, 1:surface. + double rxh(double x) const; + public: + /// \brief Detault constructor. + Sun(); + /// \brief Constructor in which the user provides, as vectors, the + /// Sun properties. + /// @param x Vector containing position nodes in cm. + /// @param rho Density, in g/cm^3, at each of the nodes. + /// @param xh Hydrogen fraction at each point. + /// \pre All input vectors must be of equal size. + Sun(std::vector x,std::vector rho,std::vector xh); + /// \brief Constructor from a user supplied solar model. + /// @param sunmodel Path to the Sun model file. + /// \details The input file should have the same columns as John Bahcalls solar model file. + /// The second column one must run from zero to one representing + /// the center and surface of the Sun respectively. The + /// fourth column must contain the Sun density in g/cm^3 at + /// a given position, while the seven column must contain + /// the hydrogen fraction which is related to the electron fraction by ye = 0.5*(1.0+rxh(r)). + /// Internally we assume that the solar radius is 695980.0 kilometers. + Sun(std::string sunmodel); + + + /// \brief Destructor + ~Sun(); + + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + + /// \class Track + /// \brief Sun trajectory + class Track: public Body::Track{ + public : + /// \brief Rectilinear trajectory from xini to xend. + /// @param x current position [eV^-1]. + /// @param xini Initial position [eV^-1]. + /// @param xend Final position [eV^-1]. + Track(double x,double xini,double xend):Body::Track(x,xini,xend){} + /// \brief Construct a trajectory between an initial position and final position. + /// @param xini Initial position in eV^-1. + /// @param xend Final position in eV^-1. + /// \details The trajectory is measured from the sun center which is set to zero. + Track(double xini,double xend):Track(xini,xini,xend){} + /// \brief Construct a trajectory from the sun center to a final position. + /// @param xend Final position in eV^-1. + /// \details The trajectory is measured from the sun center which is set to zero. + Track(double xend):Track(0.,xend){} + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + /// \brief Get track object name + static std::string GetName() {return "SunTrack";} + }; + + /// \brief Returns the body identifier. + static unsigned int GetId() {return 5;} + + /// \brief Returns the name of the body. + static std::string GetName() {return "Sun";} + + /// \brief Returns the density in g/cm^3 + double density(const GenericTrack&) const; + /// \brief Returns the electron fraction + double ye(const GenericTrack&) const; + + /// \brief Returns the radius of the Sun in natural units. + double GetRadius() const {return radius;} +}; + +/// \class SunASnu +/// \brief A model of the Sun with atmospheric solar neutrinos geometry. +class SunASnu: public Body{ + private: + /// \brief Array that contains all the Solar model parameters. + marray sun_model; + /// \brief Array of length \c arraysize containing spline position nodes. + std::vector sun_radius; + /// \brief Array of length \c arraysize containing the density at position nodes. + std::vector sun_density; + /// \brief Array of length \c arraysize containing the hydrogen fraction at position nodes. + std::vector sun_xh; + + /// \brief Size of \c sun_radius array. + unsigned int arraysize; + + /// \brief Radius of the Sun. + double radius; + + /// \brief Density spline + AkimaSpline inter_density; + /// \brief Hydrogen fraction spline + AkimaSpline inter_xh; + + /// \brief Returns the density in g/cm^3 at a given radius fraction x + /// @param x Radius fraction: 0:center, 1:surface. + double rdensity(double x) const; + /// \brief Returns the electron fraction at a given radius fraction x + /// @param x Radius fraction: 0:center, 1:surface. + double rxh(double x) const; + public: + /// \brief Detault constructor. + SunASnu(); + /// \brief Constructor from a user supplied solar model. + /// @param sunmodel Path to the Sun model file. + /// \details The input file should have the same columns as John Bahcalls solar model file. + /// The second column one must run from zero to one representing + /// the center and surface of the Sun respectively. The + /// fourth column must contain the Sun density in g/cm^3 at + /// a given position, while the seven column must contain + /// the hydrogen fraction which is related to the electron fraction by ye = 0.5*(1.0+rxh(r)). + /// Internally we assume that the solar radius is 695980.0 kilometers. + SunASnu(std::string sunmodel); + /// \brief Constructor in which the user provides, as vectors, the + /// Sun properties. + /// @param x Vector containing position nodes in cm. + /// @param rho Density, in g/cm^3, at each of the nodes. + /// @param xh Hydrogen fraction at each point. + /// \pre All input vectors must be of equal size. + SunASnu(std::vector x,std::vector rho,std::vector xh); + + ~SunASnu(); + + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + + /// \class Track + /// \brief SunASnu trajectory + class Track: public Body::Track{ + friend class SunASnu; + private: + /// \brief Radius of the Sun. + double radius_nu; + /// \brief Impact parameter. + double b_impact; + public: + /// \brief Rectilinear trajectory from xini to xend. + /// @param x current position [eV^-1]. + /// @param xini Initial position [eV^-1]. + /// @param b_impact impact parameter in eV^-1. + Track(double x,double xini,double b_impact); + /// \brief Construct a trajectory between an initial position and final position. + /// @param xini Initial position in eV^-1. + /// @param b_impact impact parameter in eV^-1. + /// \details The trajectory baseline is determined by the impact parameter and starts + /// at \c xini, and ends when the neutrino exits the sun. + Track(double xini,double b_impact):Track(xini,xini,b_impact){} + /// \brief Construct a for a given impact parameter. + /// @param b_impact_ impact parameter in eV^-1. + /// \details The trajectory baseline is determined by the impact parameter and starts + /// at \c xini = 0, and ends when the neutrino exits the sun. + Track(double b_impact_):Track(0.0,b_impact_){} + virtual void FillDerivedParams(std::vector& TrackParams) const; + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + /// \brief Get track object name + static std::string GetName() {return "SunASnuTrack";} + }; + + /// \brief Returns the body identifier. + static unsigned int GetId() {return 6;} + + /// \brief Returns the name of the body. + static std::string GetName() {return "SunASnu";} + + /// \brief Returns the density in g/cm^3 + double density(const GenericTrack&) const; + /// \brief Returns the electron fraction + double ye(const GenericTrack&) const; + + /// \brief Returns the radius of the Sun in natural units. + double GetRadius() const {return radius;} +}; + +/// \class EarthAtm +/// \brief A model of the Earth with atmospheric neutrinos geometry. +class EarthAtm: public Body{ + protected: + /// \brief Radius of the Earth. + double radius; + /// \brief Height of the atmosphere. + double atm_height; + /// \brief Radius of the Earth plus atmosphere. + double earth_with_atm_radius; + + /// \brief Earth radius position array + std::vector earth_radius; + /// \brief Earth density array + std::vector earth_density; + /// \brief Earth electron fraction array + std::vector earth_ye; + /// \brief Data arrays size + unsigned int arraysize; + + /// \brief Density spline + AkimaSpline inter_density; + /// \brief Electron fraction spline + AkimaSpline inter_ye; + + /// \brief Minimum radius. + double x_radius_min; + /// \brief Maximum radius. + double x_radius_max; + /// \brief Density at minimum radius. + double x_rho_min; + /// \brief Density at maximum radius. + double x_rho_max; + /// \brief Electron fraction at minimum radius. + double x_ye_min; + /// \brief Electron fraction at maximum radius. + double x_ye_max; + + public: + /// \brief Default constructor using supplied PREM. + EarthAtm(); + /// \brief Constructor from a user supplied Earth model. + /// @param earthmodel Path to the Earth model file. + /// \details The input file should have three columns. + /// The first one must run from zero to one representing + /// the center and surface of the Earth respectively. The + /// second column must contain the Earth density in g/cm^3 at + /// a given position, while the third column must contain + /// the electron fraction. + EarthAtm(std::string earthmodel); + /// \brief Constructor in which the user provides, as vectors, the + /// Earth properties. + /// @param x Vector containing position nodes in cm. + /// @param rho Density, in g/cm^3, at each of the nodes. + /// @param ye Electron fraction at each of the nodes. + /// \pre All input vectors must be of equal size. + EarthAtm(std::vector x,std::vector rho,std::vector ye); + ~EarthAtm(); + + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + + /// \class Track + /// \brief EarthAtm trajectory + class Track: public Body::Track{ + friend class EarthAtm; + friend class EmittingEarthAtm; + private: + /// \brief Cosine of the zenith angle. + double cosphi; + /// \brief Radius of the Earth. + double earth_radius; + /// \brief Height of the atmosphere. + double atmheight; + /// \brief Baseline. + double L; + + ///Private constructor which does not initialize all members. + ///Intended for use only by makeWithCosine. + Track(); + public : + /// \brief Construct trajectory. + /// @param x_ current position [eV^-1]. + /// @param phi Zenith angle in radians. + Track(double x_,double phi,double earth_radius_,double atmheight_): + Track(phi,earth_radius_,atmheight_){x=x_; assert(x >= 0);}; + /// \brief Construct trajectory. + /// @param phi Zenith angle in radians. + Track(double phi,double earth_radius_,double atmheight_); + /// \brief Returns the neutrino baseline in natural units. + double GetBaseline() const {return L;} + virtual void FillDerivedParams(std::vector& TrackParams) const; + ///Construct a track with the cosine of the zenith angle + static Track MakeWithCosine(double cosphi,double earth_radius_,double atmheight_); + /// \brief Serialization function + void Serialize(hid_t group) const; + /// \brief Deserialization function + static std::shared_ptr Deserialize(hid_t group); + /// \brief Get track object name + static std::string GetName() {return "EarthAtmTrack";} + }; + /// \brief Returns the body identifier. + static unsigned int GetId() {return 7;} + /// \brief Returns the name of the body. + static std::string GetName() {return "EarthAtm";} + /// \brief Returns the density in g/cm^3 + double density(const GenericTrack&) const; + /// \brief Returns the electron fraction + double ye(const GenericTrack&) const; + /// \brief Returns the radius of the Earth in km. + double GetRadius() const {return radius;} + /// \brief Returns the altitude of the top of the simulated atmosphere in km. + /// This is the altitude from which any initial neutrino flux will begin propagating. + double GetAtmosphereHeight() const {return atm_height;} + /// \brief Set the altitude of the top of the simulated atmosphere. + /// This function invalidates all tracks previously constructed for this body, so be sure to + /// use it before constructing any tracks, or to replace any previously constructed tracks. + /// \param height New height of the top of the atmosphere, in km. + void SetAtmosphereHeight(double height); + + /// \brief Construst a track with the given zenith angle, starting at the top of the atmosphere + /// \param phi Zenith angle with which the track will arrive at the surface of the Earth, + /// possibly after passing through the Earth + Track MakeTrack(double phi); + /// \brief Construst a track with the given zenith angle, starting at the top of the atmosphere + /// \param cosphi Cosine of the zenith angle with which the track will arrive at the surface of + /// the Earth, possibly after passing through the Earth + 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 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 five columns. + /// The first one must run 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 3rd column must correspond to the energy in GeV + /// in descending order over any relevant range + /// The 4th and 5th columns must correspond to the differential E_nu and Mu_nu + /// flux corresponding 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 + 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; + + // registration body + namespace detail{ + class registerBody{ + public: + registerBody(const std::string& name,std::function(hid_t)> fdeserialize); + }; + class registerTrack{ + public: + registerTrack(const std::string& name,std::function(hid_t)> fdeserialize); + }; + } + + std::function(hid_t)> GetBodyDeserializer(std::string); + std::function(hid_t)> GetTrackDeserializer(std::string); + + +// body registration +#define ASW(a, b) a ## b + +#define NUSQUIDS_REGISTER_BODY(classname) \ +namespace{ nusquids::detail::registerBody ASW(body_registerer,classname)(classname::GetName(),classname::Deserialize); nusquids::detail::registerTrack ASW(track_registerer,classname)(classname::Track::GetName(),classname::Track::Deserialize);} + +} // close namespace + +#endif diff --git a/nuSQuIDS.h b/nuSQuIDS.h new file mode 100644 index 00000000..74d08731 --- /dev/null +++ b/nuSQuIDS.h @@ -0,0 +1,1951 @@ + /****************************************************************************** + * 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 * + ******************************************************************************/ + + +#ifndef NUSQUIDS_nuSQUIDS_H +#define NUSQUIDS_nuSQUIDS_H + +#if __cplusplus < 201103L +#error C++11 compiler required. Update your compiler and use the flag -std=c++11 +#endif + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include +#include + +#include "nuSQuIDS/aligned_alloc.h" +#include "nuSQuIDS/body.h" +#include "nuSQuIDS/marray.h" +#include "nuSQuIDS/taudecay.h" +#include "nuSQuIDS/ThreadPool.h" +#include "nuSQuIDS/xsections.h" + +//#define FixCrossSections + +//#define UpdateInteractions_DEBUG + +namespace nusquids{ + +enum NeutrinoType { + neutrino=0b01, + antineutrino=0b10, + both=0b11 +}; + +enum Basis { + mass=0b01, + flavor=0b10, + interaction=0b11 +}; + + +///\class nuSQUIDS +///\brief nu-SQuIDS main class +class nuSQUIDS: public squids::SQuIDS { + // nuSQUIDSAtm is a friend class so it can use H0 to evolve operators + // and thus evaluate expectation values. + template + friend class nuSQUIDSAtm; +protected: + /// \brief Sets the basis in which the problem will be solved. + /// + /// If interaction basis is used the projectors will be evolved at + /// every time step. On the other hand, if mass basis is used no evolution + /// is performed. + Basis basis = interaction; + /// \brief number of neutrino flavors. + unsigned int numneu; + /// \brief number of energy nodes. + unsigned int ne = nx; + + ///Get the fraction of the of the number density of targets which is of each target type. + ///These fractions sum to one. + ///\return a vector of fractions, in the same order as the targets in int_struct + std::vector GetTargetNumberFractions() const; + + ///Get the number density of all nuclear targets at the current position + ///\return a vector of number densities, in the same order as the targets in int_struct + std::vector GetTargetNumberDensities() const; + + /// \brief Updates the interaction length arrays. + /// + /// Uses GetTargetNumberDensities() together with the stored cross section + /// information to update: nuSQUIDS#invlen_NC, nuSQUIDS#invlen_CC, and nuSQUIDS#invlen_INT. + void UpdateInteractions(); + + /// \brief Contains the energy nodes. + marray E_range; + + /// \brief Contains the spaces between nodes. + /// + /// It has length len(E_range)-1. + marray delE; + + /// \brief Interface that calculates and interpolates neutrino cross sections. + std::shared_ptr ncs; + + /// \brief Interface that calculate and interpolates tau decay spectral functions. + TauDecaySpectra tdc; + + /// \brief NT keeps track if the problem consists of neutrinos, antineutrinos, or both. + NeutrinoType NT = both; + + ///aligned arrays should be aligned to at least this many items + constexpr static size_t preferred_alignment=4; + constexpr static uint8_t log2(size_t n) { + return((n>1) ? log2(n/2)+1 : 0); + } + ///Round n up to the nearest multiple of preferred_alignment, doing nothing if n is already + ///such a multiple + constexpr static size_t round_up_to_aligned(size_t n){ + return((n&~(preferred_alignment-1))+((n&(preferred_alignment-1))?preferred_alignment:0)); + } + + public: + /// \brief Struct that contains all cross section information. + struct InteractionStructure { + /// \brief the configured target material types + /// \details These are recorded in the same order that they appear in dNdE_CC, + /// dNdE_NC, sigma_CC, and sigma_NC + std::vector targets; + /// \brief Neutrino charge current differential cross section with respect to + /// the outgoing lepton energy. + /// + /// This array is constructed when InitializeInteractionVectors() is called and + /// is initialized when InitializeInteractions() is called. Its dimensions are: + /// - target type (proton, neutron, etc.) + /// - neutrino type (neutrino, antineutrino) + /// - neutrino flavor (electron, muon, tau) + /// - initial energy + /// - final energy + marray> dNdE_CC; + /// \brief Neutrino neutral current differential cross section with respect to + /// the outgoing lepton energy. + /// + /// This array is constructed when InitializeInteractionVectors() is called and + /// is initialized when InitializeInteractions() is called. Its dimensions are: + /// - target type (proton, neutron, etc.) + /// - neutrino type (neutrino, antineutrino) + /// - neutrino flavor (electron, muon, tau) + /// - initial energy + /// - final energy + marray> dNdE_NC; + /// \brief Glashow resonance differential cross section (electron antineutrino only) + /// \details Dimensions are initial and final energy node + marray> dNdE_GR; + /// \brief Array that contains the neutrino charge current cross section. + /// \details The dimensions of this array are: + /// - target type (proton, neutron, etc.) + /// - neutrino type (neutrino, antineutrino) + /// - neutrino flavor (electron, muon, tau) + /// - energy + /// Its contents are in natural units, i.e. eV^-2. It is + /// initialized by InitializeInteractions() . + marray> sigma_CC; + /// \brief Array that contains the neutrino neutral current cross section. + /// \details The dimensions of this array are: + /// - target type (proton, neutron, etc.) + /// - neutrino type (neutrino, antineutrino) + /// - neutrino flavor (electron, muon, tau) + /// - energy + /// Its contents are in natural units, i.e. eV^-2. It is + /// initialized by InitializeInteractions() . + marray> sigma_NC; + /// \brief Glashow resonance cross section (electron antineutrino only) + /// \details 1 entry per energy node, in natural units + marray> sigma_GR; + /// \brief Array that contains the tau decay spectrum to all particles. + /// \details The first dimension corresponds to initial tau energy and the + /// second one to the outgoing lepton. + marray> dNdE_tau_all; + /// \brief Array that contains the tau decay spectrum to leptons. + /// \brief Array that contains the tau decay spectrum to all particles. + /// \details The first dimension corresponds to initial tau energy and the + /// second one to the outgoing lepton. + marray> dNdE_tau_lep; + /// \brief Default constructor + InteractionStructure(): + dNdE_CC(aligned_allocator{log2(preferred_alignment*sizeof(double))}), + dNdE_NC(aligned_allocator(log2(preferred_alignment*sizeof(double)))), + dNdE_GR(aligned_allocator(log2(preferred_alignment*sizeof(double)))), + sigma_CC(aligned_allocator(log2(preferred_alignment*sizeof(double)))), + sigma_NC(aligned_allocator(log2(preferred_alignment*sizeof(double)))), + sigma_GR(aligned_allocator(log2(preferred_alignment*sizeof(double)))), + dNdE_tau_all(aligned_allocator(log2(preferred_alignment*sizeof(double)))), + dNdE_tau_lep(aligned_allocator(log2(preferred_alignment*sizeof(double)))) + {} + /// \brief Assignment operator + InteractionStructure& operator=(const InteractionStructure& other){ + dNdE_CC=other.dNdE_CC; + dNdE_NC=other.dNdE_NC; + dNdE_GR=other.dNdE_GR; + sigma_CC=other.sigma_CC; + sigma_NC=other.sigma_NC; + sigma_GR=other.sigma_GR; + dNdE_tau_all=other.dNdE_tau_all; + dNdE_tau_lep=other.dNdE_tau_lep; + return(*this); + } + /// \brief Move assignment operator + InteractionStructure& operator=(InteractionStructure&& other){ + dNdE_CC=std::move(other.dNdE_CC); + dNdE_NC=std::move(other.dNdE_NC); + dNdE_GR=std::move(other.dNdE_GR); + sigma_CC=std::move(other.sigma_CC); + sigma_NC=std::move(other.sigma_NC); + sigma_GR=std::move(other.sigma_GR); + dNdE_tau_all=std::move(other.dNdE_tau_all); + dNdE_tau_lep=std::move(other.dNdE_tau_lep); + return(*this); + } + /// \brief Equality operator + bool operator==(const InteractionStructure& other) const{ + if ( + dNdE_CC.size() != other.dNdE_CC.size() or + dNdE_NC.size() != other.dNdE_NC.size() or + dNdE_GR.size() != other.dNdE_GR.size() or + sigma_CC.size() != other.sigma_CC.size() or + sigma_NC.size() != other.sigma_NC.size() or + sigma_GR.size() != other.sigma_GR.size() or + dNdE_tau_all.size() != other.dNdE_tau_all.size() or + dNdE_tau_lep.size() != other.dNdE_tau_lep.size() + ) + { + return false; + } + + if ( not std::equal(dNdE_CC.begin(),dNdE_CC.end(),other.dNdE_CC.begin()) ) + return false; + + if ( not std::equal(dNdE_NC.begin(),dNdE_NC.end(),other.dNdE_NC.begin()) ) + return false; + + if ( not std::equal(dNdE_GR.begin(),dNdE_GR.end(),other.dNdE_GR.begin()) ) + return false; + + if ( not std::equal(sigma_CC.begin(),sigma_CC.end(),other.sigma_CC.begin()) ) + return false; + + if ( not std::equal(sigma_NC.begin(),sigma_NC.end(),other.sigma_NC.begin()) ) + return false; + + if ( not std::equal(sigma_GR.begin(),sigma_GR.end(),other.sigma_GR.begin()) ) + return false; + + if ( not std::equal(dNdE_tau_all.begin(),dNdE_tau_all.end(),other.dNdE_tau_all.begin()) ) + return false; + + if ( not std::equal(dNdE_tau_lep.begin(),dNdE_tau_lep.end(),other.dNdE_tau_lep.begin()) ) + return false; + // all is good - lets roll + return true; + } + /// \brief Inequality operator + bool operator!=(const InteractionStructure& other) const { + return not (*this==other); + } + /// \brief Copy constructor operator + InteractionStructure(InteractionStructure& other): + dNdE_CC(other.dNdE_CC),dNdE_NC(other.dNdE_NC),dNdE_GR(other.dNdE_GR), + sigma_CC(other.sigma_CC),sigma_NC(other.sigma_NC),sigma_GR(other.sigma_GR), + dNdE_tau_all(other.dNdE_tau_all), + dNdE_tau_lep(other.dNdE_tau_lep) + {} + /// \brief Move constructor operator + InteractionStructure(InteractionStructure&& other): + dNdE_CC(std::move(other.dNdE_CC)),dNdE_NC(std::move(other.dNdE_NC)),dNdE_GR(std::move(other.dNdE_GR)), + sigma_CC(std::move(other.sigma_CC)),sigma_NC(std::move(other.sigma_NC)),sigma_GR(std::move(other.sigma_GR)), + dNdE_tau_all(std::move(other.dNdE_tau_all)), + dNdE_tau_lep(std::move(other.dNdE_tau_lep)) + {} + }; + + /// \brief Contains the interaction information at the current position. + struct InteractionState { + /// \brief Array that contains the inverse of the neutrino charge current mean free path. + /// \details The array contents are in natural units (i.e. eV) and is update when + /// UpdateInteractions() is called. The first dimension corresponds to the neutrino type, + /// the second to the flavor, and the last one to the energy. + marray> invlen_CC; + /// \brief Array that contains the inverse of the neutrino neutral current mean free path. + /// \details The array contents are in natural units (i.e. eV) and is update when + /// UpdateInteractions() is called. The first dimension corresponds to the neutrino type, + /// the second to the flavor, and the last one to the energy. + marray> invlen_NC; + /// \brief Glashow resonance inverse interaction length (electron antineutrino only) + /// \details 1 entry per energy node + marray> invlen_GR; + /// \brief Array that contains the inverse of the neutrino total mean free path. + /// \details The array contents are in natural units (i.e. eV) and is update when + /// UpdateInteractions() is called. Numerically it is just nuSQUIDS::invlen_NC and nuSQUIDS::invlen_CC + /// added together. + marray> invlen_INT; + + /// \brief Default constructor + InteractionState(): + invlen_CC(aligned_allocator(log2(preferred_alignment*sizeof(double)))), + invlen_NC(aligned_allocator(log2(preferred_alignment*sizeof(double)))), + invlen_GR(aligned_allocator(log2(preferred_alignment*sizeof(double)))), + invlen_INT(aligned_allocator(log2(preferred_alignment*sizeof(double)))) + {} + /// \brief Copy constructor + InteractionState(InteractionState& other): + invlen_CC(other.invlen_CC),invlen_NC(other.invlen_NC), + invlen_GR(other.invlen_GR),invlen_INT(other.invlen_INT) + {} + /// \brief Move constructor + InteractionState(InteractionState&& other): + invlen_CC(std::move(other.invlen_CC)),invlen_NC(std::move(other.invlen_NC)), + invlen_GR(std::move(other.invlen_GR)),invlen_INT(std::move(other.invlen_INT)) + {} + /// \brief Assignment operator + InteractionState& operator=(const InteractionState& other){ + invlen_NC=other.invlen_NC; + invlen_CC=other.invlen_CC; + invlen_GR=other.invlen_GR; + invlen_INT=other.invlen_INT; + return(*this); + } + /// \brief Move assignment operator + InteractionState& operator=(InteractionState&& other){ + invlen_NC=std::move(other.invlen_NC); + invlen_CC=std::move(other.invlen_CC); + invlen_GR=std::move(other.invlen_GR); + invlen_INT=std::move(other.invlen_INT); + return(*this); + } + /// \brief Equality operator + bool operator==(const InteractionState& other) const{ + if (invlen_CC.size() != other.invlen_CC.size() or + invlen_NC.size() != other.invlen_NC.size() or + invlen_GR.size() != other.invlen_GR.size() or + invlen_INT.size() != other.invlen_INT.size()) + return false; + + if ( not std::equal(invlen_CC.begin(),invlen_CC.end(),other.invlen_CC.begin()) ) + return false; + + if ( not std::equal(invlen_NC.begin(),invlen_NC.end(),other.invlen_NC.begin()) ) + return false; + + if ( not std::equal(invlen_GR.begin(),invlen_GR.end(),other.invlen_GR.begin()) ) + return false; + + if ( not std::equal(invlen_INT.begin(),invlen_INT.end(),other.invlen_INT.begin()) ) + return false; + + // all is good - lets roll + return true; + } + /// \brief Inequality operator + bool operator!=(const InteractionState& other) const { + return not (*this==other); + } + }; + + protected: + /// \brief Interaction struct that contains all cross section information tabulated. + std::shared_ptr int_struct; + /// \brief Current set of interaction information + InteractionState int_state; + + /// \brief Length upon which the neutrino fluxes will be positivized. + double positivization_scale; + + /// \brief Body where the neutrino propagation takes place. + std::shared_ptr body; + /// \brief Trayectory within the body. + /// \details Stores the position within the body and its updated every evolution + /// step. + std::shared_ptr track; + + /// \brief SU_vector that represents the neutrino square mass difference matrix in the mass basis. + /// It is used to construct nuSQUIDS#H0_array and H0() + squids::SU_vector DM2; + /// \brief Stores the time independent hamiltonian corresponding to each energy node. + marray H0_array; + + /// \brief Product of physical constants needed by HI. + double HI_constants; + /// \brief The density of the body at the current point on the track + double current_density; + /// \brief The electron fraction of the body at the current point on the track + double current_ye; + /// \brief The external flux from the body at the current point on the track + marray current_external_flux; + + /// \brief Mass basis projectors. + /// \details The i-entry corresponds to the projector in the ith mass eigenstate. + marray b0_proj; + /// \brief Flavor basis projectors. + /// \details The first dismension corresponds to the neutrino type. When NeutrinoType = both, then + /// the first dimension equal 0 means neutrinos and 1 antineutrinos. The second dimension corresponds + /// to the flavor eigenstates where 0 corresponds to nu_e, 1 to nu_mu, 2 to nu_tau, and the others + /// to sterile neutrinos. + marray b1_proj; + /// \brief Mass basis projectors in the interaction picture. + /// The index meaning are the same as nuSQUIDS#b1_proj but for mass eigenstates, + /// with an added third dimension that corresponds to the energy node index. + marray evol_b0_proj; + /// \brief Flavor basis projectors in the interaction picture. + /// The index meaning are the same as nuSQUIDS#b1_proj , + /// with an added third dimension that corresponds to the energy node index. + marray evol_b1_proj; + + ///Cache of precalculated results for InteractionsRho. + ///Used only when interactions are on and oscillations are off. + marray interaction_cache; + ///Backing storage for the SU_vectors in interaction_cache so that each can + ///be ensured to be correctly aligned for best performance. + std::unique_ptr interaction_cache_store; + ///Total size of interaction_cache_store + size_t interaction_cache_store_size; + + marray nc_factors; + marray tau_hadlep_decays; + marray tau_lep_decays; + marray gr_factors; + + ///Initialize the interaction cahce data structures to the corretc sizes for + ///the current problem. + ///Should be invoked only when interactions are on and oscillations are off. + void SetUpInteractionCache(); + + /// \brief Evolves the flavor projectors in the interaction basis to a time t. + /// \details It uses H0() to evolve SQUIDS#b0_proj and SQUIDS#b1_proj into + /// SQUIDS#evol_b0_proj and SQUIDS#evol_b1_proj. + /// \warning Since the RHS of the differential equation only involves flavor projectors + /// we do not current evolve mass projectors. + void EvolveProjectors(double t); + + // bool requirements + private: + /// \brief Boolean that signals that the full hamiltonian including extra potentials should be used for + /// the projector evolution + bool use_full_hamiltonian_for_projector_evolution = false; + /// \brief Boolean that toggles debug information to be printed + bool debug = false; + /// \brief Boolean that signals the object correct initialization. + bool inusquids = false; + /// \brief Boolean that signals that a Body object has being set. + bool ibody = false; + /// \brief Boolean that signals that the neutrino energies has being set. + bool ienergy = false; + /// \brief Boolean that signals that a Track object has being set. + bool itrack = false; + /// \brief Boolean that signals that an initial state has being set. + bool istate = false; + /// \brief Boolean that signals that interactions will be taken into account. + bool iinteraction = false; + /// \brief Whether interactions have been initialized + bool interactions_initialized = false; + /// \brief Boolean that signals that neutrino oscillations will be taken into account. + bool ioscillations = true; + /// \brief Boolean that signals that tau regeneration is being used. + bool tauregeneration = false; + /// \brief Whether the Glashsow resonance for electron antineutrinos is used. + bool iglashow = false; + /// \brief Boolean that signals that positivization will be enforced. + bool positivization = false; + /// \brief Boolean that allows to use fast constant density evolution technique. + bool allowConstantDensityOscillationOnlyEvolution = false; + /// \brief Boolean that signals that a progress bar will be printed. + bool progressbar = false; + /// \brief Integer to keep track of the progress bar evolution. + int progressbar_count = 0; + /// \brief Number of steps upon which the progress bar will be updated. + int progressbar_loop = 100; + /// \brief Time offset between SQuIDS time and Track(x). + double time_offset; + /// \brief Force flavor projections to be positive. + void PositivizeFlavors(); + /// \brief Set GSL differential cross section precision. + double gsl_int_precision = 1.e-3; + /// \brief Cutoff for low-pass filter applied during state density evolution + double evol_lowpass_cutoff = 0; + /// \brief Scale for low-pass filter applied during state density evolution + double evol_lowpass_scale = 0; + /// \brief If true the bodies can act as neutrino sources. + bool enable_neutrino_sources = false; + protected: + /// \brief Initializes flavor and mass projectors + /// \warning Antineutrinos are handle by means of the AntineutrinoCPFix() function + /// which temporary flips the CP phase. + void iniProjectors(); + /// \brief Reinitializes the flavor projectors. + /// \details It is called before evolving the system and every time the system + /// state is going to be set, since the definition of mass and flavor basis might + /// have changed. + void SetIniFlavorProyectors(); + /// \brief Initializes the time independent hamiltonian for each energy node. + /// \details It constructs nuSQUIDS#DM2 and nuSQUIDS#H0_array + void iniH0(); + /// \brief Changes the CP sign of the complex matrices stored in params. + void AntineutrinoCPFix(unsigned int irho); + /// \brief Serializes the initialization of the Body and Track objects. + /// @see ReadStateHDF5 + void SetBodyTrack(unsigned int,unsigned int,double*,unsigned int,double*); + + /// \brief General initilizer for the multi energy mode + /// @param E_vector Energy nodes [eV]. + /// @param xini The initial position of the system. + void init(marray E_vector,double xini=0.0); + + /// \brief Initilizer for the single energy mode + /// @param xini The initial position of the system. By default is set to 0. + void init(double xini=0.0); + + /// \brief Initilizes auxiliary cross section arrays. + /// \details The arrays are initialize, but not filled with contented. + /// To fill the arrays call InitializeInteractions(). + /// \param nTargets the number of target species for which to allocate space + void InitializeInteractionVectors(unsigned int nTargets); + + /// \brief Fills in auxiliary cross section arrays. + /// \details It uses nuSQUIDS#ncs and nuSQUIDS#tdc to fill in the values of + /// nuSQUIDS#dNdE_CC , nuSQUIDS#dNdE_NC , nuSQUIDS#sigma_CC , nuSQUIDS#sigma_NC , + /// nuSQUIDS#invlen_CC , nuSQUIDS#invlen_NC , nuSQUIDS#invlen_INT , + /// nuSQUIDS#dNdE_tau_all , nuSQUIDS#dNdE_tau_lep + /// in int_struct. + /// @see InitializeInteractionVectors + void GetCrossSections(); + private: + + /// \brief Prints progress bar. + /// \details To enable it call Set_ProgressBar() + void ProgressBar() const; + protected: + /// \brief Updates all quantities needed during the evaluation of the RHS. + /// @param x Position of the system. + /// \details Dependind on the basis used it evolves the flavor projectors + /// by means of EvolveProjectors(). Also, when interactions are considered, + /// it updates the interaction arrays by calling UpdateInteractions(). Finally, + /// it handles control to the user implemented updates via AddToPreDerive(). + /// @see AddToPreDerive + /// @see EvolveProjectors + /// @see UpdateInteractions + void PreDerive(double x); + /// \brief User supplied function that is called before performing a derivative. + /// @param x Position of the system. + /// @see PreDerive + virtual void AddToPreDerive(double x){} + + /// \brief Reads and constructs the object from an HDF5 file. + /// @param hdf5_filename Filename of the HDF5 to use for construction. + /// @param group Path to the group where the nuSQUIDS content will be saved. + /// @param iis InteractionStructure that contains valid cross sections. + /// \details By default all contents are assumed to be in the \c root of the HDF5 file + /// if the user wants a different location it can use \c group and \c save_cross_sections to + /// change the nuSQUIDS location and the cross section information location. + /// Finally, it calls AddToReadHDF5() to enable user defined parameters + /// to be read. + /// @see AddToReadHDF5 + /// @see WriteStateHDF5 + void ReadStateHDF5Internal(std::string hdf5_filename,std::string group = "/", std::shared_ptr iis = nullptr); + + public: + /************************************************************************************ + * CONSTRUCTORS + *************************************************************************************/ + + /// \brief Default void constructor. + nuSQUIDS(){} + + /// \brief Move constructor. + nuSQUIDS(nuSQUIDS&&); + + /// \brief Multiple energy mode constructor. + /// @param E_vector Energy nodes [eV]. + /// @param numneu Number of neutrino flavors. + /// @param NT NeutrinoType: neutrino,antineutrino, or both (simultaneous solution). + /// @param iinteraction Sets the neutrino noncoherent neutrino interactions on. + /// @param ncs Cross section object. + /// \details By default the energy scale is logarithmic and interactions are turn off. + /// \warning When interactions are present interpolation is performed to precalculate the neutrino + /// cross section which make take considertable time depending on the energy grid. + /// @see init + nuSQUIDS(marray E_vector,unsigned int numneu,NeutrinoType NT = both, + bool iinteraction = false, std::shared_ptr ncs = nullptr): + numneu(numneu),ncs(ncs),NT(NT),iinteraction(iinteraction) + {init(E_vector);} + + /// \brief Single energy mode constructor. + /// @param numneu Number of neutrino flavors. + /// @param NT NeutrinoType: neutrino or antineutrino. + /// \warning Interactions are not possible in the single energy mode, nor is simultaneous + /// neutrino-antineutrino solution (both) possible. + /// \details Constructors projectors and initializes Hamiltonian. + nuSQUIDS(unsigned int numneu, NeutrinoType NT = neutrino): + numneu(numneu),NT(NT),iinteraction(false) + {init();} + + /// \brief Constructor from a HDF5 filepath. + /// @param hdf5_filename Filename of the HDF5 to use for construction. + /// @param grp HDF5 file group path where the nuSQUIDS object is saved. + /// @param int_struct Interaction Structure to be used in object construction, + /// if this has already been obtained. May be left as NULL to cause + /// the interaction information to be read from the standard location. + /// \details Reads the HDF5 file and constructs the associated nuSQUIDS object + /// restoring all properties as well as the state. + /// @see ReadStateHDF5 + nuSQUIDS(std::string hdf5_filename, std::string grp = "/", + std::shared_ptr int_struct = nullptr): + int_struct(int_struct) + { + if(this->int_struct) + ReadStateHDF5Internal(hdf5_filename,grp,int_struct); + else + ReadStateHDF5(hdf5_filename, grp, ""); + } + + //*************************************************************** + virtual ~nuSQUIDS(); + + //*************************************************************** + ///\brief Move assigns a nuSQUIDS object from an existing object + nuSQUIDS& operator=(nuSQUIDS&&); + + protected: + /************************************************************************************ + * PHYSICS FUNCTIONS - SUPER IMPORTANT + *************************************************************************************/ + + /// \brief Returns the time independent part of the Hamiltonian + /// @param E Neutrino energy [eV] + /// @param irho Density matrix equation index. + /// \details \c irho is the index of the equation on which the Hamiltonian + /// is used. When not in NeutrinoType == \c both \c irho is only 0, but if + /// neutrinos and antineutrinos are solved simultaneously, then 0 is for + /// neutrinos and 1 for antineutrinos. + squids::SU_vector H0(double E, unsigned int irho) const; + + /// \brief Returns the time dependent part of the Hamiltonian at an energy node \c ie + /// @param ie Energy node + /// @param irho Density matrix equation index. + /// @see H0 for convention on \c irho. + virtual squids::SU_vector HI(unsigned int ie, unsigned int irho) const; + + /// \brief Returns noncoherent interaction term in the Hamiltonian at node \c ie + /// @param ie Energy node + /// @param irho Density matrix equation index. + /// @see H0 for convention on \c irho. + virtual squids::SU_vector GammaRho(unsigned int ie, unsigned int irho) const; + + /// \brief Returns the scalar quantity decay rate at an energy node \c ie + /// @param ie Energy node + /// @param iscalar scalar equation index. + /// \details When tau regenereation is considered \c iscalar = 0 corresponds + /// to the tau flux. + virtual double GammaScalar(unsigned int ie, unsigned int iscalar) const; + + /// \brief Returns interaction of the density matrix at an energy node \c ie + /// @param ie Energy node + /// @param irho Density matrix equation index. + /// @see H0 for convention on \c irho. + virtual squids::SU_vector InteractionsRho(unsigned int ie, unsigned int irho) const; + + /// \brief Returns scalar interactions at an energy node \c ie + /// @param ie Energy node + /// @param iscalar scalar equation index. + /// @see GammaScalar for convention on \c iscalar. + virtual double InteractionsScalar(unsigned int ie, unsigned int iscalar) const; + private: + /// \brief SQuIDS signature of HI + squids::SU_vector HI(unsigned int ie, unsigned int irho, double x) const {return HI(ie,irho);} + /// \brief SQuIDS signature of GammaRho + squids::SU_vector GammaRho(unsigned int ei, unsigned int irho, double x) const {return GammaRho(ei,irho);} + /// \brief SQuIDS signature of GammaScalar + double GammaScalar(unsigned int ei, unsigned int iscalar,double x) const {return GammaScalar(ei,iscalar);} + /// \brief SQuIDS signature of InteractionsRho + squids::SU_vector InteractionsRho(unsigned int ei, unsigned int irho, double x) const {return InteractionsRho(ei,irho);} + /// \brief SQuIDS signature of InteractionsScalar + double InteractionsScalar(unsigned int ei, unsigned int irho, double x) const {return InteractionsScalar(ei,irho);} + public: + /************************************************************************************ + * PUBLIC MEMBERS TO EVALUATE/SET/GET STUFF + *************************************************************************************/ + + /// \brief Sets the initial state in the single energy mode + /// @param ini_state Initial neutrino state. + /// @param basis Representation of the neutrino state either flavor or mass. + /// \details \c ini_state length has to be equal to \c numneu. If the basis is + /// flavor then the entries are interpret as nu_e, nu_mu, nu_tau, nu_sterile_1, ..., nu_sterile_n, + /// while if the mass basis is used then the first entries correspond to the active + /// mass eigenstates. + void Set_initial_state(const marray& ini_state, Basis basis = flavor); + + /// \brief Sets the initial state in the multiple energy mode when doing either neutrino or antineutrino only. + /// @param ini_state Initial neutrino state. + /// @param basis Representation of the neutrino state either flavor or mass. + /// \details \c ini_state first dimension has length equal to \c ne (number of energy nodes), while + /// the second dimension has length equal to \c numneu (number of flavors). If the basis is + /// flavor then the entries are interpret as nu_e, nu_mu, nu_tau, nu_sterile_1, ..., nu_sterile_n, + /// while if the mass basis is used then the first entries correspond to the active + /// mass eigenstates. + void Set_initial_state(const marray& ini_state, Basis basis = flavor); + + /// \brief Sets the initial state in the multiple energy mode when doing both neutrino and antineutrino. + /// @param ini_state Initial neutrino state. + /// @param basis Representation of the neutrino state either flavor or mass. + /// \details \c ini_state first dimension has length equal to \c ne (number of energy nodes). The second + /// dimension has to have size 2 and the first entry is for neutrinos while the second one for antineutrinos + /// Finally. the third dimension has length equal to \c numneu (number of flavors). If the basis is + /// flavor then the entries are interpret as nu_e, nu_mu, nu_tau, nu_sterile_1, ..., nu_sterile_n, + /// while if the mass basis is used then the first entries correspond to the active + /// mass eigenstates. + void Set_initial_state(const marray& ini_state, Basis basis = flavor); + + /// \brief Sets the body where the neutrino will propagate. + /// @param body Body object. + void Set_Body(std::shared_ptr body); + + /// \brief Sets neutrino trajectory + /// @param track Trajectory corresponding to the Body. + void Set_Track(std::shared_ptr track); + + /// \brief Sets the neutrino energy in single energy mode. + /// @param enu Neutrino energy [eV] + /// @pre NeutrinoType must be neutrino or antineutrino, not both. + void Set_E(double enu); + + /// \brief Evolves the system. + /// @pre Body, track, neutrino state, and energy must has being set + /// before calling this function. + void EvolveState(); + + /// \brief Returns the mass composition at a given node. + /// @param flv Neutrino flavor. + /// @param ie Energy node index. + /// @param rho Index of the equation, see details. + /// \details When NeutrinoType is \c both \c rho specifies wether one + /// is considering neutrinos (0) or antineutrinos (1). + + double EvalMassAtNode(unsigned int flv,unsigned int ie,unsigned int rho = 0) const; + /// \brief Returns the flavor composition at a given node. + /// @param flv Neutrino flavor. + /// @param ie Energy node index. + /// @param rho Index of the equation, see details. + /// \details When NeutrinoType is \c both \c rho specifies wether one + /// is considering neutrinos (0) or antineutrinos (1). + double EvalFlavorAtNode(unsigned int flv,unsigned int ie,unsigned int rho = 0) const; + + /// \brief Returns the mass composition at a given energy in the multiple energy mode. + /// averaging out the high frequencies. + /// @param flv Neutrino flavor. + /// @param enu Neutrino energy in natural units [eV]. + /// @param rho Index of the equation, see details. + double EvalMass(unsigned int flv,double enu,unsigned int rho) const; + + /// \brief Returns the mass composition at a given energy in the multiple energy mode. + /// averaging out the high frequencies. + /// @param flv Neutrino flavor. + /// @param enu Neutrino energy in natural units [eV]. + /// @param rho Index of the equation, see details. + /// @param scale scale upon which oscillations will be averaged out + /// @param avg bool array which is true for all scales that were averaged out + double EvalMass(unsigned int flv,double enu,unsigned int rho,double scale,std::vector& avr) const; + + /// \brief Returns the flavor composition at a given energy in the multiple energy mode. + /// @param flv Neutrino flavor. + /// @param enu Neutrino energy in natural units [eV]. + /// @param rho Index of the equation, see details. + /// \details When NeutrinoType is \c both \c rho specifies wether one + /// is considering neutrinos (0) or antineutrinos (1). + double EvalFlavor(unsigned int flv,double enu,unsigned int rho = 0) const; + + /// \brief Returns the flavor composition at a given energy in the multiple energy mode. + /// averaging out the high frequencies. + /// @param flv Neutrino flavor. + /// @param enu Neutrino energy in natural units [eV]. + /// @param rho Index of the equation, see details. + /// @param scale scale upon which oscillations will be averaged out + /// @param avg bool array which is true for all scales that were averaged out + double EvalFlavor(unsigned int flv,double enu,unsigned int rho,double scale,std::vector& avr) const; + + /// \brief Returns the mass composition in the single energy mode. + /// @param flv Neutrino flavor. + /// @param scale scale upon which oscillations will be averaged out + /// @param avg bool array which is true for all scales that were averaged out + double EvalMass(unsigned int flv, double scale, std::vector& avr) const; + + /// \brief Returns the flavor composition in the single energy mode. + /// @param flv Neutrino flavor. + /// @param scale scale upon which oscillations will be averaged out + /// @param avg bool array which is true for all scales that were averaged out + double EvalFlavor(unsigned int flv, double scale, std::vector& avr) const; + + /// \brief Returns the mass composition in the single energy mode. + /// @param flv Neutrino flavor. + double EvalMass(unsigned int flv) const; + + /// \brief Returns the flavor composition in the single energy mode. + /// @param flv Neutrino flavor. + double EvalFlavor(unsigned int flv) const; + + /// \brief Sets the cutoff for the state evolution low-pass filter. + /// @param val cutoff value + void Set_EvolLowPassCutoff(double val); + + /// \brief Sets the linear ramp size for the state evolution low-pass filter. + /// @param val Range in frequency space over which the linear ramp is applied. + void Set_EvolLowPassScale(double val); + + /// \brief Toggles tau regeneration on and off. + /// \param opt If \c true tau regeneration will be considered. + void Set_TauRegeneration(bool opt); + + /// \brief Toggles the effect of the Glashow resonance on and off. + /// \param opt If \c true resonant W^- production will be considered. + void Set_GlashowResonance(bool opt); + + /// \brief Toggles neutrino oscillations on and off. + /// \param opt If \c true neutrino oscillations will be considered. + void Set_IncludeOscillations(bool opt); + + /// \brief Toggles if non-hard processes constant density fast evolution will be used. + /// \param opt If \c true fast evolution will be attempted on constatnt density. + void Set_AllowConstantDensityOscillationOnlyEvolution(bool opt); + + /// \brief Toggles positivization of the flux. + /// @param opt If \c true the flux will be forced to be positive every \c positivization_step. + void Set_PositivityConstrain(bool opt); + + /// \brief Sets the positivization step. + /// @param step Sets the positivization step. + void Set_PositivityConstrainStep(double step); + + /// \brief Toggles the progress bar printing on and off + /// @param opt If \c true a progress bar will be printed. + void Set_ProgressBar(bool opt); + + /// \brief Returns the energy nodes values. + marray GetERange() const; + + /// \brief Returns number of energy nodes. + size_t GetNumE() const; + + /// \brief Returns the number of neutrino flavors. + unsigned int GetNumNeu() const; + + /// \brief Returns the number of rho equations. + unsigned int GetNumRho() const; + + /// \brief Returns true if noncoherent interactions are considered + bool GetUseInteractions() const {return iinteraction;} + + /// \brief Returns true if oscillations are considered + bool GetUseOscillations() const {return ioscillations;} + + /// \brinf Initiàlize interaction structure + void InitializeInteractions(); + + /// \brief Returns the interaction structure. + std::shared_ptr GetInteractionStructure() { + if(!interactions_initialized) + throw std::logic_error("Interactions not initialized"); + return int_struct; + } + + /// \brief Returns the interaction structure with constness. + std::shared_ptr GetInteractionStructure() const { + if(!interactions_initialized) + throw std::logic_error("Interactions not initialized"); + return int_struct; + } + + /// \brief Return the Hamiltonian at the current time + squids::SU_vector GetHamiltonian(unsigned int ei, unsigned int rho = 0); + /// \brief Returns the state + const squids::SU_vector& GetState(unsigned int ei,unsigned int rho = 0) const{ + return state[ei].rho[rho]; + } + /// \brief Returns the flavor projector + const squids::SU_vector& GetFlavorProj(unsigned int flv,unsigned int rho = 0) const{ + return b1_proj[rho][flv]; + } + /// \brief Returns the mass projector + const squids::SU_vector& GetMassProj(unsigned int flv,unsigned int rho = 0) const{ + return b0_proj[flv]; + } + + /// \brief Returns the trajectory object. + std::shared_ptr GetTrack(); + /// \brief Returns the body object. + std::shared_ptr GetBody(); + + const Track& GetTrackFast() const{ return(*track); } + + /// \brief Writes the object into an HDF5 file. + /// @param hdf5_filename Filename of the HDF5 to use for construction. + /// @param group Path to the group where the nuSQUIDS content will be saved. + /// @param save_cross_sections If \c true the cross section tables will also be saved. + /// @param cross_section_grp_loc Path to the group where the cross section will be saved. + /// @param overwrite If true will overwrite previous file, else it will append. + /// \details By default all contents are saved to the \c root of the HDF5 file + /// if the user wants a different location it can use \c group and \c save_cross_sections to + /// change the nuSQUIDS location and the cross section information location. Furthermore, the + /// \c save_cross_sections flag, which defaults to \c false decides if the cross section + /// will be saved. Finally, it calls AddToWriteHDF5() to enable user defined parameters + /// to be saved. + /// @see AddToWriteHDF5 + /// @see ReadStateHDF5 + void WriteStateHDF5(std::string hdf5_filename,std::string group = "/",bool save_cross_sections = true, std::string cross_section_grp_loc = "",bool overwrite = true) const; + + /// \brief User function to write user defined properties from an HDF5 file. + /// @param hdf5_loc_id HDF5 group id + /// \details This function will be called by WriteStateHDF5() which will provide + /// the hdf5_loc_id enabling the user to read user defined contents. + /// @see WriteStateHDF5 + /// @see AddToReadHDF5 + virtual void AddToWriteHDF5(hid_t hdf5_loc_id) const; + + /// \brief Reads and constructs the object from an HDF5 file. + /// @param hdf5_filename Filename of the HDF5 to use for construction. + /// @param group Path to the group where the nuSQUIDS content will be saved. + /// @param cross_section_grp_loc Path to the group where the cross section will be saved. + /// \details By default all contents are assumed to be in the \c root of the HDF5 file + /// if the user wants a different location it can use \c group and \c save_cross_sections to + /// change the nuSQUIDS location and the cross section information location. + /// Finally, it calls AddToReadHDF5() to enable user defined parameters + /// to be read. + /// @see AddToReadHDF5 + /// @see WriteStateHDF5 + void ReadStateHDF5(std::string hdf5_filename,std::string group = "/", std::string cross_section_grp_loc = ""); + + + /// \brief User function to read user defined properties from an HDF5 file. + /// @param hdf5_loc_id HDF5 group id + /// \details This function will be called by ReadStateHDF5() which will provide + /// the hdf5_loc_id enabling the user to read user defined contents. + /// @see ReadStateHDF5 + /// @see AddToWriteHDF5 + virtual void AddToReadHDF5(hid_t hdf5_loc_id); + + /// \brief Sets the mixing angle th_ij. + /// @param i the (zero-based) index of the first state + /// @param j the (zero-based) index of the second state + /// must be larger than \c i. + /// @param angle Angle to use in radians. + /// \details Sets the neutrino mixing angle. In our zero-based convention, e.g., the th_12 is i = 0, j = 1.,etc. + virtual void Set_MixingAngle(unsigned int i, unsigned int j,double angle); + + /// \brief Returns the mixing angle th_ij in radians. + /// @param i the (zero-based) index of the first state + /// @param j the (zero-based) index of the second state + /// must be larger than \c i. + /// \details Gets the neutrino mixing angle. In our zero-based convention, e.g., the th_12 is i = 0, j = 1.,etc. + double Get_MixingAngle(unsigned int i, unsigned int j) const; + + /// \brief Sets the CP phase for the ij-rotation. + /// @param i the (zero-based) index of the first state + /// @param j the (zero-based) index of the second state + /// must be larger than \c i. + /// @param angle Phase to use in radians. + /// \details Sets the CP phase for the ij-rotation. In our zero-based convention, e.g., the delta_13 = delta_CP is i = 0, j = 2.,etc. + virtual void Set_CPPhase(unsigned int i, unsigned int j,double angle); + + /// \brief Returns the CP phase of the ij-rotation in radians. + /// @param i the (zero-based) index of the first state + /// @param j the (zero-based) index of the second state + /// must be larger than \c i. + /// \details Gets the CP phase for the ij-rotation. In our zero-based convention, e.g., the delta_13 = delta_CP is i = 0, j = 2.,etc. + double Get_CPPhase(unsigned int i, unsigned int j) const; + + /// \brief Sets the square mass difference with respect to first mass eigenstate + /// @param i the (zero-based) index of the second state + /// must be larger than \c 0. + /// @param sq Square mass difference in eV^2. + /// \details Sets square mass difference with respect to the first mass eigenstate. In our zero-based convention, e.g., the \f$\Delta m^2_{12}\f$ corresponds to (i = 1),etc. + virtual void Set_SquareMassDifference(unsigned int i, double sq); + + /// \brief Returns the square mass difference between state-i and the first mass eigenstate. + /// @param i the (zero-based) index of the first state + /// must be larger than \c i. + /// \details Returns square mass difference with respect to the first mass eigenstate. In our zero-based convention, e.g., the \f$\Delta m^2_{12}\f$ corresponds to (i = 1),etc. + double Get_SquareMassDifference(unsigned int i) const; + + /// \brief Sets the mixing parameters to default. + void Set_MixingParametersToDefault(); + + /// \brief Sets the basis on which the evolution will be perfomed. + /// @param basis Evolution basis can be either \c mass or \c interaction. + /// \details By detault the solution is perfomed in the interaction basis, which + /// is recommended for most problems. One might consider the user of the mass + /// basis when collective effects, defined by new + /// interactions introduce phases that cannot be cancelled as one goes to the + /// interaction basis. + void Set_Basis(Basis basis); + + /// \brief Sets GSL differential cross section integrator precision + /// @param gsl_int_precision Double that specifies the differential cross section integrator precision. + void Set_GSL_DifferentialCrossSection_Integrator_Precision(double gsl_int_precision_){ + gsl_int_precision = gsl_int_precision_; + } + + /// \brief Returns the GSL differential cross section integrator precision + double Get_GSL_DifferentialCrossSection_Integrator_Precision() const { + return gsl_int_precision; + } + + /// \brief Sets the debug information on/off + /// @param debug_ Boolean that toggles debuging on and off. + void Set_Debug(bool debug_){ + debug = debug_; + } + + int Get_DebugCounter() const{ return progressbar_count; } + + /// \brief Returns lepton mixing matrix + std::unique_ptr GetTransformationMatrix() const { + return params.GetTransformationMatrix(GetNumNeu()); + } + + /// \brief Returns the neutrino interaction cross sections + std::shared_ptr GetNeutrinoCrossSections() { + return(ncs); + } + + /// \brief Returns the neutrino interaction cross sections + void SetNeutrinoCrossSections(std::shared_ptr xs) { + ncs=xs; + interactions_initialized=false; + } + + /// \brief Enables neutrino flux emission from bodies + void Set_NeutrinoSources(bool enable_neutrino_sources_){ + if(enable_neutrino_sources_){ + current_external_flux.resize(std::vector{ne,nrhos,numneu}); + std::fill(current_external_flux.begin(),current_external_flux.end(),0.0); + } + if(not iinteraction and enable_neutrino_sources_) + Set_OtherRhoTerms(true); + enable_neutrino_sources = enable_neutrino_sources_; + } + + /// \brief Returns true if neutrino flux emission from bodies is considered + bool Get_NeutrinoSources(){ + return enable_neutrino_sources; + } +}; + +/** + * The following class provides functionalities + * for atmospheric neutrino experiments + * where a collection of trajectories is explored. + */ + +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 = BaseNusType; + private: + /// \brief Random number generator + gsl_rng * r_gsl; + // /// \brief Mixing matrix + // std::unique_ptr U{nullptr,[](gsl_matrix_complex*m){delete m;}}; + /// \brief Internal units + const squids::Const units; + /// \brief Boolean that signals that a progress bar will be printed. + bool progressbar = false; + /// \brief Linear interpolator. + double LinInter(double x,double xM, double xP, double yM, double yP) const { + double f2=(x-xM)/(xP-xM); + double f1=1-f2; + return f1*yM + f2*yP; + } + /// \brief Boolean that signals that an initial state has being set. + bool iinistate; + /// \brief Boolean that signals the object correct initialization. + bool inusquidsatm; + /// \brief Contains the cos(zenith) nodes. + marray costh_array; + /// \brief Contains the energy nodes [eV]. + marray enu_array; + /// \brief Contains the log of energy nodes (in log of eV). + marray log_enu_array; + /// \brief Contains the nuSQUIDS objects for each zenith. + std::vector nusq_array; + + /// \brief Contains the Earth in atmospheric configuration. + 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 + std::shared_ptr ncs; + /// \brief Contains the interaction information structure. + std::shared_ptr int_struct; + /// \brief the number of threads to use when performing the evolution + unsigned int evalThreads; + public: + /************************************************************************************ + * CONSTRUCTORS + *************************************************************************************/ + + /// \brief Basic constructor. + /// @param costh_array One dimensional array containing zenith angles to be calculated. + /// @param energy_min Minimum neutrino energy value [ev]. + /// @param energy_max Maximum neutrino energy value [eV]. + /// @param energy_div Number of energy divisions. + /// @param numneu Number of neutrino flavors. + /// @param NT Signals the neutrino type : neutrino, antineutrion or both (simultaneous solution) + /// @param iinteraction Sets the neutrino noncoherent neutrino interactions on. + /// \details By defaults interactions are not considered and the neutrino energy scale is assume logarithmic. + template + nuSQUIDSAtm(marray costh_array, ArgTypes&&... args): + costh_array(costh_array), + evalThreads(1) + { + gsl_rng_env_setup(); + const gsl_rng_type * T_gsl = gsl_rng_default; + r_gsl = gsl_rng_alloc (T_gsl); + + 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]); + } + + 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); }); + for(unsigned int ie = 0 ; ie < enu_array.size(); ie++) + log_enu_array[ie] = log(enu_array[ie]); + + inusquidsatm = true; + } + + /// \brief Constructor from a HDF5 filepath. + /// @param hdf5_filename Filename of the HDF5 to use for construction. + /// \details Reads the HDF5 file and construct the associated nuSQUIDSAtm object + /// restoring all properties as well as the state. + /// @see ReadStateHDF5 + nuSQUIDSAtm(std::string hdf5_filename) { + gsl_rng_env_setup(); + const gsl_rng_type * T_gsl = gsl_rng_default; + r_gsl = gsl_rng_alloc (T_gsl); + ReadStateHDF5(hdf5_filename); + //U = nusq_array[0].GetTransformationMatrix(); + } + + /// \brief Move constructor. + nuSQUIDSAtm(nuSQUIDSAtm&& other): + progressbar(other.progressbar), + iinistate(other.iinistate), + inusquidsatm(other.inusquidsatm), + costh_array(std::move(other.costh_array)), + enu_array(std::move(other.enu_array)), + log_enu_array(std::move(other.log_enu_array)), + nusq_array(std::move(other.nusq_array)), + earth_atm(std::move(other.earth_atm)), + track_array(std::move(other.track_array)), + ncs(std::move(other.ncs)), + int_struct(std::move(other.int_struct)), + evalThreads(other.evalThreads) + { + other.inusquidsatm = false; + } + + //*************************************************************** + ///\brief Move assigns a nuSQUIDSAtm object from an existing object + nuSQUIDSAtm& operator=(nuSQUIDSAtm&& other){ + if(&other==this) + return(*this); + + progressbar = other.progressbar; + iinistate = other.iinistate; + inusquidsatm = other.inusquidsatm; + costh_array = std::move(other.costh_array); + enu_array = std::move(other.enu_array); + log_enu_array = std::move(other.log_enu_array); + nusq_array = std::move(other.nusq_array); + earth_atm = std::move(other.earth_atm); + track_array = std::move(other.track_array); + ncs = std::move(other.ncs); + int_struct = std::move(other.int_struct); + evalThreads = other.evalThreads; + + // initial nusquids object render useless + other.inusquidsatm = false; + + return(*this); + } + /************************************************************************************ + * 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. + /// @param basis Representation of the neutrino state either flavor or mass. + /// \details \c ini_state first dimension lenght has to be equal to the number of + /// zenith bins, the second dimension is the number of energy bins, and the third + /// dimension corresponds to the number of neutrino flavors. If the basis is + /// flavor then the entries are interpret as nu_e, nu_mu, nu_tau, nu_sterile_1, ..., nu_sterile_n, + /// while if the mass basis is used then the first entries correspond to the active + /// mass eigenstates. + void Set_initial_state(const marray& ini_flux, Basis basis=flavor){ + if(ini_flux.extent(0) != costh_array.extent(0)) + throw std::runtime_error("nuSQUIDSAtm::Error::First dimension of input array is incorrect."); + if(ini_flux.extent(1) != enu_array.extent(0)) + throw std::runtime_error("nuSQUIDSAtm::Error::Second dimension of input array is incorrect."); + unsigned int i = 0; + for(BaseSQUIDS& nsq : nusq_array){ + marray slice{ini_flux.extent(1),ini_flux.extent(2)}; + for(size_t j=0; j& ini_flux, Basis basis=flavor){ + if(ini_flux.extent(0) != costh_array.extent(0)) + throw std::runtime_error( + "nuSQUIDSAtm::Error::First dimension of input array is incorrect."); + unsigned int i = 0; + for(BaseSQUIDS& nsq : nusq_array){ + marray slice{ini_flux.extent(1),ini_flux.extent(2),ini_flux.extent(3)}; + for(size_t j=0; jSet_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(){ + if(not iinistate) + throw std::runtime_error("nuSQUIDSAtm::Error::State not initialized."); + if(not inusquidsatm) + throw std::runtime_error("nuSQUIDSAtm::Error::nuSQUIDSAtm not initialized."); + + if(nusq_array.front().GetUseInteractions()){ + // setting the interaction structure also on the main object + nusq_array.front().InitializeInteractions(); + int_struct = nusq_array.front().GetInteractionStructure(); + for(BaseSQUIDS& nsq : nusq_array){ + nsq.int_struct=int_struct; //copy identical cross section data + nsq.InitializeInteractions(); //still need to initialize intermediate buffers + } + } + + if(evalThreads==1){ + unsigned int i = 0; + for(BaseSQUIDS& nsq : nusq_array){ + if(progressbar){ + std::cout << "Calculating cos(th) = " + std::to_string(costh_array[i]) << std::endl; + i++; + } + nsq.EvolveState(); + if(progressbar){ + std::cout << std::endl; + } + } + } + else{ + bool overrideProgressBar=progressbar; + Set_ProgressBar(false); + ThreadPool tpool(evalThreads); + std::vector> tasks; + tasks.reserve(nusq_array.size()); + for(BaseSQUIDS& nsq : nusq_array) + tasks.emplace_back(tpool.enqueue([&](){ nsq.EvolveState(); })); + unsigned int i = 0; + for(const auto& task : tasks){ + task.wait(); + if(overrideProgressBar) + std::cout << "cos(th) = " + std::to_string(costh_array[i++]) << " complete" << std::endl; + } + Set_ProgressBar(overrideProgressBar); + } + } + + /// \brief Returns the flavor composition at a given energy and zenith. + /// @param flv Neutrino flavor. + /// @param costh Cosine of the zenith. + /// @param enu Neutrino energy [eV]. + /// @param rho Index of the equation, see details. + /// \details When NeutrinoType is \c both \c rho specifies wether one + /// is considering neutrinos (0) or antineutrinos (1). Bilinear interpolation + /// is done in the logarithm of the energy and cos(zenith). + double EvalFlavor(unsigned int flv,double costh,double enu,unsigned int rho = 0, bool randomize_production_height = false) const { + // here the energy enters in eV + if(not iinistate) + throw std::runtime_error("nuSQUIDSAtm::Error::State not initialized."); + if(not inusquidsatm) + throw std::runtime_error("nuSQUIDSAtm::Error::nuSQUIDSAtm not initialized."); + if(not ( rho == 0 or rho == 1)) + throw std::runtime_error("nuSQUIDSAtm::Error::EvalFlavor rho has to be 0 or 1."); + + if( costh < *costh_array.begin() or costh > *costh_array.rbegin()) + throw std::runtime_error("nuSQUIDSAtm::Error::EvalFlavor::cos(th) out of bounds."); + if( enu < *enu_array.begin() or enu > *enu_array.rbegin() ) + throw std::runtime_error("nuSQUIDSAtm::Error::EvalFlavor::neutrino energy out of bounds.(Emin = " + + std::to_string(*enu_array.begin()) + + ",Emax = " + + std::to_string(*enu_array.rbegin()) + + ", Enu = " + std::to_string(enu) + ")"); + + auto cthit=std::lower_bound(costh_array.begin(),costh_array.end(),costh); + if(cthit==costh_array.end()) + throw std::runtime_error("SQUIDS::GetExpectationValueD : x value not in the array."); + if(cthit!=costh_array.begin()) + cthit--; + size_t cth_M=std::distance(costh_array.begin(),cthit); + + /*double logE = log(enu); + auto logeit=std::lower_bound(log_enu_array.begin(),log_enu_array.end(),logE); + if(logeit==log_enu_array.end()) + throw std::runtime_error("SQUIDS::GetExpectationValueD : x value not in the array."); + if(logeit!=log_enu_array.begin()) + logeit--; + size_t loge_M=std::distance(log_enu_array.begin(),logeit);*/ + + auto eit=std::lower_bound(enu_array.begin(),enu_array.end(),enu); + if(eit==enu_array.end()) + throw std::runtime_error("SQUIDS::GetExpectationValueD : x value not in the array."); + if(eit!=enu_array.begin()) + eit--; + size_t loge_M=std::distance(enu_array.begin(),eit); + + EarthAtm::Track track=earth_atm->MakeTrackWithCosine(costh); + double delta_t_final = track.GetFinalX()-track.GetInitialX(); + if (randomize_production_height){ + double production_height = gsl_ran_flat(r_gsl,-15*units.km,15*units.km); + delta_t_final += production_height; + } + + // assuming offsets are zero + double delta_t_1 = nusq_array[cth_M].Get_t() - nusq_array[cth_M].Get_t_initial(); + double delta_t_2 = nusq_array[cth_M+1].Get_t() - nusq_array[cth_M+1].Get_t_initial(); + double delta_t_final_1 = nusq_array[cth_M].GetTrackFast().GetFinalX() - nusq_array[cth_M].GetTrackFast().GetInitialX(); + double delta_t_final_2 = nusq_array[cth_M+1].GetTrackFast().GetFinalX() - nusq_array[cth_M+1].GetTrackFast().GetInitialX(); + double t_inter = 0.5*(delta_t_final*delta_t_1/delta_t_final_1 + delta_t_final*delta_t_2/delta_t_final_2); + + squids::SU_vector H0_at_enu = nusq_array[0].H0(enu,rho); + + struct storage_type{ + squids::SU_vector evol_proj, temp1; + storage_type(unsigned int dim):evol_proj(dim),temp1(dim){} + }; + #ifdef SQUIDS_THREAD_LOCAL + static SQUIDS_THREAD_LOCAL + #endif + storage_type storage(H0_at_enu.Dim()); + + storage.evol_proj = nusq_array[0].GetFlavorProj(flv,rho).Evolve(H0_at_enu,t_inter); + + //coefficients for energy interpolation + double f2=(enu-enu_array[loge_M])/(enu_array[loge_M+1]-enu_array[loge_M]); + double f1=1-f2; + + storage.temp1 =f1*nusq_array[cth_M ].GetState(loge_M ,rho); + storage.temp1+=f2*nusq_array[cth_M ].GetState(loge_M+1,rho); + double phiM=squids::SUTrace(storage.temp1,storage.evol_proj); + + storage.temp1 =f1*nusq_array[cth_M+1].GetState(loge_M ,rho); + storage.temp1+=f2*nusq_array[cth_M+1].GetState(loge_M+1,rho); + double phiP=squids::SUTrace(storage.temp1,storage.evol_proj); + + //perform angular interpolation + return(LinInter(costh,costh_array[cth_M],costh_array[cth_M+1],phiM,phiP)); + } + + /// \brief Returns the flavor composition at a given energy and zenith. + /// @param flv Neutrino flavor. + /// @param costh Cosine of the zenith. + /// @param enu Neutrino energy [eV]. + /// @param rho Index of the equation, see details. + /// @param scale Scale is a float that sets the averaged oscillation scaled. + /// @param avr Is a bool vector of size of the number of neutrinos that signal which scaled has been averaged. + /// \details When NeutrinoType is \c both \c rho specifies wether one + /// is considering neutrinos (0) or antineutrinos (1). Bilinear interpolation + /// is done in the logarithm of the energy and cos(zenith). + double EvalFlavor(unsigned int flv,double costh,double enu,unsigned int rho,double scale,std::vector avr) const { + // here the energy enters in eV + if(not iinistate) + throw std::runtime_error("nuSQUIDSAtm::Error::State not initialized."); + if(not inusquidsatm) + throw std::runtime_error("nuSQUIDSAtm::Error::nuSQUIDSAtm not initialized."); + if(not ( rho == 0 or rho == 1)) + throw std::runtime_error("nuSQUIDSAtm::Error::EvalFlavor rho has to be 0 or 1."); + + if( costh < *costh_array.begin() or costh > *costh_array.rbegin()) + throw std::runtime_error("nuSQUIDSAtm::Error::EvalFlavor::cos(th) out of bounds."); + if( enu < *enu_array.begin() or enu > *enu_array.rbegin() ) + throw std::runtime_error("nuSQUIDSAtm::Error::EvalFlavor::neutrino energy out of bounds.(Emin = " + + std::to_string(*enu_array.begin()) + + ",Emax = " + + std::to_string(*enu_array.rbegin()) + + ", Enu = " + std::to_string(enu) + ")"); + + auto cthit=std::lower_bound(costh_array.begin(),costh_array.end(),costh); + if(cthit==costh_array.end()) + throw std::runtime_error("SQUIDS::GetExpectationValueD : x value not in the array."); + if(cthit!=costh_array.begin()) + cthit--; + size_t cth_M=std::distance(costh_array.begin(),cthit); + + /*double logE = log(enu); + auto logeit=std::lower_bound(log_enu_array.begin(),log_enu_array.end(),logE); + if(logeit==log_enu_array.end()) + throw std::runtime_error("SQUIDS::GetExpectationValueD : x value not in the array."); + if(logeit!=log_enu_array.begin()) + logeit--; + size_t loge_M=std::distance(log_enu_array.begin(),logeit);*/ + + auto eit=std::lower_bound(enu_array.begin(),enu_array.end(),enu); + if(eit==enu_array.end()) + throw std::runtime_error("SQUIDS::GetExpectationValueD : x value not in the array."); + if(eit!=enu_array.begin()) + eit--; + size_t loge_M=std::distance(enu_array.begin(),eit); + + EarthAtm::Track track=earth_atm->MakeTrackWithCosine(costh); + double delta_t_final = track.GetFinalX()-track.GetInitialX(); + + // assuming offsets are zero + double delta_t_1 = nusq_array[cth_M].Get_t() - nusq_array[cth_M].Get_t_initial(); + double delta_t_2 = nusq_array[cth_M+1].Get_t() - nusq_array[cth_M+1].Get_t_initial(); + double delta_t_final_1 = nusq_array[cth_M].GetTrackFast().GetFinalX() - nusq_array[cth_M].GetTrackFast().GetInitialX(); + double delta_t_final_2 = nusq_array[cth_M+1].GetTrackFast().GetFinalX() - nusq_array[cth_M+1].GetTrackFast().GetInitialX(); + double t_inter = 0.5*(delta_t_final*delta_t_1/delta_t_final_1 + delta_t_final*delta_t_2/delta_t_final_2); + + squids::SU_vector H0_at_enu = nusq_array[0].H0(enu,rho); + + struct storage_type{ + squids::SU_vector evol_proj, temp1; + storage_type(unsigned int dim):evol_proj(dim),temp1(dim){} + }; + #ifdef SQUIDS_THREAD_LOCAL + static SQUIDS_THREAD_LOCAL + #endif + storage_type storage(H0_at_enu.Dim()); + // preevolution buffer + std::unique_ptr evol_buffer(new double[H0_at_enu.GetEvolveBufferSize()]); + H0_at_enu.PrepareEvolve(evol_buffer.get(),t_inter,scale,avr); + storage.evol_proj = nusq_array[0].GetFlavorProj(flv,rho).Evolve(evol_buffer.get()); + + //coefficients for energy interpolation + double f2=(enu-enu_array[loge_M])/(enu_array[loge_M+1]-enu_array[loge_M]); + double f1=1-f2; + + storage.temp1 =f1*nusq_array[cth_M ].GetState(loge_M ,rho); + storage.temp1+=f2*nusq_array[cth_M ].GetState(loge_M+1,rho); + + double phiM=squids::SUTrace(storage.temp1,storage.evol_proj); + + storage.temp1 =f1*nusq_array[cth_M+1].GetState(loge_M ,rho); + storage.temp1+=f2*nusq_array[cth_M+1].GetState(loge_M+1,rho); + double phiP=squids::SUTrace(storage.temp1,storage.evol_proj); + + //perform angular interpolation + return(LinInter(costh,costh_array[cth_M],costh_array[cth_M+1],phiM,phiP)); + } + + /// \brief Writes the object into an HDF5 file. + /// @param hdf5_filename Filename of the HDF5 into which save the object. + /// \details All contents are saved to the \c root of the HDF5 file. + /// @see ReadStateHDF5 + void WriteStateHDF5(std::string filename, bool overwrite = true) const{ + if(not iinistate) + throw std::runtime_error("nuSQUIDSAtm::Error::State not initialized."); + if(not inusquidsatm) + throw std::runtime_error("nuSQUIDSAtm::Error::nuSQUIDSAtm not initialized."); + + hid_t file_id,root_id; + hid_t dset_id; + // create HDF5 file + if(overwrite) + file_id = H5Fcreate (filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + else + file_id = H5Fcreate (filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT, H5P_DEFAULT); + if (file_id < 0) + throw std::runtime_error("nuSQUIDS::Error::Cannot create file at " + filename + "."); + root_id = H5Gopen(file_id, "/",H5P_DEFAULT); + + // write the zenith range + hsize_t costhdims[1]={costh_array.extent(0)}; + dset_id = H5LTmake_dataset(root_id,"zenith_angles",1,costhdims,H5T_NATIVE_DOUBLE,costh_array.get_data()); + hsize_t energydims[1]={enu_array.extent(0)}; + dset_id = H5LTmake_dataset(root_id,"energy_range",1,energydims,H5T_NATIVE_DOUBLE,enu_array.get_data()); + + H5Gclose (root_id); + H5Fclose (file_id); + + unsigned int i = 0; + for(const BaseSQUIDS& nsq : nusq_array){ + // use only the first one to write the cross sections on /crosssections + // the last false is because the HDF5 file has already been created, should not overwrite + nsq.WriteStateHDF5(filename,"/costh_"+std::to_string(costh_array[i]),i==0,"crosssections",false); + i++; + } + } + + /// \brief Reads the object from an HDF5 file. + /// @param hdf5_filename Filename of the HDF5 to use for construction. + /// \details All contents are assumed to be saved to the \c root of the HDF5 file. + /// @see WriteStateHDF5 + void ReadStateHDF5(std::string hdf5_filename){ + hid_t file_id,group_id,root_id; + // open HDF5 file + H5File file(H5Fopen(hdf5_filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)); + if (file < 0) + throw std::runtime_error("nuSQUIDSAtm::ReadStateHDF5: Unable to open file: " + hdf5_filename + ". No such file or directory"); + root_id = H5Gopen(file, "/",H5P_DEFAULT); + group_id = root_id; + + // read the zenith range dimension + hsize_t costhdims[1]; + H5LTget_dataset_info(group_id, "zenith_angles", costhdims, nullptr, nullptr); + + double data[costhdims[0]]; + H5LTread_dataset_double(group_id, "zenith_angles", data); + costh_array.resize(std::vector {static_cast(costhdims[0])}); + for (unsigned int i = 0; i < costhdims[0]; i ++) + costh_array[i] = data[i]; + + hsize_t energydims[1]; + H5LTget_dataset_info(group_id, "energy_range", energydims, nullptr, nullptr); + + double enu_data[energydims[0]]; + H5LTread_dataset_double(group_id, "energy_range", enu_data); + enu_array.resize(std::vector{static_cast(energydims[0])});log_enu_array.resize(std::vector{static_cast(energydims[0])}); + for (unsigned int i = 0; i < energydims[0]; i ++){ + enu_array[i] = enu_data[i]; + log_enu_array[i] = log(enu_data[i]); + } + + H5Gclose(root_id); + + // resize apropiately the nuSQUIDSAtm container vector + nusq_array.clear(); + nusq_array = std::vector(costhdims[0]); + + unsigned int i = 0; + for(BaseSQUIDS& nsq : nusq_array){ + if(i==0){ + // read the cross sections stored in /crosssections + nsq.ReadStateHDF5(hdf5_filename,"/costh_"+std::to_string(costh_array[i]),"crosssections"); + if(nsq.iinteraction) + int_struct = nsq.GetInteractionStructure(); + } else { + // re-use the shared cross sections + nsq.ReadStateHDF5Internal(hdf5_filename,"/costh_"+std::to_string(costh_array[i]),int_struct); + } + i++; + } + earth_atm = std::dynamic_pointer_cast(nusq_array.front().GetBody()); + assert(earth_atm); + + iinistate = true; + inusquidsatm = true; + } + + /// \brief Sets the mixing parameters to default. + void Set_MixingParametersToDefault(){ + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_MixingParametersToDefault(); + } + } + + /// \brief Sets the mixing angle th_ij. + /// @param i the (zero-based) index of the first state + /// @param j the (zero-based) index of the second state + /// must be larger than \c i. + /// @param angle Angle to use in radians. + /// \details Sets the neutrino mixing angle. In our zero-based convention, e.g., the th_12 is i = 0, j = 1.,etc. + void Set_MixingAngle(unsigned int i, unsigned int j,double angle){ + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_MixingAngle(i,j,angle); + } + } + + /// \brief Returns the mixing angle th_ij in radians. + /// @param i the (zero-based) index of the first state + /// @param j the (zero-based) index of the second state + /// must be larger than \c i. + /// \details Gets the neutrino mixing angle. In our zero-based convention, e.g., the th_12 is i = 0, j = 1.,etc. + double Get_MixingAngle(unsigned int i, unsigned int j) const{ + return nusq_array[0].Get_MixingAngle(i,j); + } + + /// \brief Sets the CP phase for the ij-rotation. + /// @param i the (zero-based) index of the first state + /// @param j the (zero-based) index of the second state + /// must be larger than \c i. + /// @param angle Phase to use in radians. + /// \details Sets the CP phase for the ij-rotation. In our zero-based convention, e.g., the delta_13 = delta_CP is i = 0, j = 2.,etc. + void Set_CPPhase(unsigned int i, unsigned int j,double angle){ + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_CPPhase(i,j,angle); + } + } + + /// \brief Returns the CP phase of the ij-rotation in radians. + /// @param i the (zero-based) index of the first state + /// @param j the (zero-based) index of the second state + /// must be larger than \c i. + /// \details Gets the CP phase for the ij-rotation. In our zero-based convention, e.g., the delta_13 = delta_CP is i = 0, j = 2.,etc. + double Get_CPPhase(unsigned int i, unsigned int j) const{ + return nusq_array[0].Get_CPPhase(i,j); + } + + /// \brief Sets the square mass difference with respect to first mass eigenstate + /// @param i the (zero-based) index of the second state + /// must be larger than \c 0. + /// @param sq Square mass difference in eV^2. + /// \details Sets square mass difference with respect to the first mass eigenstate. In our zero-based convention, e.g., the \f$\Delta m^2_{12}\f$ corresponds to (i = 1),etc. + void Set_SquareMassDifference(unsigned int i,double sq){ + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_SquareMassDifference(i,sq); + } + } + + /// \brief Returns the square mass difference between state-i and the first mass eigenstate. + /// @param i the (zero-based) index of the first state + /// must be larger than \c i. + /// \details Returns square mass difference with respect to the first mass eigenstate. In our zero-based convention, e.g., the \f$\Delta m^2_{12}\f$ corresponds to (i = 1),etc. + double Get_SquareMassDifference(unsigned int i) const{ + return nusq_array[0].Get_SquareMassDifference(i); + } + + /// \brief Set the initial integration step size + /// \param h initial step size + void Set_h(double h){ + for(BaseSQUIDS& nsq : nusq_array) + nsq.Set_h(h); + } + + /// \brief Set the initial integration step size for a single energy node + /// \param h initial step size + /// \param idx energy node index + void Set_h(double h, unsigned int idx){ + nusq_array[idx].Set_h(h); + } + + /// \brief Set the maximum integration step size + /// \param h maximum step size + void Set_h_max(double h){ + for(BaseSQUIDS& nsq : nusq_array) + nsq.Set_h_max(h); + } + + /// \brief Set the maximum integration step size for a single energy node + /// \param h maximum step size + /// \param idx energy node index + void Set_h_max(double h, unsigned int idx){ + nusq_array[idx].Set_h_max(h); + } + + /// \brief Set the minimum integration step size + /// \param h minimum step size + void Set_h_min(double h){ + for(BaseSQUIDS& nsq : nusq_array) + nsq.Set_h_min(h); + } + + /// \brief Set the minimum integration step size for a single energy node + /// \param h minimum step size + /// \param idx energy node index + void Set_h_min(double h, unsigned int idx){ + nusq_array[idx].Set_h_min(h); + } + + /// \brief Set the allowed absolute numerical error + /// \param eps maximum allowed error + void Set_abs_error(double eps){ + for(BaseSQUIDS& nsq : nusq_array) + nsq.Set_abs_error(eps); + } + + /// \brief Set the allowed absolute numerical error for a single energy node + /// \param eps maximum allowed error + /// \param idx energy node index + void Set_abs_error(double eps, unsigned int idx){ + nusq_array[idx].Set_abs_error(eps); + } + + /// \brief Set the allowed relative numerical error + /// \param eps maximum allowed error + void Set_rel_error(double eps){ + for(BaseSQUIDS& nsq : nusq_array) + nsq.Set_rel_error(eps); + } + + /// \brief Set the allowed relative numerical error for a single energy node + /// \param eps maximum allowed error + /// \param idx energy node index + void Set_rel_error(double eps, unsigned int idx){ + nusq_array[idx].Set_rel_error(eps); + } + + /// \brief Sets the GSL solver + /// @param opt GSL stepper function. + void Set_GSL_step(gsl_odeiv2_step_type const * opt){ + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_GSL_step(opt); + } + } + + /// \brief Toggles the progress bar printing on and off + /// @param opt If \c true a progress bar will be printed. + void Set_ProgressBar(bool opt){ + progressbar = opt; + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_ProgressBar(opt); + } + } + + /// \brief Returns the interaction structure. + std::shared_ptr GetInteractionStructure() const { + return int_struct; + } + + /// \brief Returns number of energy nodes. + size_t GetNumE() const{ + return enu_array.extent(0); + } + + /// \brief Returns number of zenith nodes. + size_t GetNumCos() const{ + return costh_array.extent(0); + } + + /// \brief Returns the number of neutrino flavors. + unsigned int GetNumNeu() const{ + return nusq_array[0].GetNumNeu(); + } + + /// \brief Returns the number of rho equations. + unsigned int GetNumRho() const{ + return nusq_array[0].GetNumRho(); + } + + /// \brief Returns number of nodes. + size_t GetNumNodes() const{ + return nusq_array.size(); + } + + /// \brief Returns the (evolved) states of all nodes + /// @param rho Index of equation. In `both` propagation mode, neutrinos = 0 and + /// antineutrinos = 1. + marray GetStates(unsigned int rho = 0){ + marray states {GetNumNodes(),GetNumNeu()*GetNumNeu()}; + for(int i=0; i state_vec = nusq_array[i].GetState(0, rho).GetComponents(); + for(int j=0; j GetERange() const{ + return enu_array; + } + + /// \brief Returns the cos(zenith) nodes values. + marray GetCosthRange() const{ + return costh_array; + } + + /// \brief Returns the nuSQUIDSBase object for ci-th zenith. + BaseSQUIDS& GetnuSQuIDS(unsigned int ci) { + if(ci >= nusq_array.size()) + throw std::runtime_error("nuSQUIDSAtm::GetnuSQuIDS: input index out of range."); + return nusq_array[ci]; + } + + /// \brief Returns the vector of nuSQUIDSBase objects. + std::vector& GetnuSQuIDS() { + return nusq_array; + } + + /// \brief Toggles tau regeneration on and off. + /// @param opt If \c true tau regeneration will be considered. + void Set_TauRegeneration(bool opt){ + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_TauRegeneration(opt); + } + } + + void Set_GlashowResonance(bool opt){ + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_GlashowResonance(opt); + } + } + + /// \brief Toggles neutrino oscillations on and off. + /// @param opt If \c true neutrino oscillations will be considered. + void Set_IncludeOscillations(bool opt){ + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_IncludeOscillations(opt); + } + } + + /// \brief Toggles if non-hard processes constant density fast evolution will be used. + /// \param opt If \c true fast evolution will be attempted on constatnt density. + void Set_AllowConstantDensityOscillationOnlyEvolution(bool opt){ + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_AllowConstantDensityOscillationOnlyEvolution(opt); + } + } + + + /// \brief Toggles positivization of the flux. + /// @param opt If \c true the flux will be forced to be positive every \c positivization_step. + void Set_PositivityConstrain(bool opt){ + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_PositivityConstrain(opt); + } + } + /// \brief Stes the step upon which the positivity correction would be apply. + /// @param step The step upon which the positivization will take place. + void Set_PositivityConstrainStep(double step){ + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_PositivityConstrainStep(step); + } + } + + /// \brief Set the number of threads used for evolution. + /// + /// \param nThreads the number of execution threads EvolveState should use. + /// If zero, use an automatic, system defined number. + void Set_EvalThreads(unsigned int nThreads){ + evalThreads=nThreads; + if(evalThreads==0){ //if the user wants automatic selection + evalThreads=std::thread::hardware_concurrency(); + if(evalThreads==0) //if autodetection failed + evalThreads=1; + } + } + + /// \brief Get the number of threads used for evolution. + /// + /// \return The number of execution threads EvolveState should use. + unsigned int Get_EvalThreads() const{ + return(evalThreads); + } + + /// \brief Stes Earth object to be use. + /// @param earth Shared pointer to Earth object. + void Set_EarthModel(std::shared_ptr earth){ + earth_atm=earth; + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_Body(earth); + } + } + + /// \brief Returns the neutrino interaction cross sections + std::shared_ptr GetNeutrinoCrossSections(){ + return(ncs); + } + + /// \brief Sets the neutrino interaction cross sections + void SetNeutrinoCrossSections(std::shared_ptr xs){ + ncs=xs; + for(BaseSQUIDS& nsq : nusq_array) + nsq.SetNeutrinoCrossSections(ncs); + } + + void Set_EvolLowPassCutoff(double val){ + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_EvolLowPassCutoff(val); + } + } + + void Set_EvolLowPassScale(double val){ + for(BaseSQUIDS& nsq : nusq_array){ + nsq.Set_EvolLowPassScale(val); + } + } + + void Set_NeutrinoSources(bool opt){ + for(BaseSQUIDS& nsq : nusq_array) + nsq.Set_NeutrinoSources(opt); + } +}; + + +} // close namespace +#endif From f4eaa0b688ae7e7b0c8476a69b8f5eacd4999d67 Mon Sep 17 00:00:00 2001 From: ConnorSponsler <33155897+ConnorSponsler@users.noreply.github.com> Date: Sun, 19 Nov 2023 16:35:26 -0500 Subject: [PATCH 08/12] Add files via upload --- src/body.cpp | 26 ++++++++------------------ 1 file changed, 8 insertions(+), 18 deletions(-) diff --git a/src/body.cpp b/src/body.cpp index 70b78f31..a873a474 100644 --- a/src/body.cpp +++ b/src/body.cpp @@ -946,7 +946,7 @@ EmittingEarthAtm::EmittingEarthAtm(std::string filepath):EarthAtm() unique_sort(atm_coszen, czen_min, czen_max); - int czen_number = atm_coszen.size(); + long unsigned int czen_number = atm_coszen.size(); marray coszens; coszens.resize(0, czen_number); for (unsigned int i=0; i < czen_number;i++){ @@ -954,7 +954,7 @@ EmittingEarthAtm::EmittingEarthAtm(std::string filepath):EarthAtm() } unique_sort(atm_heights, height_min, height_max); - int height_number = atm_heights.size(); + long unsigned int height_number = atm_heights.size(); marray heights; heights.resize(0, height_number); for (unsigned int i=0; i < height_number; i++){ @@ -962,7 +962,7 @@ EmittingEarthAtm::EmittingEarthAtm(std::string filepath):EarthAtm() } unique_sort(atm_energy, energy_min, energy_max); - int energy_number = atm_energy.size(); + long unsigned int energy_number = atm_energy.size(); marray energies; energies.resize(0, energy_number); for (unsigned int i=0; i < energy_number; i++){ @@ -973,25 +973,15 @@ EmittingEarthAtm::EmittingEarthAtm(std::string filepath):EarthAtm() std::vector>> flux_matrices(num_prod_flav, std::vector>(2)); - marray nu_e_flux({static_cast(czen_number), - static_cast(height_number), - static_cast(energy_number)}); - marray nu_e_bar_flux({static_cast(czen_number), - static_cast(height_number), - static_cast(energy_number)}); - marray nu_mu_flux({static_cast(czen_number), - static_cast(height_number), - static_cast(energy_number)}); - marray nu_mu_bar_flux({static_cast(czen_number), - static_cast(height_number), - static_cast(energy_number)}); + 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({static_cast(czen_number), - static_cast(height_number), - static_cast(energy_number)}); + flux_matrices[flav][rho] = marray({czen_number, height_number, energy_number}); } } From 8886b60eac4352f364511394eb65ee83cd3c168e Mon Sep 17 00:00:00 2001 From: ConnorSponsler <33155897+ConnorSponsler@users.noreply.github.com> Date: Sun, 19 Nov 2023 16:41:15 -0500 Subject: [PATCH 09/12] Delete duplicate nuSQuIDS.h --- nuSQuIDS.h | 1951 ---------------------------------------------------- 1 file changed, 1951 deletions(-) delete mode 100644 nuSQuIDS.h diff --git a/nuSQuIDS.h b/nuSQuIDS.h deleted file mode 100644 index 74d08731..00000000 --- a/nuSQuIDS.h +++ /dev/null @@ -1,1951 +0,0 @@ - /****************************************************************************** - * 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 * - ******************************************************************************/ - - -#ifndef NUSQUIDS_nuSQUIDS_H -#define NUSQUIDS_nuSQUIDS_H - -#if __cplusplus < 201103L -#error C++11 compiler required. Update your compiler and use the flag -std=c++11 -#endif - -#include -#include -#include -#include -#include - -#include - -#include -#include -#include -#include - -#include -#include - -#include "nuSQuIDS/aligned_alloc.h" -#include "nuSQuIDS/body.h" -#include "nuSQuIDS/marray.h" -#include "nuSQuIDS/taudecay.h" -#include "nuSQuIDS/ThreadPool.h" -#include "nuSQuIDS/xsections.h" - -//#define FixCrossSections - -//#define UpdateInteractions_DEBUG - -namespace nusquids{ - -enum NeutrinoType { - neutrino=0b01, - antineutrino=0b10, - both=0b11 -}; - -enum Basis { - mass=0b01, - flavor=0b10, - interaction=0b11 -}; - - -///\class nuSQUIDS -///\brief nu-SQuIDS main class -class nuSQUIDS: public squids::SQuIDS { - // nuSQUIDSAtm is a friend class so it can use H0 to evolve operators - // and thus evaluate expectation values. - template - friend class nuSQUIDSAtm; -protected: - /// \brief Sets the basis in which the problem will be solved. - /// - /// If interaction basis is used the projectors will be evolved at - /// every time step. On the other hand, if mass basis is used no evolution - /// is performed. - Basis basis = interaction; - /// \brief number of neutrino flavors. - unsigned int numneu; - /// \brief number of energy nodes. - unsigned int ne = nx; - - ///Get the fraction of the of the number density of targets which is of each target type. - ///These fractions sum to one. - ///\return a vector of fractions, in the same order as the targets in int_struct - std::vector GetTargetNumberFractions() const; - - ///Get the number density of all nuclear targets at the current position - ///\return a vector of number densities, in the same order as the targets in int_struct - std::vector GetTargetNumberDensities() const; - - /// \brief Updates the interaction length arrays. - /// - /// Uses GetTargetNumberDensities() together with the stored cross section - /// information to update: nuSQUIDS#invlen_NC, nuSQUIDS#invlen_CC, and nuSQUIDS#invlen_INT. - void UpdateInteractions(); - - /// \brief Contains the energy nodes. - marray E_range; - - /// \brief Contains the spaces between nodes. - /// - /// It has length len(E_range)-1. - marray delE; - - /// \brief Interface that calculates and interpolates neutrino cross sections. - std::shared_ptr ncs; - - /// \brief Interface that calculate and interpolates tau decay spectral functions. - TauDecaySpectra tdc; - - /// \brief NT keeps track if the problem consists of neutrinos, antineutrinos, or both. - NeutrinoType NT = both; - - ///aligned arrays should be aligned to at least this many items - constexpr static size_t preferred_alignment=4; - constexpr static uint8_t log2(size_t n) { - return((n>1) ? log2(n/2)+1 : 0); - } - ///Round n up to the nearest multiple of preferred_alignment, doing nothing if n is already - ///such a multiple - constexpr static size_t round_up_to_aligned(size_t n){ - return((n&~(preferred_alignment-1))+((n&(preferred_alignment-1))?preferred_alignment:0)); - } - - public: - /// \brief Struct that contains all cross section information. - struct InteractionStructure { - /// \brief the configured target material types - /// \details These are recorded in the same order that they appear in dNdE_CC, - /// dNdE_NC, sigma_CC, and sigma_NC - std::vector targets; - /// \brief Neutrino charge current differential cross section with respect to - /// the outgoing lepton energy. - /// - /// This array is constructed when InitializeInteractionVectors() is called and - /// is initialized when InitializeInteractions() is called. Its dimensions are: - /// - target type (proton, neutron, etc.) - /// - neutrino type (neutrino, antineutrino) - /// - neutrino flavor (electron, muon, tau) - /// - initial energy - /// - final energy - marray> dNdE_CC; - /// \brief Neutrino neutral current differential cross section with respect to - /// the outgoing lepton energy. - /// - /// This array is constructed when InitializeInteractionVectors() is called and - /// is initialized when InitializeInteractions() is called. Its dimensions are: - /// - target type (proton, neutron, etc.) - /// - neutrino type (neutrino, antineutrino) - /// - neutrino flavor (electron, muon, tau) - /// - initial energy - /// - final energy - marray> dNdE_NC; - /// \brief Glashow resonance differential cross section (electron antineutrino only) - /// \details Dimensions are initial and final energy node - marray> dNdE_GR; - /// \brief Array that contains the neutrino charge current cross section. - /// \details The dimensions of this array are: - /// - target type (proton, neutron, etc.) - /// - neutrino type (neutrino, antineutrino) - /// - neutrino flavor (electron, muon, tau) - /// - energy - /// Its contents are in natural units, i.e. eV^-2. It is - /// initialized by InitializeInteractions() . - marray> sigma_CC; - /// \brief Array that contains the neutrino neutral current cross section. - /// \details The dimensions of this array are: - /// - target type (proton, neutron, etc.) - /// - neutrino type (neutrino, antineutrino) - /// - neutrino flavor (electron, muon, tau) - /// - energy - /// Its contents are in natural units, i.e. eV^-2. It is - /// initialized by InitializeInteractions() . - marray> sigma_NC; - /// \brief Glashow resonance cross section (electron antineutrino only) - /// \details 1 entry per energy node, in natural units - marray> sigma_GR; - /// \brief Array that contains the tau decay spectrum to all particles. - /// \details The first dimension corresponds to initial tau energy and the - /// second one to the outgoing lepton. - marray> dNdE_tau_all; - /// \brief Array that contains the tau decay spectrum to leptons. - /// \brief Array that contains the tau decay spectrum to all particles. - /// \details The first dimension corresponds to initial tau energy and the - /// second one to the outgoing lepton. - marray> dNdE_tau_lep; - /// \brief Default constructor - InteractionStructure(): - dNdE_CC(aligned_allocator{log2(preferred_alignment*sizeof(double))}), - dNdE_NC(aligned_allocator(log2(preferred_alignment*sizeof(double)))), - dNdE_GR(aligned_allocator(log2(preferred_alignment*sizeof(double)))), - sigma_CC(aligned_allocator(log2(preferred_alignment*sizeof(double)))), - sigma_NC(aligned_allocator(log2(preferred_alignment*sizeof(double)))), - sigma_GR(aligned_allocator(log2(preferred_alignment*sizeof(double)))), - dNdE_tau_all(aligned_allocator(log2(preferred_alignment*sizeof(double)))), - dNdE_tau_lep(aligned_allocator(log2(preferred_alignment*sizeof(double)))) - {} - /// \brief Assignment operator - InteractionStructure& operator=(const InteractionStructure& other){ - dNdE_CC=other.dNdE_CC; - dNdE_NC=other.dNdE_NC; - dNdE_GR=other.dNdE_GR; - sigma_CC=other.sigma_CC; - sigma_NC=other.sigma_NC; - sigma_GR=other.sigma_GR; - dNdE_tau_all=other.dNdE_tau_all; - dNdE_tau_lep=other.dNdE_tau_lep; - return(*this); - } - /// \brief Move assignment operator - InteractionStructure& operator=(InteractionStructure&& other){ - dNdE_CC=std::move(other.dNdE_CC); - dNdE_NC=std::move(other.dNdE_NC); - dNdE_GR=std::move(other.dNdE_GR); - sigma_CC=std::move(other.sigma_CC); - sigma_NC=std::move(other.sigma_NC); - sigma_GR=std::move(other.sigma_GR); - dNdE_tau_all=std::move(other.dNdE_tau_all); - dNdE_tau_lep=std::move(other.dNdE_tau_lep); - return(*this); - } - /// \brief Equality operator - bool operator==(const InteractionStructure& other) const{ - if ( - dNdE_CC.size() != other.dNdE_CC.size() or - dNdE_NC.size() != other.dNdE_NC.size() or - dNdE_GR.size() != other.dNdE_GR.size() or - sigma_CC.size() != other.sigma_CC.size() or - sigma_NC.size() != other.sigma_NC.size() or - sigma_GR.size() != other.sigma_GR.size() or - dNdE_tau_all.size() != other.dNdE_tau_all.size() or - dNdE_tau_lep.size() != other.dNdE_tau_lep.size() - ) - { - return false; - } - - if ( not std::equal(dNdE_CC.begin(),dNdE_CC.end(),other.dNdE_CC.begin()) ) - return false; - - if ( not std::equal(dNdE_NC.begin(),dNdE_NC.end(),other.dNdE_NC.begin()) ) - return false; - - if ( not std::equal(dNdE_GR.begin(),dNdE_GR.end(),other.dNdE_GR.begin()) ) - return false; - - if ( not std::equal(sigma_CC.begin(),sigma_CC.end(),other.sigma_CC.begin()) ) - return false; - - if ( not std::equal(sigma_NC.begin(),sigma_NC.end(),other.sigma_NC.begin()) ) - return false; - - if ( not std::equal(sigma_GR.begin(),sigma_GR.end(),other.sigma_GR.begin()) ) - return false; - - if ( not std::equal(dNdE_tau_all.begin(),dNdE_tau_all.end(),other.dNdE_tau_all.begin()) ) - return false; - - if ( not std::equal(dNdE_tau_lep.begin(),dNdE_tau_lep.end(),other.dNdE_tau_lep.begin()) ) - return false; - // all is good - lets roll - return true; - } - /// \brief Inequality operator - bool operator!=(const InteractionStructure& other) const { - return not (*this==other); - } - /// \brief Copy constructor operator - InteractionStructure(InteractionStructure& other): - dNdE_CC(other.dNdE_CC),dNdE_NC(other.dNdE_NC),dNdE_GR(other.dNdE_GR), - sigma_CC(other.sigma_CC),sigma_NC(other.sigma_NC),sigma_GR(other.sigma_GR), - dNdE_tau_all(other.dNdE_tau_all), - dNdE_tau_lep(other.dNdE_tau_lep) - {} - /// \brief Move constructor operator - InteractionStructure(InteractionStructure&& other): - dNdE_CC(std::move(other.dNdE_CC)),dNdE_NC(std::move(other.dNdE_NC)),dNdE_GR(std::move(other.dNdE_GR)), - sigma_CC(std::move(other.sigma_CC)),sigma_NC(std::move(other.sigma_NC)),sigma_GR(std::move(other.sigma_GR)), - dNdE_tau_all(std::move(other.dNdE_tau_all)), - dNdE_tau_lep(std::move(other.dNdE_tau_lep)) - {} - }; - - /// \brief Contains the interaction information at the current position. - struct InteractionState { - /// \brief Array that contains the inverse of the neutrino charge current mean free path. - /// \details The array contents are in natural units (i.e. eV) and is update when - /// UpdateInteractions() is called. The first dimension corresponds to the neutrino type, - /// the second to the flavor, and the last one to the energy. - marray> invlen_CC; - /// \brief Array that contains the inverse of the neutrino neutral current mean free path. - /// \details The array contents are in natural units (i.e. eV) and is update when - /// UpdateInteractions() is called. The first dimension corresponds to the neutrino type, - /// the second to the flavor, and the last one to the energy. - marray> invlen_NC; - /// \brief Glashow resonance inverse interaction length (electron antineutrino only) - /// \details 1 entry per energy node - marray> invlen_GR; - /// \brief Array that contains the inverse of the neutrino total mean free path. - /// \details The array contents are in natural units (i.e. eV) and is update when - /// UpdateInteractions() is called. Numerically it is just nuSQUIDS::invlen_NC and nuSQUIDS::invlen_CC - /// added together. - marray> invlen_INT; - - /// \brief Default constructor - InteractionState(): - invlen_CC(aligned_allocator(log2(preferred_alignment*sizeof(double)))), - invlen_NC(aligned_allocator(log2(preferred_alignment*sizeof(double)))), - invlen_GR(aligned_allocator(log2(preferred_alignment*sizeof(double)))), - invlen_INT(aligned_allocator(log2(preferred_alignment*sizeof(double)))) - {} - /// \brief Copy constructor - InteractionState(InteractionState& other): - invlen_CC(other.invlen_CC),invlen_NC(other.invlen_NC), - invlen_GR(other.invlen_GR),invlen_INT(other.invlen_INT) - {} - /// \brief Move constructor - InteractionState(InteractionState&& other): - invlen_CC(std::move(other.invlen_CC)),invlen_NC(std::move(other.invlen_NC)), - invlen_GR(std::move(other.invlen_GR)),invlen_INT(std::move(other.invlen_INT)) - {} - /// \brief Assignment operator - InteractionState& operator=(const InteractionState& other){ - invlen_NC=other.invlen_NC; - invlen_CC=other.invlen_CC; - invlen_GR=other.invlen_GR; - invlen_INT=other.invlen_INT; - return(*this); - } - /// \brief Move assignment operator - InteractionState& operator=(InteractionState&& other){ - invlen_NC=std::move(other.invlen_NC); - invlen_CC=std::move(other.invlen_CC); - invlen_GR=std::move(other.invlen_GR); - invlen_INT=std::move(other.invlen_INT); - return(*this); - } - /// \brief Equality operator - bool operator==(const InteractionState& other) const{ - if (invlen_CC.size() != other.invlen_CC.size() or - invlen_NC.size() != other.invlen_NC.size() or - invlen_GR.size() != other.invlen_GR.size() or - invlen_INT.size() != other.invlen_INT.size()) - return false; - - if ( not std::equal(invlen_CC.begin(),invlen_CC.end(),other.invlen_CC.begin()) ) - return false; - - if ( not std::equal(invlen_NC.begin(),invlen_NC.end(),other.invlen_NC.begin()) ) - return false; - - if ( not std::equal(invlen_GR.begin(),invlen_GR.end(),other.invlen_GR.begin()) ) - return false; - - if ( not std::equal(invlen_INT.begin(),invlen_INT.end(),other.invlen_INT.begin()) ) - return false; - - // all is good - lets roll - return true; - } - /// \brief Inequality operator - bool operator!=(const InteractionState& other) const { - return not (*this==other); - } - }; - - protected: - /// \brief Interaction struct that contains all cross section information tabulated. - std::shared_ptr int_struct; - /// \brief Current set of interaction information - InteractionState int_state; - - /// \brief Length upon which the neutrino fluxes will be positivized. - double positivization_scale; - - /// \brief Body where the neutrino propagation takes place. - std::shared_ptr body; - /// \brief Trayectory within the body. - /// \details Stores the position within the body and its updated every evolution - /// step. - std::shared_ptr track; - - /// \brief SU_vector that represents the neutrino square mass difference matrix in the mass basis. - /// It is used to construct nuSQUIDS#H0_array and H0() - squids::SU_vector DM2; - /// \brief Stores the time independent hamiltonian corresponding to each energy node. - marray H0_array; - - /// \brief Product of physical constants needed by HI. - double HI_constants; - /// \brief The density of the body at the current point on the track - double current_density; - /// \brief The electron fraction of the body at the current point on the track - double current_ye; - /// \brief The external flux from the body at the current point on the track - marray current_external_flux; - - /// \brief Mass basis projectors. - /// \details The i-entry corresponds to the projector in the ith mass eigenstate. - marray b0_proj; - /// \brief Flavor basis projectors. - /// \details The first dismension corresponds to the neutrino type. When NeutrinoType = both, then - /// the first dimension equal 0 means neutrinos and 1 antineutrinos. The second dimension corresponds - /// to the flavor eigenstates where 0 corresponds to nu_e, 1 to nu_mu, 2 to nu_tau, and the others - /// to sterile neutrinos. - marray b1_proj; - /// \brief Mass basis projectors in the interaction picture. - /// The index meaning are the same as nuSQUIDS#b1_proj but for mass eigenstates, - /// with an added third dimension that corresponds to the energy node index. - marray evol_b0_proj; - /// \brief Flavor basis projectors in the interaction picture. - /// The index meaning are the same as nuSQUIDS#b1_proj , - /// with an added third dimension that corresponds to the energy node index. - marray evol_b1_proj; - - ///Cache of precalculated results for InteractionsRho. - ///Used only when interactions are on and oscillations are off. - marray interaction_cache; - ///Backing storage for the SU_vectors in interaction_cache so that each can - ///be ensured to be correctly aligned for best performance. - std::unique_ptr interaction_cache_store; - ///Total size of interaction_cache_store - size_t interaction_cache_store_size; - - marray nc_factors; - marray tau_hadlep_decays; - marray tau_lep_decays; - marray gr_factors; - - ///Initialize the interaction cahce data structures to the corretc sizes for - ///the current problem. - ///Should be invoked only when interactions are on and oscillations are off. - void SetUpInteractionCache(); - - /// \brief Evolves the flavor projectors in the interaction basis to a time t. - /// \details It uses H0() to evolve SQUIDS#b0_proj and SQUIDS#b1_proj into - /// SQUIDS#evol_b0_proj and SQUIDS#evol_b1_proj. - /// \warning Since the RHS of the differential equation only involves flavor projectors - /// we do not current evolve mass projectors. - void EvolveProjectors(double t); - - // bool requirements - private: - /// \brief Boolean that signals that the full hamiltonian including extra potentials should be used for - /// the projector evolution - bool use_full_hamiltonian_for_projector_evolution = false; - /// \brief Boolean that toggles debug information to be printed - bool debug = false; - /// \brief Boolean that signals the object correct initialization. - bool inusquids = false; - /// \brief Boolean that signals that a Body object has being set. - bool ibody = false; - /// \brief Boolean that signals that the neutrino energies has being set. - bool ienergy = false; - /// \brief Boolean that signals that a Track object has being set. - bool itrack = false; - /// \brief Boolean that signals that an initial state has being set. - bool istate = false; - /// \brief Boolean that signals that interactions will be taken into account. - bool iinteraction = false; - /// \brief Whether interactions have been initialized - bool interactions_initialized = false; - /// \brief Boolean that signals that neutrino oscillations will be taken into account. - bool ioscillations = true; - /// \brief Boolean that signals that tau regeneration is being used. - bool tauregeneration = false; - /// \brief Whether the Glashsow resonance for electron antineutrinos is used. - bool iglashow = false; - /// \brief Boolean that signals that positivization will be enforced. - bool positivization = false; - /// \brief Boolean that allows to use fast constant density evolution technique. - bool allowConstantDensityOscillationOnlyEvolution = false; - /// \brief Boolean that signals that a progress bar will be printed. - bool progressbar = false; - /// \brief Integer to keep track of the progress bar evolution. - int progressbar_count = 0; - /// \brief Number of steps upon which the progress bar will be updated. - int progressbar_loop = 100; - /// \brief Time offset between SQuIDS time and Track(x). - double time_offset; - /// \brief Force flavor projections to be positive. - void PositivizeFlavors(); - /// \brief Set GSL differential cross section precision. - double gsl_int_precision = 1.e-3; - /// \brief Cutoff for low-pass filter applied during state density evolution - double evol_lowpass_cutoff = 0; - /// \brief Scale for low-pass filter applied during state density evolution - double evol_lowpass_scale = 0; - /// \brief If true the bodies can act as neutrino sources. - bool enable_neutrino_sources = false; - protected: - /// \brief Initializes flavor and mass projectors - /// \warning Antineutrinos are handle by means of the AntineutrinoCPFix() function - /// which temporary flips the CP phase. - void iniProjectors(); - /// \brief Reinitializes the flavor projectors. - /// \details It is called before evolving the system and every time the system - /// state is going to be set, since the definition of mass and flavor basis might - /// have changed. - void SetIniFlavorProyectors(); - /// \brief Initializes the time independent hamiltonian for each energy node. - /// \details It constructs nuSQUIDS#DM2 and nuSQUIDS#H0_array - void iniH0(); - /// \brief Changes the CP sign of the complex matrices stored in params. - void AntineutrinoCPFix(unsigned int irho); - /// \brief Serializes the initialization of the Body and Track objects. - /// @see ReadStateHDF5 - void SetBodyTrack(unsigned int,unsigned int,double*,unsigned int,double*); - - /// \brief General initilizer for the multi energy mode - /// @param E_vector Energy nodes [eV]. - /// @param xini The initial position of the system. - void init(marray E_vector,double xini=0.0); - - /// \brief Initilizer for the single energy mode - /// @param xini The initial position of the system. By default is set to 0. - void init(double xini=0.0); - - /// \brief Initilizes auxiliary cross section arrays. - /// \details The arrays are initialize, but not filled with contented. - /// To fill the arrays call InitializeInteractions(). - /// \param nTargets the number of target species for which to allocate space - void InitializeInteractionVectors(unsigned int nTargets); - - /// \brief Fills in auxiliary cross section arrays. - /// \details It uses nuSQUIDS#ncs and nuSQUIDS#tdc to fill in the values of - /// nuSQUIDS#dNdE_CC , nuSQUIDS#dNdE_NC , nuSQUIDS#sigma_CC , nuSQUIDS#sigma_NC , - /// nuSQUIDS#invlen_CC , nuSQUIDS#invlen_NC , nuSQUIDS#invlen_INT , - /// nuSQUIDS#dNdE_tau_all , nuSQUIDS#dNdE_tau_lep - /// in int_struct. - /// @see InitializeInteractionVectors - void GetCrossSections(); - private: - - /// \brief Prints progress bar. - /// \details To enable it call Set_ProgressBar() - void ProgressBar() const; - protected: - /// \brief Updates all quantities needed during the evaluation of the RHS. - /// @param x Position of the system. - /// \details Dependind on the basis used it evolves the flavor projectors - /// by means of EvolveProjectors(). Also, when interactions are considered, - /// it updates the interaction arrays by calling UpdateInteractions(). Finally, - /// it handles control to the user implemented updates via AddToPreDerive(). - /// @see AddToPreDerive - /// @see EvolveProjectors - /// @see UpdateInteractions - void PreDerive(double x); - /// \brief User supplied function that is called before performing a derivative. - /// @param x Position of the system. - /// @see PreDerive - virtual void AddToPreDerive(double x){} - - /// \brief Reads and constructs the object from an HDF5 file. - /// @param hdf5_filename Filename of the HDF5 to use for construction. - /// @param group Path to the group where the nuSQUIDS content will be saved. - /// @param iis InteractionStructure that contains valid cross sections. - /// \details By default all contents are assumed to be in the \c root of the HDF5 file - /// if the user wants a different location it can use \c group and \c save_cross_sections to - /// change the nuSQUIDS location and the cross section information location. - /// Finally, it calls AddToReadHDF5() to enable user defined parameters - /// to be read. - /// @see AddToReadHDF5 - /// @see WriteStateHDF5 - void ReadStateHDF5Internal(std::string hdf5_filename,std::string group = "/", std::shared_ptr iis = nullptr); - - public: - /************************************************************************************ - * CONSTRUCTORS - *************************************************************************************/ - - /// \brief Default void constructor. - nuSQUIDS(){} - - /// \brief Move constructor. - nuSQUIDS(nuSQUIDS&&); - - /// \brief Multiple energy mode constructor. - /// @param E_vector Energy nodes [eV]. - /// @param numneu Number of neutrino flavors. - /// @param NT NeutrinoType: neutrino,antineutrino, or both (simultaneous solution). - /// @param iinteraction Sets the neutrino noncoherent neutrino interactions on. - /// @param ncs Cross section object. - /// \details By default the energy scale is logarithmic and interactions are turn off. - /// \warning When interactions are present interpolation is performed to precalculate the neutrino - /// cross section which make take considertable time depending on the energy grid. - /// @see init - nuSQUIDS(marray E_vector,unsigned int numneu,NeutrinoType NT = both, - bool iinteraction = false, std::shared_ptr ncs = nullptr): - numneu(numneu),ncs(ncs),NT(NT),iinteraction(iinteraction) - {init(E_vector);} - - /// \brief Single energy mode constructor. - /// @param numneu Number of neutrino flavors. - /// @param NT NeutrinoType: neutrino or antineutrino. - /// \warning Interactions are not possible in the single energy mode, nor is simultaneous - /// neutrino-antineutrino solution (both) possible. - /// \details Constructors projectors and initializes Hamiltonian. - nuSQUIDS(unsigned int numneu, NeutrinoType NT = neutrino): - numneu(numneu),NT(NT),iinteraction(false) - {init();} - - /// \brief Constructor from a HDF5 filepath. - /// @param hdf5_filename Filename of the HDF5 to use for construction. - /// @param grp HDF5 file group path where the nuSQUIDS object is saved. - /// @param int_struct Interaction Structure to be used in object construction, - /// if this has already been obtained. May be left as NULL to cause - /// the interaction information to be read from the standard location. - /// \details Reads the HDF5 file and constructs the associated nuSQUIDS object - /// restoring all properties as well as the state. - /// @see ReadStateHDF5 - nuSQUIDS(std::string hdf5_filename, std::string grp = "/", - std::shared_ptr int_struct = nullptr): - int_struct(int_struct) - { - if(this->int_struct) - ReadStateHDF5Internal(hdf5_filename,grp,int_struct); - else - ReadStateHDF5(hdf5_filename, grp, ""); - } - - //*************************************************************** - virtual ~nuSQUIDS(); - - //*************************************************************** - ///\brief Move assigns a nuSQUIDS object from an existing object - nuSQUIDS& operator=(nuSQUIDS&&); - - protected: - /************************************************************************************ - * PHYSICS FUNCTIONS - SUPER IMPORTANT - *************************************************************************************/ - - /// \brief Returns the time independent part of the Hamiltonian - /// @param E Neutrino energy [eV] - /// @param irho Density matrix equation index. - /// \details \c irho is the index of the equation on which the Hamiltonian - /// is used. When not in NeutrinoType == \c both \c irho is only 0, but if - /// neutrinos and antineutrinos are solved simultaneously, then 0 is for - /// neutrinos and 1 for antineutrinos. - squids::SU_vector H0(double E, unsigned int irho) const; - - /// \brief Returns the time dependent part of the Hamiltonian at an energy node \c ie - /// @param ie Energy node - /// @param irho Density matrix equation index. - /// @see H0 for convention on \c irho. - virtual squids::SU_vector HI(unsigned int ie, unsigned int irho) const; - - /// \brief Returns noncoherent interaction term in the Hamiltonian at node \c ie - /// @param ie Energy node - /// @param irho Density matrix equation index. - /// @see H0 for convention on \c irho. - virtual squids::SU_vector GammaRho(unsigned int ie, unsigned int irho) const; - - /// \brief Returns the scalar quantity decay rate at an energy node \c ie - /// @param ie Energy node - /// @param iscalar scalar equation index. - /// \details When tau regenereation is considered \c iscalar = 0 corresponds - /// to the tau flux. - virtual double GammaScalar(unsigned int ie, unsigned int iscalar) const; - - /// \brief Returns interaction of the density matrix at an energy node \c ie - /// @param ie Energy node - /// @param irho Density matrix equation index. - /// @see H0 for convention on \c irho. - virtual squids::SU_vector InteractionsRho(unsigned int ie, unsigned int irho) const; - - /// \brief Returns scalar interactions at an energy node \c ie - /// @param ie Energy node - /// @param iscalar scalar equation index. - /// @see GammaScalar for convention on \c iscalar. - virtual double InteractionsScalar(unsigned int ie, unsigned int iscalar) const; - private: - /// \brief SQuIDS signature of HI - squids::SU_vector HI(unsigned int ie, unsigned int irho, double x) const {return HI(ie,irho);} - /// \brief SQuIDS signature of GammaRho - squids::SU_vector GammaRho(unsigned int ei, unsigned int irho, double x) const {return GammaRho(ei,irho);} - /// \brief SQuIDS signature of GammaScalar - double GammaScalar(unsigned int ei, unsigned int iscalar,double x) const {return GammaScalar(ei,iscalar);} - /// \brief SQuIDS signature of InteractionsRho - squids::SU_vector InteractionsRho(unsigned int ei, unsigned int irho, double x) const {return InteractionsRho(ei,irho);} - /// \brief SQuIDS signature of InteractionsScalar - double InteractionsScalar(unsigned int ei, unsigned int irho, double x) const {return InteractionsScalar(ei,irho);} - public: - /************************************************************************************ - * PUBLIC MEMBERS TO EVALUATE/SET/GET STUFF - *************************************************************************************/ - - /// \brief Sets the initial state in the single energy mode - /// @param ini_state Initial neutrino state. - /// @param basis Representation of the neutrino state either flavor or mass. - /// \details \c ini_state length has to be equal to \c numneu. If the basis is - /// flavor then the entries are interpret as nu_e, nu_mu, nu_tau, nu_sterile_1, ..., nu_sterile_n, - /// while if the mass basis is used then the first entries correspond to the active - /// mass eigenstates. - void Set_initial_state(const marray& ini_state, Basis basis = flavor); - - /// \brief Sets the initial state in the multiple energy mode when doing either neutrino or antineutrino only. - /// @param ini_state Initial neutrino state. - /// @param basis Representation of the neutrino state either flavor or mass. - /// \details \c ini_state first dimension has length equal to \c ne (number of energy nodes), while - /// the second dimension has length equal to \c numneu (number of flavors). If the basis is - /// flavor then the entries are interpret as nu_e, nu_mu, nu_tau, nu_sterile_1, ..., nu_sterile_n, - /// while if the mass basis is used then the first entries correspond to the active - /// mass eigenstates. - void Set_initial_state(const marray& ini_state, Basis basis = flavor); - - /// \brief Sets the initial state in the multiple energy mode when doing both neutrino and antineutrino. - /// @param ini_state Initial neutrino state. - /// @param basis Representation of the neutrino state either flavor or mass. - /// \details \c ini_state first dimension has length equal to \c ne (number of energy nodes). The second - /// dimension has to have size 2 and the first entry is for neutrinos while the second one for antineutrinos - /// Finally. the third dimension has length equal to \c numneu (number of flavors). If the basis is - /// flavor then the entries are interpret as nu_e, nu_mu, nu_tau, nu_sterile_1, ..., nu_sterile_n, - /// while if the mass basis is used then the first entries correspond to the active - /// mass eigenstates. - void Set_initial_state(const marray& ini_state, Basis basis = flavor); - - /// \brief Sets the body where the neutrino will propagate. - /// @param body Body object. - void Set_Body(std::shared_ptr body); - - /// \brief Sets neutrino trajectory - /// @param track Trajectory corresponding to the Body. - void Set_Track(std::shared_ptr track); - - /// \brief Sets the neutrino energy in single energy mode. - /// @param enu Neutrino energy [eV] - /// @pre NeutrinoType must be neutrino or antineutrino, not both. - void Set_E(double enu); - - /// \brief Evolves the system. - /// @pre Body, track, neutrino state, and energy must has being set - /// before calling this function. - void EvolveState(); - - /// \brief Returns the mass composition at a given node. - /// @param flv Neutrino flavor. - /// @param ie Energy node index. - /// @param rho Index of the equation, see details. - /// \details When NeutrinoType is \c both \c rho specifies wether one - /// is considering neutrinos (0) or antineutrinos (1). - - double EvalMassAtNode(unsigned int flv,unsigned int ie,unsigned int rho = 0) const; - /// \brief Returns the flavor composition at a given node. - /// @param flv Neutrino flavor. - /// @param ie Energy node index. - /// @param rho Index of the equation, see details. - /// \details When NeutrinoType is \c both \c rho specifies wether one - /// is considering neutrinos (0) or antineutrinos (1). - double EvalFlavorAtNode(unsigned int flv,unsigned int ie,unsigned int rho = 0) const; - - /// \brief Returns the mass composition at a given energy in the multiple energy mode. - /// averaging out the high frequencies. - /// @param flv Neutrino flavor. - /// @param enu Neutrino energy in natural units [eV]. - /// @param rho Index of the equation, see details. - double EvalMass(unsigned int flv,double enu,unsigned int rho) const; - - /// \brief Returns the mass composition at a given energy in the multiple energy mode. - /// averaging out the high frequencies. - /// @param flv Neutrino flavor. - /// @param enu Neutrino energy in natural units [eV]. - /// @param rho Index of the equation, see details. - /// @param scale scale upon which oscillations will be averaged out - /// @param avg bool array which is true for all scales that were averaged out - double EvalMass(unsigned int flv,double enu,unsigned int rho,double scale,std::vector& avr) const; - - /// \brief Returns the flavor composition at a given energy in the multiple energy mode. - /// @param flv Neutrino flavor. - /// @param enu Neutrino energy in natural units [eV]. - /// @param rho Index of the equation, see details. - /// \details When NeutrinoType is \c both \c rho specifies wether one - /// is considering neutrinos (0) or antineutrinos (1). - double EvalFlavor(unsigned int flv,double enu,unsigned int rho = 0) const; - - /// \brief Returns the flavor composition at a given energy in the multiple energy mode. - /// averaging out the high frequencies. - /// @param flv Neutrino flavor. - /// @param enu Neutrino energy in natural units [eV]. - /// @param rho Index of the equation, see details. - /// @param scale scale upon which oscillations will be averaged out - /// @param avg bool array which is true for all scales that were averaged out - double EvalFlavor(unsigned int flv,double enu,unsigned int rho,double scale,std::vector& avr) const; - - /// \brief Returns the mass composition in the single energy mode. - /// @param flv Neutrino flavor. - /// @param scale scale upon which oscillations will be averaged out - /// @param avg bool array which is true for all scales that were averaged out - double EvalMass(unsigned int flv, double scale, std::vector& avr) const; - - /// \brief Returns the flavor composition in the single energy mode. - /// @param flv Neutrino flavor. - /// @param scale scale upon which oscillations will be averaged out - /// @param avg bool array which is true for all scales that were averaged out - double EvalFlavor(unsigned int flv, double scale, std::vector& avr) const; - - /// \brief Returns the mass composition in the single energy mode. - /// @param flv Neutrino flavor. - double EvalMass(unsigned int flv) const; - - /// \brief Returns the flavor composition in the single energy mode. - /// @param flv Neutrino flavor. - double EvalFlavor(unsigned int flv) const; - - /// \brief Sets the cutoff for the state evolution low-pass filter. - /// @param val cutoff value - void Set_EvolLowPassCutoff(double val); - - /// \brief Sets the linear ramp size for the state evolution low-pass filter. - /// @param val Range in frequency space over which the linear ramp is applied. - void Set_EvolLowPassScale(double val); - - /// \brief Toggles tau regeneration on and off. - /// \param opt If \c true tau regeneration will be considered. - void Set_TauRegeneration(bool opt); - - /// \brief Toggles the effect of the Glashow resonance on and off. - /// \param opt If \c true resonant W^- production will be considered. - void Set_GlashowResonance(bool opt); - - /// \brief Toggles neutrino oscillations on and off. - /// \param opt If \c true neutrino oscillations will be considered. - void Set_IncludeOscillations(bool opt); - - /// \brief Toggles if non-hard processes constant density fast evolution will be used. - /// \param opt If \c true fast evolution will be attempted on constatnt density. - void Set_AllowConstantDensityOscillationOnlyEvolution(bool opt); - - /// \brief Toggles positivization of the flux. - /// @param opt If \c true the flux will be forced to be positive every \c positivization_step. - void Set_PositivityConstrain(bool opt); - - /// \brief Sets the positivization step. - /// @param step Sets the positivization step. - void Set_PositivityConstrainStep(double step); - - /// \brief Toggles the progress bar printing on and off - /// @param opt If \c true a progress bar will be printed. - void Set_ProgressBar(bool opt); - - /// \brief Returns the energy nodes values. - marray GetERange() const; - - /// \brief Returns number of energy nodes. - size_t GetNumE() const; - - /// \brief Returns the number of neutrino flavors. - unsigned int GetNumNeu() const; - - /// \brief Returns the number of rho equations. - unsigned int GetNumRho() const; - - /// \brief Returns true if noncoherent interactions are considered - bool GetUseInteractions() const {return iinteraction;} - - /// \brief Returns true if oscillations are considered - bool GetUseOscillations() const {return ioscillations;} - - /// \brinf Initiàlize interaction structure - void InitializeInteractions(); - - /// \brief Returns the interaction structure. - std::shared_ptr GetInteractionStructure() { - if(!interactions_initialized) - throw std::logic_error("Interactions not initialized"); - return int_struct; - } - - /// \brief Returns the interaction structure with constness. - std::shared_ptr GetInteractionStructure() const { - if(!interactions_initialized) - throw std::logic_error("Interactions not initialized"); - return int_struct; - } - - /// \brief Return the Hamiltonian at the current time - squids::SU_vector GetHamiltonian(unsigned int ei, unsigned int rho = 0); - /// \brief Returns the state - const squids::SU_vector& GetState(unsigned int ei,unsigned int rho = 0) const{ - return state[ei].rho[rho]; - } - /// \brief Returns the flavor projector - const squids::SU_vector& GetFlavorProj(unsigned int flv,unsigned int rho = 0) const{ - return b1_proj[rho][flv]; - } - /// \brief Returns the mass projector - const squids::SU_vector& GetMassProj(unsigned int flv,unsigned int rho = 0) const{ - return b0_proj[flv]; - } - - /// \brief Returns the trajectory object. - std::shared_ptr GetTrack(); - /// \brief Returns the body object. - std::shared_ptr GetBody(); - - const Track& GetTrackFast() const{ return(*track); } - - /// \brief Writes the object into an HDF5 file. - /// @param hdf5_filename Filename of the HDF5 to use for construction. - /// @param group Path to the group where the nuSQUIDS content will be saved. - /// @param save_cross_sections If \c true the cross section tables will also be saved. - /// @param cross_section_grp_loc Path to the group where the cross section will be saved. - /// @param overwrite If true will overwrite previous file, else it will append. - /// \details By default all contents are saved to the \c root of the HDF5 file - /// if the user wants a different location it can use \c group and \c save_cross_sections to - /// change the nuSQUIDS location and the cross section information location. Furthermore, the - /// \c save_cross_sections flag, which defaults to \c false decides if the cross section - /// will be saved. Finally, it calls AddToWriteHDF5() to enable user defined parameters - /// to be saved. - /// @see AddToWriteHDF5 - /// @see ReadStateHDF5 - void WriteStateHDF5(std::string hdf5_filename,std::string group = "/",bool save_cross_sections = true, std::string cross_section_grp_loc = "",bool overwrite = true) const; - - /// \brief User function to write user defined properties from an HDF5 file. - /// @param hdf5_loc_id HDF5 group id - /// \details This function will be called by WriteStateHDF5() which will provide - /// the hdf5_loc_id enabling the user to read user defined contents. - /// @see WriteStateHDF5 - /// @see AddToReadHDF5 - virtual void AddToWriteHDF5(hid_t hdf5_loc_id) const; - - /// \brief Reads and constructs the object from an HDF5 file. - /// @param hdf5_filename Filename of the HDF5 to use for construction. - /// @param group Path to the group where the nuSQUIDS content will be saved. - /// @param cross_section_grp_loc Path to the group where the cross section will be saved. - /// \details By default all contents are assumed to be in the \c root of the HDF5 file - /// if the user wants a different location it can use \c group and \c save_cross_sections to - /// change the nuSQUIDS location and the cross section information location. - /// Finally, it calls AddToReadHDF5() to enable user defined parameters - /// to be read. - /// @see AddToReadHDF5 - /// @see WriteStateHDF5 - void ReadStateHDF5(std::string hdf5_filename,std::string group = "/", std::string cross_section_grp_loc = ""); - - - /// \brief User function to read user defined properties from an HDF5 file. - /// @param hdf5_loc_id HDF5 group id - /// \details This function will be called by ReadStateHDF5() which will provide - /// the hdf5_loc_id enabling the user to read user defined contents. - /// @see ReadStateHDF5 - /// @see AddToWriteHDF5 - virtual void AddToReadHDF5(hid_t hdf5_loc_id); - - /// \brief Sets the mixing angle th_ij. - /// @param i the (zero-based) index of the first state - /// @param j the (zero-based) index of the second state - /// must be larger than \c i. - /// @param angle Angle to use in radians. - /// \details Sets the neutrino mixing angle. In our zero-based convention, e.g., the th_12 is i = 0, j = 1.,etc. - virtual void Set_MixingAngle(unsigned int i, unsigned int j,double angle); - - /// \brief Returns the mixing angle th_ij in radians. - /// @param i the (zero-based) index of the first state - /// @param j the (zero-based) index of the second state - /// must be larger than \c i. - /// \details Gets the neutrino mixing angle. In our zero-based convention, e.g., the th_12 is i = 0, j = 1.,etc. - double Get_MixingAngle(unsigned int i, unsigned int j) const; - - /// \brief Sets the CP phase for the ij-rotation. - /// @param i the (zero-based) index of the first state - /// @param j the (zero-based) index of the second state - /// must be larger than \c i. - /// @param angle Phase to use in radians. - /// \details Sets the CP phase for the ij-rotation. In our zero-based convention, e.g., the delta_13 = delta_CP is i = 0, j = 2.,etc. - virtual void Set_CPPhase(unsigned int i, unsigned int j,double angle); - - /// \brief Returns the CP phase of the ij-rotation in radians. - /// @param i the (zero-based) index of the first state - /// @param j the (zero-based) index of the second state - /// must be larger than \c i. - /// \details Gets the CP phase for the ij-rotation. In our zero-based convention, e.g., the delta_13 = delta_CP is i = 0, j = 2.,etc. - double Get_CPPhase(unsigned int i, unsigned int j) const; - - /// \brief Sets the square mass difference with respect to first mass eigenstate - /// @param i the (zero-based) index of the second state - /// must be larger than \c 0. - /// @param sq Square mass difference in eV^2. - /// \details Sets square mass difference with respect to the first mass eigenstate. In our zero-based convention, e.g., the \f$\Delta m^2_{12}\f$ corresponds to (i = 1),etc. - virtual void Set_SquareMassDifference(unsigned int i, double sq); - - /// \brief Returns the square mass difference between state-i and the first mass eigenstate. - /// @param i the (zero-based) index of the first state - /// must be larger than \c i. - /// \details Returns square mass difference with respect to the first mass eigenstate. In our zero-based convention, e.g., the \f$\Delta m^2_{12}\f$ corresponds to (i = 1),etc. - double Get_SquareMassDifference(unsigned int i) const; - - /// \brief Sets the mixing parameters to default. - void Set_MixingParametersToDefault(); - - /// \brief Sets the basis on which the evolution will be perfomed. - /// @param basis Evolution basis can be either \c mass or \c interaction. - /// \details By detault the solution is perfomed in the interaction basis, which - /// is recommended for most problems. One might consider the user of the mass - /// basis when collective effects, defined by new - /// interactions introduce phases that cannot be cancelled as one goes to the - /// interaction basis. - void Set_Basis(Basis basis); - - /// \brief Sets GSL differential cross section integrator precision - /// @param gsl_int_precision Double that specifies the differential cross section integrator precision. - void Set_GSL_DifferentialCrossSection_Integrator_Precision(double gsl_int_precision_){ - gsl_int_precision = gsl_int_precision_; - } - - /// \brief Returns the GSL differential cross section integrator precision - double Get_GSL_DifferentialCrossSection_Integrator_Precision() const { - return gsl_int_precision; - } - - /// \brief Sets the debug information on/off - /// @param debug_ Boolean that toggles debuging on and off. - void Set_Debug(bool debug_){ - debug = debug_; - } - - int Get_DebugCounter() const{ return progressbar_count; } - - /// \brief Returns lepton mixing matrix - std::unique_ptr GetTransformationMatrix() const { - return params.GetTransformationMatrix(GetNumNeu()); - } - - /// \brief Returns the neutrino interaction cross sections - std::shared_ptr GetNeutrinoCrossSections() { - return(ncs); - } - - /// \brief Returns the neutrino interaction cross sections - void SetNeutrinoCrossSections(std::shared_ptr xs) { - ncs=xs; - interactions_initialized=false; - } - - /// \brief Enables neutrino flux emission from bodies - void Set_NeutrinoSources(bool enable_neutrino_sources_){ - if(enable_neutrino_sources_){ - current_external_flux.resize(std::vector{ne,nrhos,numneu}); - std::fill(current_external_flux.begin(),current_external_flux.end(),0.0); - } - if(not iinteraction and enable_neutrino_sources_) - Set_OtherRhoTerms(true); - enable_neutrino_sources = enable_neutrino_sources_; - } - - /// \brief Returns true if neutrino flux emission from bodies is considered - bool Get_NeutrinoSources(){ - return enable_neutrino_sources; - } -}; - -/** - * The following class provides functionalities - * for atmospheric neutrino experiments - * where a collection of trajectories is explored. - */ - -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 = BaseNusType; - private: - /// \brief Random number generator - gsl_rng * r_gsl; - // /// \brief Mixing matrix - // std::unique_ptr U{nullptr,[](gsl_matrix_complex*m){delete m;}}; - /// \brief Internal units - const squids::Const units; - /// \brief Boolean that signals that a progress bar will be printed. - bool progressbar = false; - /// \brief Linear interpolator. - double LinInter(double x,double xM, double xP, double yM, double yP) const { - double f2=(x-xM)/(xP-xM); - double f1=1-f2; - return f1*yM + f2*yP; - } - /// \brief Boolean that signals that an initial state has being set. - bool iinistate; - /// \brief Boolean that signals the object correct initialization. - bool inusquidsatm; - /// \brief Contains the cos(zenith) nodes. - marray costh_array; - /// \brief Contains the energy nodes [eV]. - marray enu_array; - /// \brief Contains the log of energy nodes (in log of eV). - marray log_enu_array; - /// \brief Contains the nuSQUIDS objects for each zenith. - std::vector nusq_array; - - /// \brief Contains the Earth in atmospheric configuration. - 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 - std::shared_ptr ncs; - /// \brief Contains the interaction information structure. - std::shared_ptr int_struct; - /// \brief the number of threads to use when performing the evolution - unsigned int evalThreads; - public: - /************************************************************************************ - * CONSTRUCTORS - *************************************************************************************/ - - /// \brief Basic constructor. - /// @param costh_array One dimensional array containing zenith angles to be calculated. - /// @param energy_min Minimum neutrino energy value [ev]. - /// @param energy_max Maximum neutrino energy value [eV]. - /// @param energy_div Number of energy divisions. - /// @param numneu Number of neutrino flavors. - /// @param NT Signals the neutrino type : neutrino, antineutrion or both (simultaneous solution) - /// @param iinteraction Sets the neutrino noncoherent neutrino interactions on. - /// \details By defaults interactions are not considered and the neutrino energy scale is assume logarithmic. - template - nuSQUIDSAtm(marray costh_array, ArgTypes&&... args): - costh_array(costh_array), - evalThreads(1) - { - gsl_rng_env_setup(); - const gsl_rng_type * T_gsl = gsl_rng_default; - r_gsl = gsl_rng_alloc (T_gsl); - - 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]); - } - - 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); }); - for(unsigned int ie = 0 ; ie < enu_array.size(); ie++) - log_enu_array[ie] = log(enu_array[ie]); - - inusquidsatm = true; - } - - /// \brief Constructor from a HDF5 filepath. - /// @param hdf5_filename Filename of the HDF5 to use for construction. - /// \details Reads the HDF5 file and construct the associated nuSQUIDSAtm object - /// restoring all properties as well as the state. - /// @see ReadStateHDF5 - nuSQUIDSAtm(std::string hdf5_filename) { - gsl_rng_env_setup(); - const gsl_rng_type * T_gsl = gsl_rng_default; - r_gsl = gsl_rng_alloc (T_gsl); - ReadStateHDF5(hdf5_filename); - //U = nusq_array[0].GetTransformationMatrix(); - } - - /// \brief Move constructor. - nuSQUIDSAtm(nuSQUIDSAtm&& other): - progressbar(other.progressbar), - iinistate(other.iinistate), - inusquidsatm(other.inusquidsatm), - costh_array(std::move(other.costh_array)), - enu_array(std::move(other.enu_array)), - log_enu_array(std::move(other.log_enu_array)), - nusq_array(std::move(other.nusq_array)), - earth_atm(std::move(other.earth_atm)), - track_array(std::move(other.track_array)), - ncs(std::move(other.ncs)), - int_struct(std::move(other.int_struct)), - evalThreads(other.evalThreads) - { - other.inusquidsatm = false; - } - - //*************************************************************** - ///\brief Move assigns a nuSQUIDSAtm object from an existing object - nuSQUIDSAtm& operator=(nuSQUIDSAtm&& other){ - if(&other==this) - return(*this); - - progressbar = other.progressbar; - iinistate = other.iinistate; - inusquidsatm = other.inusquidsatm; - costh_array = std::move(other.costh_array); - enu_array = std::move(other.enu_array); - log_enu_array = std::move(other.log_enu_array); - nusq_array = std::move(other.nusq_array); - earth_atm = std::move(other.earth_atm); - track_array = std::move(other.track_array); - ncs = std::move(other.ncs); - int_struct = std::move(other.int_struct); - evalThreads = other.evalThreads; - - // initial nusquids object render useless - other.inusquidsatm = false; - - return(*this); - } - /************************************************************************************ - * 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. - /// @param basis Representation of the neutrino state either flavor or mass. - /// \details \c ini_state first dimension lenght has to be equal to the number of - /// zenith bins, the second dimension is the number of energy bins, and the third - /// dimension corresponds to the number of neutrino flavors. If the basis is - /// flavor then the entries are interpret as nu_e, nu_mu, nu_tau, nu_sterile_1, ..., nu_sterile_n, - /// while if the mass basis is used then the first entries correspond to the active - /// mass eigenstates. - void Set_initial_state(const marray& ini_flux, Basis basis=flavor){ - if(ini_flux.extent(0) != costh_array.extent(0)) - throw std::runtime_error("nuSQUIDSAtm::Error::First dimension of input array is incorrect."); - if(ini_flux.extent(1) != enu_array.extent(0)) - throw std::runtime_error("nuSQUIDSAtm::Error::Second dimension of input array is incorrect."); - unsigned int i = 0; - for(BaseSQUIDS& nsq : nusq_array){ - marray slice{ini_flux.extent(1),ini_flux.extent(2)}; - for(size_t j=0; j& ini_flux, Basis basis=flavor){ - if(ini_flux.extent(0) != costh_array.extent(0)) - throw std::runtime_error( - "nuSQUIDSAtm::Error::First dimension of input array is incorrect."); - unsigned int i = 0; - for(BaseSQUIDS& nsq : nusq_array){ - marray slice{ini_flux.extent(1),ini_flux.extent(2),ini_flux.extent(3)}; - for(size_t j=0; jSet_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(){ - if(not iinistate) - throw std::runtime_error("nuSQUIDSAtm::Error::State not initialized."); - if(not inusquidsatm) - throw std::runtime_error("nuSQUIDSAtm::Error::nuSQUIDSAtm not initialized."); - - if(nusq_array.front().GetUseInteractions()){ - // setting the interaction structure also on the main object - nusq_array.front().InitializeInteractions(); - int_struct = nusq_array.front().GetInteractionStructure(); - for(BaseSQUIDS& nsq : nusq_array){ - nsq.int_struct=int_struct; //copy identical cross section data - nsq.InitializeInteractions(); //still need to initialize intermediate buffers - } - } - - if(evalThreads==1){ - unsigned int i = 0; - for(BaseSQUIDS& nsq : nusq_array){ - if(progressbar){ - std::cout << "Calculating cos(th) = " + std::to_string(costh_array[i]) << std::endl; - i++; - } - nsq.EvolveState(); - if(progressbar){ - std::cout << std::endl; - } - } - } - else{ - bool overrideProgressBar=progressbar; - Set_ProgressBar(false); - ThreadPool tpool(evalThreads); - std::vector> tasks; - tasks.reserve(nusq_array.size()); - for(BaseSQUIDS& nsq : nusq_array) - tasks.emplace_back(tpool.enqueue([&](){ nsq.EvolveState(); })); - unsigned int i = 0; - for(const auto& task : tasks){ - task.wait(); - if(overrideProgressBar) - std::cout << "cos(th) = " + std::to_string(costh_array[i++]) << " complete" << std::endl; - } - Set_ProgressBar(overrideProgressBar); - } - } - - /// \brief Returns the flavor composition at a given energy and zenith. - /// @param flv Neutrino flavor. - /// @param costh Cosine of the zenith. - /// @param enu Neutrino energy [eV]. - /// @param rho Index of the equation, see details. - /// \details When NeutrinoType is \c both \c rho specifies wether one - /// is considering neutrinos (0) or antineutrinos (1). Bilinear interpolation - /// is done in the logarithm of the energy and cos(zenith). - double EvalFlavor(unsigned int flv,double costh,double enu,unsigned int rho = 0, bool randomize_production_height = false) const { - // here the energy enters in eV - if(not iinistate) - throw std::runtime_error("nuSQUIDSAtm::Error::State not initialized."); - if(not inusquidsatm) - throw std::runtime_error("nuSQUIDSAtm::Error::nuSQUIDSAtm not initialized."); - if(not ( rho == 0 or rho == 1)) - throw std::runtime_error("nuSQUIDSAtm::Error::EvalFlavor rho has to be 0 or 1."); - - if( costh < *costh_array.begin() or costh > *costh_array.rbegin()) - throw std::runtime_error("nuSQUIDSAtm::Error::EvalFlavor::cos(th) out of bounds."); - if( enu < *enu_array.begin() or enu > *enu_array.rbegin() ) - throw std::runtime_error("nuSQUIDSAtm::Error::EvalFlavor::neutrino energy out of bounds.(Emin = " + - std::to_string(*enu_array.begin()) + - ",Emax = " + - std::to_string(*enu_array.rbegin()) + - ", Enu = " + std::to_string(enu) + ")"); - - auto cthit=std::lower_bound(costh_array.begin(),costh_array.end(),costh); - if(cthit==costh_array.end()) - throw std::runtime_error("SQUIDS::GetExpectationValueD : x value not in the array."); - if(cthit!=costh_array.begin()) - cthit--; - size_t cth_M=std::distance(costh_array.begin(),cthit); - - /*double logE = log(enu); - auto logeit=std::lower_bound(log_enu_array.begin(),log_enu_array.end(),logE); - if(logeit==log_enu_array.end()) - throw std::runtime_error("SQUIDS::GetExpectationValueD : x value not in the array."); - if(logeit!=log_enu_array.begin()) - logeit--; - size_t loge_M=std::distance(log_enu_array.begin(),logeit);*/ - - auto eit=std::lower_bound(enu_array.begin(),enu_array.end(),enu); - if(eit==enu_array.end()) - throw std::runtime_error("SQUIDS::GetExpectationValueD : x value not in the array."); - if(eit!=enu_array.begin()) - eit--; - size_t loge_M=std::distance(enu_array.begin(),eit); - - EarthAtm::Track track=earth_atm->MakeTrackWithCosine(costh); - double delta_t_final = track.GetFinalX()-track.GetInitialX(); - if (randomize_production_height){ - double production_height = gsl_ran_flat(r_gsl,-15*units.km,15*units.km); - delta_t_final += production_height; - } - - // assuming offsets are zero - double delta_t_1 = nusq_array[cth_M].Get_t() - nusq_array[cth_M].Get_t_initial(); - double delta_t_2 = nusq_array[cth_M+1].Get_t() - nusq_array[cth_M+1].Get_t_initial(); - double delta_t_final_1 = nusq_array[cth_M].GetTrackFast().GetFinalX() - nusq_array[cth_M].GetTrackFast().GetInitialX(); - double delta_t_final_2 = nusq_array[cth_M+1].GetTrackFast().GetFinalX() - nusq_array[cth_M+1].GetTrackFast().GetInitialX(); - double t_inter = 0.5*(delta_t_final*delta_t_1/delta_t_final_1 + delta_t_final*delta_t_2/delta_t_final_2); - - squids::SU_vector H0_at_enu = nusq_array[0].H0(enu,rho); - - struct storage_type{ - squids::SU_vector evol_proj, temp1; - storage_type(unsigned int dim):evol_proj(dim),temp1(dim){} - }; - #ifdef SQUIDS_THREAD_LOCAL - static SQUIDS_THREAD_LOCAL - #endif - storage_type storage(H0_at_enu.Dim()); - - storage.evol_proj = nusq_array[0].GetFlavorProj(flv,rho).Evolve(H0_at_enu,t_inter); - - //coefficients for energy interpolation - double f2=(enu-enu_array[loge_M])/(enu_array[loge_M+1]-enu_array[loge_M]); - double f1=1-f2; - - storage.temp1 =f1*nusq_array[cth_M ].GetState(loge_M ,rho); - storage.temp1+=f2*nusq_array[cth_M ].GetState(loge_M+1,rho); - double phiM=squids::SUTrace(storage.temp1,storage.evol_proj); - - storage.temp1 =f1*nusq_array[cth_M+1].GetState(loge_M ,rho); - storage.temp1+=f2*nusq_array[cth_M+1].GetState(loge_M+1,rho); - double phiP=squids::SUTrace(storage.temp1,storage.evol_proj); - - //perform angular interpolation - return(LinInter(costh,costh_array[cth_M],costh_array[cth_M+1],phiM,phiP)); - } - - /// \brief Returns the flavor composition at a given energy and zenith. - /// @param flv Neutrino flavor. - /// @param costh Cosine of the zenith. - /// @param enu Neutrino energy [eV]. - /// @param rho Index of the equation, see details. - /// @param scale Scale is a float that sets the averaged oscillation scaled. - /// @param avr Is a bool vector of size of the number of neutrinos that signal which scaled has been averaged. - /// \details When NeutrinoType is \c both \c rho specifies wether one - /// is considering neutrinos (0) or antineutrinos (1). Bilinear interpolation - /// is done in the logarithm of the energy and cos(zenith). - double EvalFlavor(unsigned int flv,double costh,double enu,unsigned int rho,double scale,std::vector avr) const { - // here the energy enters in eV - if(not iinistate) - throw std::runtime_error("nuSQUIDSAtm::Error::State not initialized."); - if(not inusquidsatm) - throw std::runtime_error("nuSQUIDSAtm::Error::nuSQUIDSAtm not initialized."); - if(not ( rho == 0 or rho == 1)) - throw std::runtime_error("nuSQUIDSAtm::Error::EvalFlavor rho has to be 0 or 1."); - - if( costh < *costh_array.begin() or costh > *costh_array.rbegin()) - throw std::runtime_error("nuSQUIDSAtm::Error::EvalFlavor::cos(th) out of bounds."); - if( enu < *enu_array.begin() or enu > *enu_array.rbegin() ) - throw std::runtime_error("nuSQUIDSAtm::Error::EvalFlavor::neutrino energy out of bounds.(Emin = " + - std::to_string(*enu_array.begin()) + - ",Emax = " + - std::to_string(*enu_array.rbegin()) + - ", Enu = " + std::to_string(enu) + ")"); - - auto cthit=std::lower_bound(costh_array.begin(),costh_array.end(),costh); - if(cthit==costh_array.end()) - throw std::runtime_error("SQUIDS::GetExpectationValueD : x value not in the array."); - if(cthit!=costh_array.begin()) - cthit--; - size_t cth_M=std::distance(costh_array.begin(),cthit); - - /*double logE = log(enu); - auto logeit=std::lower_bound(log_enu_array.begin(),log_enu_array.end(),logE); - if(logeit==log_enu_array.end()) - throw std::runtime_error("SQUIDS::GetExpectationValueD : x value not in the array."); - if(logeit!=log_enu_array.begin()) - logeit--; - size_t loge_M=std::distance(log_enu_array.begin(),logeit);*/ - - auto eit=std::lower_bound(enu_array.begin(),enu_array.end(),enu); - if(eit==enu_array.end()) - throw std::runtime_error("SQUIDS::GetExpectationValueD : x value not in the array."); - if(eit!=enu_array.begin()) - eit--; - size_t loge_M=std::distance(enu_array.begin(),eit); - - EarthAtm::Track track=earth_atm->MakeTrackWithCosine(costh); - double delta_t_final = track.GetFinalX()-track.GetInitialX(); - - // assuming offsets are zero - double delta_t_1 = nusq_array[cth_M].Get_t() - nusq_array[cth_M].Get_t_initial(); - double delta_t_2 = nusq_array[cth_M+1].Get_t() - nusq_array[cth_M+1].Get_t_initial(); - double delta_t_final_1 = nusq_array[cth_M].GetTrackFast().GetFinalX() - nusq_array[cth_M].GetTrackFast().GetInitialX(); - double delta_t_final_2 = nusq_array[cth_M+1].GetTrackFast().GetFinalX() - nusq_array[cth_M+1].GetTrackFast().GetInitialX(); - double t_inter = 0.5*(delta_t_final*delta_t_1/delta_t_final_1 + delta_t_final*delta_t_2/delta_t_final_2); - - squids::SU_vector H0_at_enu = nusq_array[0].H0(enu,rho); - - struct storage_type{ - squids::SU_vector evol_proj, temp1; - storage_type(unsigned int dim):evol_proj(dim),temp1(dim){} - }; - #ifdef SQUIDS_THREAD_LOCAL - static SQUIDS_THREAD_LOCAL - #endif - storage_type storage(H0_at_enu.Dim()); - // preevolution buffer - std::unique_ptr evol_buffer(new double[H0_at_enu.GetEvolveBufferSize()]); - H0_at_enu.PrepareEvolve(evol_buffer.get(),t_inter,scale,avr); - storage.evol_proj = nusq_array[0].GetFlavorProj(flv,rho).Evolve(evol_buffer.get()); - - //coefficients for energy interpolation - double f2=(enu-enu_array[loge_M])/(enu_array[loge_M+1]-enu_array[loge_M]); - double f1=1-f2; - - storage.temp1 =f1*nusq_array[cth_M ].GetState(loge_M ,rho); - storage.temp1+=f2*nusq_array[cth_M ].GetState(loge_M+1,rho); - - double phiM=squids::SUTrace(storage.temp1,storage.evol_proj); - - storage.temp1 =f1*nusq_array[cth_M+1].GetState(loge_M ,rho); - storage.temp1+=f2*nusq_array[cth_M+1].GetState(loge_M+1,rho); - double phiP=squids::SUTrace(storage.temp1,storage.evol_proj); - - //perform angular interpolation - return(LinInter(costh,costh_array[cth_M],costh_array[cth_M+1],phiM,phiP)); - } - - /// \brief Writes the object into an HDF5 file. - /// @param hdf5_filename Filename of the HDF5 into which save the object. - /// \details All contents are saved to the \c root of the HDF5 file. - /// @see ReadStateHDF5 - void WriteStateHDF5(std::string filename, bool overwrite = true) const{ - if(not iinistate) - throw std::runtime_error("nuSQUIDSAtm::Error::State not initialized."); - if(not inusquidsatm) - throw std::runtime_error("nuSQUIDSAtm::Error::nuSQUIDSAtm not initialized."); - - hid_t file_id,root_id; - hid_t dset_id; - // create HDF5 file - if(overwrite) - file_id = H5Fcreate (filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - else - file_id = H5Fcreate (filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT, H5P_DEFAULT); - if (file_id < 0) - throw std::runtime_error("nuSQUIDS::Error::Cannot create file at " + filename + "."); - root_id = H5Gopen(file_id, "/",H5P_DEFAULT); - - // write the zenith range - hsize_t costhdims[1]={costh_array.extent(0)}; - dset_id = H5LTmake_dataset(root_id,"zenith_angles",1,costhdims,H5T_NATIVE_DOUBLE,costh_array.get_data()); - hsize_t energydims[1]={enu_array.extent(0)}; - dset_id = H5LTmake_dataset(root_id,"energy_range",1,energydims,H5T_NATIVE_DOUBLE,enu_array.get_data()); - - H5Gclose (root_id); - H5Fclose (file_id); - - unsigned int i = 0; - for(const BaseSQUIDS& nsq : nusq_array){ - // use only the first one to write the cross sections on /crosssections - // the last false is because the HDF5 file has already been created, should not overwrite - nsq.WriteStateHDF5(filename,"/costh_"+std::to_string(costh_array[i]),i==0,"crosssections",false); - i++; - } - } - - /// \brief Reads the object from an HDF5 file. - /// @param hdf5_filename Filename of the HDF5 to use for construction. - /// \details All contents are assumed to be saved to the \c root of the HDF5 file. - /// @see WriteStateHDF5 - void ReadStateHDF5(std::string hdf5_filename){ - hid_t file_id,group_id,root_id; - // open HDF5 file - H5File file(H5Fopen(hdf5_filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)); - if (file < 0) - throw std::runtime_error("nuSQUIDSAtm::ReadStateHDF5: Unable to open file: " + hdf5_filename + ". No such file or directory"); - root_id = H5Gopen(file, "/",H5P_DEFAULT); - group_id = root_id; - - // read the zenith range dimension - hsize_t costhdims[1]; - H5LTget_dataset_info(group_id, "zenith_angles", costhdims, nullptr, nullptr); - - double data[costhdims[0]]; - H5LTread_dataset_double(group_id, "zenith_angles", data); - costh_array.resize(std::vector {static_cast(costhdims[0])}); - for (unsigned int i = 0; i < costhdims[0]; i ++) - costh_array[i] = data[i]; - - hsize_t energydims[1]; - H5LTget_dataset_info(group_id, "energy_range", energydims, nullptr, nullptr); - - double enu_data[energydims[0]]; - H5LTread_dataset_double(group_id, "energy_range", enu_data); - enu_array.resize(std::vector{static_cast(energydims[0])});log_enu_array.resize(std::vector{static_cast(energydims[0])}); - for (unsigned int i = 0; i < energydims[0]; i ++){ - enu_array[i] = enu_data[i]; - log_enu_array[i] = log(enu_data[i]); - } - - H5Gclose(root_id); - - // resize apropiately the nuSQUIDSAtm container vector - nusq_array.clear(); - nusq_array = std::vector(costhdims[0]); - - unsigned int i = 0; - for(BaseSQUIDS& nsq : nusq_array){ - if(i==0){ - // read the cross sections stored in /crosssections - nsq.ReadStateHDF5(hdf5_filename,"/costh_"+std::to_string(costh_array[i]),"crosssections"); - if(nsq.iinteraction) - int_struct = nsq.GetInteractionStructure(); - } else { - // re-use the shared cross sections - nsq.ReadStateHDF5Internal(hdf5_filename,"/costh_"+std::to_string(costh_array[i]),int_struct); - } - i++; - } - earth_atm = std::dynamic_pointer_cast(nusq_array.front().GetBody()); - assert(earth_atm); - - iinistate = true; - inusquidsatm = true; - } - - /// \brief Sets the mixing parameters to default. - void Set_MixingParametersToDefault(){ - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_MixingParametersToDefault(); - } - } - - /// \brief Sets the mixing angle th_ij. - /// @param i the (zero-based) index of the first state - /// @param j the (zero-based) index of the second state - /// must be larger than \c i. - /// @param angle Angle to use in radians. - /// \details Sets the neutrino mixing angle. In our zero-based convention, e.g., the th_12 is i = 0, j = 1.,etc. - void Set_MixingAngle(unsigned int i, unsigned int j,double angle){ - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_MixingAngle(i,j,angle); - } - } - - /// \brief Returns the mixing angle th_ij in radians. - /// @param i the (zero-based) index of the first state - /// @param j the (zero-based) index of the second state - /// must be larger than \c i. - /// \details Gets the neutrino mixing angle. In our zero-based convention, e.g., the th_12 is i = 0, j = 1.,etc. - double Get_MixingAngle(unsigned int i, unsigned int j) const{ - return nusq_array[0].Get_MixingAngle(i,j); - } - - /// \brief Sets the CP phase for the ij-rotation. - /// @param i the (zero-based) index of the first state - /// @param j the (zero-based) index of the second state - /// must be larger than \c i. - /// @param angle Phase to use in radians. - /// \details Sets the CP phase for the ij-rotation. In our zero-based convention, e.g., the delta_13 = delta_CP is i = 0, j = 2.,etc. - void Set_CPPhase(unsigned int i, unsigned int j,double angle){ - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_CPPhase(i,j,angle); - } - } - - /// \brief Returns the CP phase of the ij-rotation in radians. - /// @param i the (zero-based) index of the first state - /// @param j the (zero-based) index of the second state - /// must be larger than \c i. - /// \details Gets the CP phase for the ij-rotation. In our zero-based convention, e.g., the delta_13 = delta_CP is i = 0, j = 2.,etc. - double Get_CPPhase(unsigned int i, unsigned int j) const{ - return nusq_array[0].Get_CPPhase(i,j); - } - - /// \brief Sets the square mass difference with respect to first mass eigenstate - /// @param i the (zero-based) index of the second state - /// must be larger than \c 0. - /// @param sq Square mass difference in eV^2. - /// \details Sets square mass difference with respect to the first mass eigenstate. In our zero-based convention, e.g., the \f$\Delta m^2_{12}\f$ corresponds to (i = 1),etc. - void Set_SquareMassDifference(unsigned int i,double sq){ - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_SquareMassDifference(i,sq); - } - } - - /// \brief Returns the square mass difference between state-i and the first mass eigenstate. - /// @param i the (zero-based) index of the first state - /// must be larger than \c i. - /// \details Returns square mass difference with respect to the first mass eigenstate. In our zero-based convention, e.g., the \f$\Delta m^2_{12}\f$ corresponds to (i = 1),etc. - double Get_SquareMassDifference(unsigned int i) const{ - return nusq_array[0].Get_SquareMassDifference(i); - } - - /// \brief Set the initial integration step size - /// \param h initial step size - void Set_h(double h){ - for(BaseSQUIDS& nsq : nusq_array) - nsq.Set_h(h); - } - - /// \brief Set the initial integration step size for a single energy node - /// \param h initial step size - /// \param idx energy node index - void Set_h(double h, unsigned int idx){ - nusq_array[idx].Set_h(h); - } - - /// \brief Set the maximum integration step size - /// \param h maximum step size - void Set_h_max(double h){ - for(BaseSQUIDS& nsq : nusq_array) - nsq.Set_h_max(h); - } - - /// \brief Set the maximum integration step size for a single energy node - /// \param h maximum step size - /// \param idx energy node index - void Set_h_max(double h, unsigned int idx){ - nusq_array[idx].Set_h_max(h); - } - - /// \brief Set the minimum integration step size - /// \param h minimum step size - void Set_h_min(double h){ - for(BaseSQUIDS& nsq : nusq_array) - nsq.Set_h_min(h); - } - - /// \brief Set the minimum integration step size for a single energy node - /// \param h minimum step size - /// \param idx energy node index - void Set_h_min(double h, unsigned int idx){ - nusq_array[idx].Set_h_min(h); - } - - /// \brief Set the allowed absolute numerical error - /// \param eps maximum allowed error - void Set_abs_error(double eps){ - for(BaseSQUIDS& nsq : nusq_array) - nsq.Set_abs_error(eps); - } - - /// \brief Set the allowed absolute numerical error for a single energy node - /// \param eps maximum allowed error - /// \param idx energy node index - void Set_abs_error(double eps, unsigned int idx){ - nusq_array[idx].Set_abs_error(eps); - } - - /// \brief Set the allowed relative numerical error - /// \param eps maximum allowed error - void Set_rel_error(double eps){ - for(BaseSQUIDS& nsq : nusq_array) - nsq.Set_rel_error(eps); - } - - /// \brief Set the allowed relative numerical error for a single energy node - /// \param eps maximum allowed error - /// \param idx energy node index - void Set_rel_error(double eps, unsigned int idx){ - nusq_array[idx].Set_rel_error(eps); - } - - /// \brief Sets the GSL solver - /// @param opt GSL stepper function. - void Set_GSL_step(gsl_odeiv2_step_type const * opt){ - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_GSL_step(opt); - } - } - - /// \brief Toggles the progress bar printing on and off - /// @param opt If \c true a progress bar will be printed. - void Set_ProgressBar(bool opt){ - progressbar = opt; - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_ProgressBar(opt); - } - } - - /// \brief Returns the interaction structure. - std::shared_ptr GetInteractionStructure() const { - return int_struct; - } - - /// \brief Returns number of energy nodes. - size_t GetNumE() const{ - return enu_array.extent(0); - } - - /// \brief Returns number of zenith nodes. - size_t GetNumCos() const{ - return costh_array.extent(0); - } - - /// \brief Returns the number of neutrino flavors. - unsigned int GetNumNeu() const{ - return nusq_array[0].GetNumNeu(); - } - - /// \brief Returns the number of rho equations. - unsigned int GetNumRho() const{ - return nusq_array[0].GetNumRho(); - } - - /// \brief Returns number of nodes. - size_t GetNumNodes() const{ - return nusq_array.size(); - } - - /// \brief Returns the (evolved) states of all nodes - /// @param rho Index of equation. In `both` propagation mode, neutrinos = 0 and - /// antineutrinos = 1. - marray GetStates(unsigned int rho = 0){ - marray states {GetNumNodes(),GetNumNeu()*GetNumNeu()}; - for(int i=0; i state_vec = nusq_array[i].GetState(0, rho).GetComponents(); - for(int j=0; j GetERange() const{ - return enu_array; - } - - /// \brief Returns the cos(zenith) nodes values. - marray GetCosthRange() const{ - return costh_array; - } - - /// \brief Returns the nuSQUIDSBase object for ci-th zenith. - BaseSQUIDS& GetnuSQuIDS(unsigned int ci) { - if(ci >= nusq_array.size()) - throw std::runtime_error("nuSQUIDSAtm::GetnuSQuIDS: input index out of range."); - return nusq_array[ci]; - } - - /// \brief Returns the vector of nuSQUIDSBase objects. - std::vector& GetnuSQuIDS() { - return nusq_array; - } - - /// \brief Toggles tau regeneration on and off. - /// @param opt If \c true tau regeneration will be considered. - void Set_TauRegeneration(bool opt){ - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_TauRegeneration(opt); - } - } - - void Set_GlashowResonance(bool opt){ - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_GlashowResonance(opt); - } - } - - /// \brief Toggles neutrino oscillations on and off. - /// @param opt If \c true neutrino oscillations will be considered. - void Set_IncludeOscillations(bool opt){ - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_IncludeOscillations(opt); - } - } - - /// \brief Toggles if non-hard processes constant density fast evolution will be used. - /// \param opt If \c true fast evolution will be attempted on constatnt density. - void Set_AllowConstantDensityOscillationOnlyEvolution(bool opt){ - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_AllowConstantDensityOscillationOnlyEvolution(opt); - } - } - - - /// \brief Toggles positivization of the flux. - /// @param opt If \c true the flux will be forced to be positive every \c positivization_step. - void Set_PositivityConstrain(bool opt){ - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_PositivityConstrain(opt); - } - } - /// \brief Stes the step upon which the positivity correction would be apply. - /// @param step The step upon which the positivization will take place. - void Set_PositivityConstrainStep(double step){ - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_PositivityConstrainStep(step); - } - } - - /// \brief Set the number of threads used for evolution. - /// - /// \param nThreads the number of execution threads EvolveState should use. - /// If zero, use an automatic, system defined number. - void Set_EvalThreads(unsigned int nThreads){ - evalThreads=nThreads; - if(evalThreads==0){ //if the user wants automatic selection - evalThreads=std::thread::hardware_concurrency(); - if(evalThreads==0) //if autodetection failed - evalThreads=1; - } - } - - /// \brief Get the number of threads used for evolution. - /// - /// \return The number of execution threads EvolveState should use. - unsigned int Get_EvalThreads() const{ - return(evalThreads); - } - - /// \brief Stes Earth object to be use. - /// @param earth Shared pointer to Earth object. - void Set_EarthModel(std::shared_ptr earth){ - earth_atm=earth; - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_Body(earth); - } - } - - /// \brief Returns the neutrino interaction cross sections - std::shared_ptr GetNeutrinoCrossSections(){ - return(ncs); - } - - /// \brief Sets the neutrino interaction cross sections - void SetNeutrinoCrossSections(std::shared_ptr xs){ - ncs=xs; - for(BaseSQUIDS& nsq : nusq_array) - nsq.SetNeutrinoCrossSections(ncs); - } - - void Set_EvolLowPassCutoff(double val){ - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_EvolLowPassCutoff(val); - } - } - - void Set_EvolLowPassScale(double val){ - for(BaseSQUIDS& nsq : nusq_array){ - nsq.Set_EvolLowPassScale(val); - } - } - - void Set_NeutrinoSources(bool opt){ - for(BaseSQUIDS& nsq : nusq_array) - nsq.Set_NeutrinoSources(opt); - } -}; - - -} // close namespace -#endif From 7beb7c2ee5994bc936e92481169f8b46755ba20b Mon Sep 17 00:00:00 2001 From: ConnorSponsler <33155897+ConnorSponsler@users.noreply.github.com> Date: Sun, 19 Nov 2023 16:41:29 -0500 Subject: [PATCH 10/12] Delete duplicate body.h --- body.h | 873 --------------------------------------------------------- 1 file changed, 873 deletions(-) delete mode 100644 body.h diff --git a/body.h b/body.h deleted file mode 100644 index c35d7aeb..00000000 --- a/body.h +++ /dev/null @@ -1,873 +0,0 @@ - /****************************************************************************** - * 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 * - ******************************************************************************/ - - -#ifndef NUSQUIDS_BODY_H -#define NUSQUIDS_BODY_H - -#if __cplusplus < 201103L -#error C++11 compiler required. Update your compiler and use the flag -std=c++11 -#endif - -#include -#include -#include -#include -#include -#include -#include "nuSQuIDS/marray.h" -#include "nuSQuIDS/tools.h" - -namespace nusquids{ - -// Forward declaration of nuSQuIDS class -class nuSQUIDS; - -/// \class Body -/// \brief Abstract body class. -/// \details This abstract class serves as a prototype -/// of neutrino propagation environments. When implementing -/// the subclasses one must define the density() and ye() -class Body{ - protected: - /// \brief Body parameters. - std::vector BodyParams; - /// \brief bool that signal if the body is constant density or not - bool is_constant_density = false; - public: - /// \brief Body Constructor; - Body(){} - virtual ~Body(){} - - /// \brief Serialization function - virtual void Serialize(hid_t group) const=0; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - - /// \class Track - /// \brief Trajectory subclass. Specifies the trajectory - /// to take within the body. - class Track{ - friend class Body; - protected: - /// \brief Current position. - double x; - /// \brief Initial position. - double xini; - /// \brief Final position. - double xend; - public: - /// \brief Default constructor. - Track(double x,double xini, double xend): x(x), xini(xini),xend(xend) {} - /// \brief Default constructor. - Track(double xini, double xend): Track(xini,xini,xend) {} - /// \brief Serialization function - virtual void Serialize(hid_t group) const=0; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - /// \brief Get track object name - static std::string GetName() {return "BodyTrack";}; - /// \brief Sets the current position along the trajectory. - void SetX(double y){ x = y; } - /// \brief Returns the current position along the trajectory. - double GetX() const { return x; } - /// \brief Returns the initial position measured along the trajectory. - double GetInitialX() const { return xini; } - /// \brief Returns the final position measured along the trajectory. - double GetFinalX() const { return xend; } - /// \brief Reverses the track initial and final position - void ReverseTrack() { - double xtmp = xend; - xend = xini; - xini = xtmp; - } - /// \brief Returns parameters that define the trajectory. - std::vector GetTrackParams() const { - std::vector TrackParams{xini,xend}; - FillDerivedParams(TrackParams); - return TrackParams; - } - /// Should be implemented by derived classes to append their - /// additional parameters to TrackParams - virtual void FillDerivedParams(std::vector& TrackParams) const{} - }; - /// \brief Return the density at a given trajectory object. - virtual double density(const Track&) const {return 0.0;} - /// \brief Return the electron fraction at a given trajectory object. - virtual double ye(const Track&) const {return 1.0;} - /// \brief Returns parameters that define the body. - const std::vector& GetBodyParams() const { return BodyParams;} - /// \brief Returns the body identifier. - static unsigned int GetId() {return 0;} - /// \brief Returns the name of the body. - static std::string GetName() {return "Body";} - /// \brief Return true if the body is a constant density. - virtual bool IsConstantDensity() const {return is_constant_density;} - /// \brief Return true if the body is a constant density. - virtual void SetIsConstantDensity(bool icd) {is_constant_density = icd;} - /// \brief Sets the injected neutrino flux from the Body at a given Track location. - /// \details This function is called by nuSQuIDS during the calculation to obtain the injected flux. - /// @param flux Is a three dimensional array with dimensions: energy, rho, and flavor, which has previously been allocated by nuSQuIDS. - /// @param Track Trajectory object. - /// @param nuSQuIDS nuSQuIDS object that queries the flux. - virtual void injected_neutrino_flux(marray& flux, const Track&, const nuSQUIDS&) {} -}; - -// type defining -typedef Body::Track GenericTrack; - -/// \class Vacuum -/// \brief Vacuum body. -/// \details A zero density body. -class Vacuum: public Body { - public: - /// \brief Constructor. - Vacuum():Body(){} - - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - - /// \class Track - /// \brief Vacuum trajectory - class Track: public Body::Track { - public : - /// \brief Rectilinear trajectory from xini to xend. - /// @param x current position [eV^-1]. - /// @param xini Initial position [eV^-1]. - /// @param xend Final position [eV^-1]. - Track(double x,double xini,double xend):Body::Track(x,xini,xend){} - /// \brief Rectilinear trajectory from xini to xend. - /// @param xini Initial position [eV^-1]. - /// @param xend Final position [eV^-1]. - Track(double xini,double xend):Track(xini,xini,xend){} - /// \brief Construct a trajectory of distance xend. - /// @param xend Final position in eV^-1. - /// \details In this case initial position is assumed 0. - Track(double xend):Track(0.0,xend){} - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - /// \brief Get track object name - static std::string GetName() {return "VacuumTrack";} - }; - /// \brief Returns the body identifier. - static unsigned int GetId() {return 1;} - /// \brief Returns the name of the body. - static std::string GetName() {return "Vacuum";} - /// \brief Returns the density in g/cm^3 - double density(const GenericTrack&) const; - /// \brief Returns the electron fraction - double ye(const GenericTrack&) const; - /// \brief Returns true as this body has constant density by definition - bool IsConstantDensity() const; -}; - -/// \class ConstantDensity -/// \brief Constant density body. -/// \details A body with constant density -class ConstantDensity: public Body{ - private: - /// \brief Constant density in g/cm^3 - const double constant_density; - /// \brief Constant electron fraction - const double constant_ye; - public: - /// \brief Constructor - /// @param density Density in g/cm^3. - /// @param ye electron fraction. - ConstantDensity(double density,double ye); - - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - - /// \class Track - /// \brief Constant density trajectory - class Track: public Body::Track{ - public : - /// \brief Rectilinear trajectory from xini to xend. - /// @param x current position [eV^-1]. - /// @param xini Initial position [eV^-1]. - /// @param xend Final position [eV^-1]. - Track(double x,double xini,double xend):Body::Track(x,xini,xend){} - /// \brief Construct a trajectory between an initial position and final position. - /// @param xini Initial position in eV^-1. - /// @param xend Final position in eV^-1. - Track(double xini,double xend):Track(xini,xini,xend){} - /// \brief Construct a trajectory of distance xend. - /// @param xend Final position in eV^-1. - /// \details In this case initial position is assumed 0. - Track(double xend):Track(0.0,xend){} - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - /// \brief Get track object name - static std::string GetName() {return "ConstantDensityTrack";} - }; - - /// \brief Returns the body identifier. - static unsigned int GetId() {return 2;} - /// \brief Returns the name of the body. - static std::string GetName() {return "ConstantDensity";} - /// \brief Returns the density in g/cm^3 - double density(const GenericTrack&) const; - /// \brief Returns the electron fraction - double ye(const GenericTrack&) const; - /// \brief Returns true as this body has constant density by definition - bool IsConstantDensity() const; -}; - -/// \class VariableDensity -/// \brief Variable user defined density body. -/// \details The user provides arrays which contain the body -/// information. -class VariableDensity: public Body{ - private: - /// \brief Array of length \c arraysize containing spline position nodes. - std::vector x_arr; - /// \brief Array of length \c arraysize containing the density at each position nodes. - std::vector density_arr; - /// \brief Array of length \c arraysize containing the electron fraction at each position nodes. - std::vector ye_arr; - /// \brief Size of array used to create spline. - unsigned int arraysize; - - /// \brief Minimum value of \c x_arr. - double x_max; - /// \brief Maximum value of \c x_arr. - double x_min; - - /// \brief Density spline - AkimaSpline inter_density; - /// \brief Electron fraction spline - AkimaSpline inter_ye; - public: - /// \brief Constructor. - /// @param x Vector containing position nodes in cm. - /// @param density Density, in g/cm^3, at each of the nodes. - /// @param ye Electron fraction at each of the nodes. - /// \pre All input vectors must be of equal size. - VariableDensity(std::vector x,std::vector density,std::vector ye); - /// \brief Destructor. - ~VariableDensity(); - - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - - /// \class Track - /// \brief Variable density trajectory - class Track: public Body::Track{ - public : - /// \brief Rectilinear trajectory from xini to xend. - /// @param x current position [eV^-1]. - /// @param xini Initial position [eV^-1]. - /// @param xend Final position [eV^-1]. - Track(double x,double xini,double xend):Body::Track(x,xini,xend){} - /// \brief Construct a trajectory between an initial position and final position. - /// @param xini Initial position in eV^-1. - /// @param xend Final position in eV^-1. - Track(double xini,double xend):Track(xini,xini,xend){} - /// \brief Construct a trajectory of distance xend. - /// @param xend Final position in eV^-1. - /// \details In this case initial position is assumed 0. - Track(double xend):Track(0.0,xend){} - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - /// \brief Get track object name - static std::string GetName() {return "VariableDensityTrack";} - - }; - /// \brief Returns the body identifier. - static unsigned int GetId() {return 3;} - /// \brief Returns the name of the body. - static std::string GetName() {return "VariableDensity";} - /// \brief Returns the density in g/cm^3 - double density(const GenericTrack&) const; - /// \brief Returns the electron fraction - double ye(const GenericTrack&) const; -}; - -/// \class Earth -/// \brief Earth model based on PREM. -class Earth: public Body{ - private: - /// \brief Radius of the Earth. - double radius; - /// \brief Earth radius position array - std::vector earth_radius; - /// \brief Earth density array - std::vector earth_density; - /// \brief Earth electron fraction array - std::vector earth_ye; - /// \brief Data arrays size - unsigned int arraysize; - - /// \brief Density spline - AkimaSpline inter_density; - /// \brief Electron fraction spline - AkimaSpline inter_ye; - - /// \brief Minimum radius. - double x_radius_min; - /// \brief Maximum radius. - double x_radius_max; - /// \brief Density at minimum radius. - double x_rho_min; - /// \brief Density at maximum radius. - double x_rho_max; - /// \brief Electron fraction at minimum radius. - double x_ye_min; - /// \brief Electron fraction at maximum radius. - double x_ye_max; - public: - /// \brief Default constructor using supplied PREM. - Earth(); - /// \brief Constructor from a user supplied Earth model. - /// @param earthmodel Path to the Earth model file. - /// \details The input file should have three columns. - /// The first one must run from zero to one representing - /// the center and surface of the Earth respectively. The - /// second column must contain the Earth density in g/cm^3 at - /// a given position, while the third column must contain - /// the electron fraction. - Earth(std::string earthmodel); - /// \brief Constructor in which the user provides, as vectors, the - /// Earth properties. - /// @param x Vector containing position nodes in cm. - /// @param rho Density, in g/cm^3, at each of the nodes. - /// @param ye Electron fraction at each of the nodes. - /// \pre All input vectors must be of equal size. - Earth(std::vector x,std::vector rho,std::vector ye); - /// \brief Destructor. - ~Earth(); - - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - - /// \class Track - /// \brief Earth trajectory - class Track: public Body::Track{ - private: - /// \brief Baseline of the neutrino experiment in natural units. - const double baseline; - public : - /// \brief Rectilinear trajectory from xini to xend. - /// @param x current position [eV^-1]. - /// @param xini Initial position [eV^-1]. - /// @param xend Final position [eV^-1]. - Track(double x,double xini,double xend,double baseline):Body::Track(x,xini,xend),baseline(baseline){} - /// \brief Construct a trajectory between an initial position and final position. - /// @param xini Initial position in eV^-1. - /// @param xend Final position in eV^-1. - /// @param baseline Baseline of experiment in eV^-1. - Track(double xini,double xend,double baseline):Track(xini,xini,xend,baseline){} - /// \brief Construct a trajectory where the neutrino travels a distance given by baseline. - /// @param baseline Traverse distance in eV^-1. - /// \details In this case \c xini = 0, \c xend = \c baseline. - Track(double baseline):Track(0.,baseline,baseline){} - /// \brief Returns the neutrino baseline in natural units. - double GetBaseline() const {return baseline;} - virtual void FillDerivedParams(std::vector& TrackParams) const; - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - /// \brief Get track object name - static std::string GetName() {return "EarthTrack";} - }; - /// \brief Returns the body identifier. - static unsigned int GetId() {return 4;} - - /// \brief Returns the name of the body. - static std::string GetName() {return "Earth";} - - /// \brief Returns the density in g/cm^3 - double density(const GenericTrack&) const; - /// \brief Returns the electron fraction - double ye(const GenericTrack&) const; - - /// \brief Returns the radius of the Earth in natural units. - double GetRadius() const {return radius;} -}; - -/// \class Sun -/// \brief A model of the Sun. -class Sun: public Body{ - private: - /// \brief Array that contains all the Solar model parameters. - marray sun_model; - /// \brief Array that contains all the Solar electron density parameters. - marray sun_model_nele; - /// \brief Array of length \c arraysize containing spline position nodes. - std::vector sun_radius; - /// \brief Array of length \c arraysize containing the density at position nodes. - std::vector sun_density; - /// \brief Array of length \c arraysize containing the hydrogen fraction at position nodes. - std::vector sun_xh; - - /// \brief Size of \c sun_radius array. - unsigned int arraysize; - - /// \brief Radius of the Sun. - double radius; - - /// \brief Density spline - AkimaSpline inter_density; - /// \brief Hydrogen fraction spline - AkimaSpline inter_xh; - - /// \brief Returns the density in g/cm^3 at a given radius fraction x - /// @param x Radius fraction: 0:center, 1:surface. - double rdensity(double x) const; - /// \brief Returns the electron fraction at a given radius fraction x - /// @param x Radius fraction: 0:center, 1:surface. - double rxh(double x) const; - public: - /// \brief Detault constructor. - Sun(); - /// \brief Constructor in which the user provides, as vectors, the - /// Sun properties. - /// @param x Vector containing position nodes in cm. - /// @param rho Density, in g/cm^3, at each of the nodes. - /// @param xh Hydrogen fraction at each point. - /// \pre All input vectors must be of equal size. - Sun(std::vector x,std::vector rho,std::vector xh); - /// \brief Constructor from a user supplied solar model. - /// @param sunmodel Path to the Sun model file. - /// \details The input file should have the same columns as John Bahcalls solar model file. - /// The second column one must run from zero to one representing - /// the center and surface of the Sun respectively. The - /// fourth column must contain the Sun density in g/cm^3 at - /// a given position, while the seven column must contain - /// the hydrogen fraction which is related to the electron fraction by ye = 0.5*(1.0+rxh(r)). - /// Internally we assume that the solar radius is 695980.0 kilometers. - Sun(std::string sunmodel); - - - /// \brief Destructor - ~Sun(); - - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - - /// \class Track - /// \brief Sun trajectory - class Track: public Body::Track{ - public : - /// \brief Rectilinear trajectory from xini to xend. - /// @param x current position [eV^-1]. - /// @param xini Initial position [eV^-1]. - /// @param xend Final position [eV^-1]. - Track(double x,double xini,double xend):Body::Track(x,xini,xend){} - /// \brief Construct a trajectory between an initial position and final position. - /// @param xini Initial position in eV^-1. - /// @param xend Final position in eV^-1. - /// \details The trajectory is measured from the sun center which is set to zero. - Track(double xini,double xend):Track(xini,xini,xend){} - /// \brief Construct a trajectory from the sun center to a final position. - /// @param xend Final position in eV^-1. - /// \details The trajectory is measured from the sun center which is set to zero. - Track(double xend):Track(0.,xend){} - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - /// \brief Get track object name - static std::string GetName() {return "SunTrack";} - }; - - /// \brief Returns the body identifier. - static unsigned int GetId() {return 5;} - - /// \brief Returns the name of the body. - static std::string GetName() {return "Sun";} - - /// \brief Returns the density in g/cm^3 - double density(const GenericTrack&) const; - /// \brief Returns the electron fraction - double ye(const GenericTrack&) const; - - /// \brief Returns the radius of the Sun in natural units. - double GetRadius() const {return radius;} -}; - -/// \class SunASnu -/// \brief A model of the Sun with atmospheric solar neutrinos geometry. -class SunASnu: public Body{ - private: - /// \brief Array that contains all the Solar model parameters. - marray sun_model; - /// \brief Array of length \c arraysize containing spline position nodes. - std::vector sun_radius; - /// \brief Array of length \c arraysize containing the density at position nodes. - std::vector sun_density; - /// \brief Array of length \c arraysize containing the hydrogen fraction at position nodes. - std::vector sun_xh; - - /// \brief Size of \c sun_radius array. - unsigned int arraysize; - - /// \brief Radius of the Sun. - double radius; - - /// \brief Density spline - AkimaSpline inter_density; - /// \brief Hydrogen fraction spline - AkimaSpline inter_xh; - - /// \brief Returns the density in g/cm^3 at a given radius fraction x - /// @param x Radius fraction: 0:center, 1:surface. - double rdensity(double x) const; - /// \brief Returns the electron fraction at a given radius fraction x - /// @param x Radius fraction: 0:center, 1:surface. - double rxh(double x) const; - public: - /// \brief Detault constructor. - SunASnu(); - /// \brief Constructor from a user supplied solar model. - /// @param sunmodel Path to the Sun model file. - /// \details The input file should have the same columns as John Bahcalls solar model file. - /// The second column one must run from zero to one representing - /// the center and surface of the Sun respectively. The - /// fourth column must contain the Sun density in g/cm^3 at - /// a given position, while the seven column must contain - /// the hydrogen fraction which is related to the electron fraction by ye = 0.5*(1.0+rxh(r)). - /// Internally we assume that the solar radius is 695980.0 kilometers. - SunASnu(std::string sunmodel); - /// \brief Constructor in which the user provides, as vectors, the - /// Sun properties. - /// @param x Vector containing position nodes in cm. - /// @param rho Density, in g/cm^3, at each of the nodes. - /// @param xh Hydrogen fraction at each point. - /// \pre All input vectors must be of equal size. - SunASnu(std::vector x,std::vector rho,std::vector xh); - - ~SunASnu(); - - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - - /// \class Track - /// \brief SunASnu trajectory - class Track: public Body::Track{ - friend class SunASnu; - private: - /// \brief Radius of the Sun. - double radius_nu; - /// \brief Impact parameter. - double b_impact; - public: - /// \brief Rectilinear trajectory from xini to xend. - /// @param x current position [eV^-1]. - /// @param xini Initial position [eV^-1]. - /// @param b_impact impact parameter in eV^-1. - Track(double x,double xini,double b_impact); - /// \brief Construct a trajectory between an initial position and final position. - /// @param xini Initial position in eV^-1. - /// @param b_impact impact parameter in eV^-1. - /// \details The trajectory baseline is determined by the impact parameter and starts - /// at \c xini, and ends when the neutrino exits the sun. - Track(double xini,double b_impact):Track(xini,xini,b_impact){} - /// \brief Construct a for a given impact parameter. - /// @param b_impact_ impact parameter in eV^-1. - /// \details The trajectory baseline is determined by the impact parameter and starts - /// at \c xini = 0, and ends when the neutrino exits the sun. - Track(double b_impact_):Track(0.0,b_impact_){} - virtual void FillDerivedParams(std::vector& TrackParams) const; - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - /// \brief Get track object name - static std::string GetName() {return "SunASnuTrack";} - }; - - /// \brief Returns the body identifier. - static unsigned int GetId() {return 6;} - - /// \brief Returns the name of the body. - static std::string GetName() {return "SunASnu";} - - /// \brief Returns the density in g/cm^3 - double density(const GenericTrack&) const; - /// \brief Returns the electron fraction - double ye(const GenericTrack&) const; - - /// \brief Returns the radius of the Sun in natural units. - double GetRadius() const {return radius;} -}; - -/// \class EarthAtm -/// \brief A model of the Earth with atmospheric neutrinos geometry. -class EarthAtm: public Body{ - protected: - /// \brief Radius of the Earth. - double radius; - /// \brief Height of the atmosphere. - double atm_height; - /// \brief Radius of the Earth plus atmosphere. - double earth_with_atm_radius; - - /// \brief Earth radius position array - std::vector earth_radius; - /// \brief Earth density array - std::vector earth_density; - /// \brief Earth electron fraction array - std::vector earth_ye; - /// \brief Data arrays size - unsigned int arraysize; - - /// \brief Density spline - AkimaSpline inter_density; - /// \brief Electron fraction spline - AkimaSpline inter_ye; - - /// \brief Minimum radius. - double x_radius_min; - /// \brief Maximum radius. - double x_radius_max; - /// \brief Density at minimum radius. - double x_rho_min; - /// \brief Density at maximum radius. - double x_rho_max; - /// \brief Electron fraction at minimum radius. - double x_ye_min; - /// \brief Electron fraction at maximum radius. - double x_ye_max; - - public: - /// \brief Default constructor using supplied PREM. - EarthAtm(); - /// \brief Constructor from a user supplied Earth model. - /// @param earthmodel Path to the Earth model file. - /// \details The input file should have three columns. - /// The first one must run from zero to one representing - /// the center and surface of the Earth respectively. The - /// second column must contain the Earth density in g/cm^3 at - /// a given position, while the third column must contain - /// the electron fraction. - EarthAtm(std::string earthmodel); - /// \brief Constructor in which the user provides, as vectors, the - /// Earth properties. - /// @param x Vector containing position nodes in cm. - /// @param rho Density, in g/cm^3, at each of the nodes. - /// @param ye Electron fraction at each of the nodes. - /// \pre All input vectors must be of equal size. - EarthAtm(std::vector x,std::vector rho,std::vector ye); - ~EarthAtm(); - - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - - /// \class Track - /// \brief EarthAtm trajectory - class Track: public Body::Track{ - friend class EarthAtm; - friend class EmittingEarthAtm; - private: - /// \brief Cosine of the zenith angle. - double cosphi; - /// \brief Radius of the Earth. - double earth_radius; - /// \brief Height of the atmosphere. - double atmheight; - /// \brief Baseline. - double L; - - ///Private constructor which does not initialize all members. - ///Intended for use only by makeWithCosine. - Track(); - public : - /// \brief Construct trajectory. - /// @param x_ current position [eV^-1]. - /// @param phi Zenith angle in radians. - Track(double x_,double phi,double earth_radius_,double atmheight_): - Track(phi,earth_radius_,atmheight_){x=x_; assert(x >= 0);}; - /// \brief Construct trajectory. - /// @param phi Zenith angle in radians. - Track(double phi,double earth_radius_,double atmheight_); - /// \brief Returns the neutrino baseline in natural units. - double GetBaseline() const {return L;} - virtual void FillDerivedParams(std::vector& TrackParams) const; - ///Construct a track with the cosine of the zenith angle - static Track MakeWithCosine(double cosphi,double earth_radius_,double atmheight_); - /// \brief Serialization function - void Serialize(hid_t group) const; - /// \brief Deserialization function - static std::shared_ptr Deserialize(hid_t group); - /// \brief Get track object name - static std::string GetName() {return "EarthAtmTrack";} - }; - /// \brief Returns the body identifier. - static unsigned int GetId() {return 7;} - /// \brief Returns the name of the body. - static std::string GetName() {return "EarthAtm";} - /// \brief Returns the density in g/cm^3 - double density(const GenericTrack&) const; - /// \brief Returns the electron fraction - double ye(const GenericTrack&) const; - /// \brief Returns the radius of the Earth in km. - double GetRadius() const {return radius;} - /// \brief Returns the altitude of the top of the simulated atmosphere in km. - /// This is the altitude from which any initial neutrino flux will begin propagating. - double GetAtmosphereHeight() const {return atm_height;} - /// \brief Set the altitude of the top of the simulated atmosphere. - /// This function invalidates all tracks previously constructed for this body, so be sure to - /// use it before constructing any tracks, or to replace any previously constructed tracks. - /// \param height New height of the top of the atmosphere, in km. - void SetAtmosphereHeight(double height); - - /// \brief Construst a track with the given zenith angle, starting at the top of the atmosphere - /// \param phi Zenith angle with which the track will arrive at the surface of the Earth, - /// possibly after passing through the Earth - Track MakeTrack(double phi); - /// \brief Construst a track with the given zenith angle, starting at the top of the atmosphere - /// \param cosphi Cosine of the zenith angle with which the track will arrive at the surface of - /// the Earth, possibly after passing through the Earth - 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 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 five columns. - /// The first one must run 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 3rd column must correspond to the energy in GeV - /// in descending order over any relevant range - /// The 4th and 5th columns must correspond to the differential E_nu and Mu_nu - /// flux corresponding 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 - 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; - - // registration body - namespace detail{ - class registerBody{ - public: - registerBody(const std::string& name,std::function(hid_t)> fdeserialize); - }; - class registerTrack{ - public: - registerTrack(const std::string& name,std::function(hid_t)> fdeserialize); - }; - } - - std::function(hid_t)> GetBodyDeserializer(std::string); - std::function(hid_t)> GetTrackDeserializer(std::string); - - -// body registration -#define ASW(a, b) a ## b - -#define NUSQUIDS_REGISTER_BODY(classname) \ -namespace{ nusquids::detail::registerBody ASW(body_registerer,classname)(classname::GetName(),classname::Deserialize); nusquids::detail::registerTrack ASW(track_registerer,classname)(classname::Track::GetName(),classname::Track::Deserialize);} - -} // close namespace - -#endif From 22fc70d4a7184eae0fb1a2df60e82743bc130f93 Mon Sep 17 00:00:00 2001 From: ConnorSponsler <33155897+ConnorSponsler@users.noreply.github.com> Date: Sun, 19 Nov 2023 16:56:52 -0500 Subject: [PATCH 11/12] Add files via upload update description of requirements for user supplied data file --- include/nuSQuIDS/body.h | 91 +++++++++++++---------------------------- 1 file changed, 28 insertions(+), 63 deletions(-) diff --git a/include/nuSQuIDS/body.h b/include/nuSQuIDS/body.h index 1848595a..9f575731 100644 --- a/include/nuSQuIDS/body.h +++ b/include/nuSQuIDS/body.h @@ -669,8 +669,6 @@ class EarthAtm: public Body{ double x_ye_min; /// \brief Electron fraction at maximum radius. double x_ye_max; - - marray ESpace; public: /// \brief Default constructor using supplied PREM. @@ -772,85 +770,52 @@ class EarthAtm: public Body{ /// \brief A model of the production of neutrino cosmic ray flux production within the atmosphere. class EmittingEarthAtm: public EarthAtm{ protected: - /// \brief Atm differential nu electron flux array - std::vector diff_enu_prod; - marray enu_flux; - /// \brief Atm differential nu muon flux array - std::vector diff_munu_prod; - marray munu_flux; - /// \brief Atm coszen array - std::vector atm_coszen; - marray coszens; - /// \brief Atm heights array - std::vector atm_heights; - marray heights; - /// \brief Atm energies array - std::vector atm_energy; - marray energies; - /// \brief Minimum coszen. + /// \brief Min/Max coszen. double czen_min; - /// \brief Maximum coszen. double czen_max; - /// \brief number of coszens - long unsigned int czen_number; - /// \brief Minimum height. + /// \brief Min/Max height. double height_min; - /// \brief Maximum height. double height_max; - /// \brief number of heights - long unsigned int height_number; - /// \brief Minimum energy. + /// \brief Min/Max energy. double energy_min; - /// \brief Maximum energy. double energy_max; - /// \brief number of energies - long unsigned int energy_number; + /// \brief vector containing interpolators for the differential flux - std::vector Interpolators; - /// \brief indices for the electron and muon flavors - /// defaulted to the first and second flavor states - int E_nu_index = 0; - int Mu_nu_index = 1; + std::vector> Interpolators; - bool init_emit_energies = false; + /// \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 with user passed energy space - EmittingEarthAtm(marray ESpace); - /// \brief Constructor from a user supplied production model. /// @param prodmodel Path to the production model data file. - /// \details The input file should have five columns. - /// The first one must run in descending order from or between 1 and -1 - /// representing the range of coszen values corresponding to tracks - /// relevant to your propagation + /// \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 3rd column must correspond to the energy in GeV + /// The third column must correspond to the energy in GeV /// in descending order over any relevant range - /// The 4th and 5th columns must correspond to the differential E_nu and Mu_nu - /// flux corresponding to the coszen, height, and energy in the row + /// 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 passes the user's list of energies to the class - /// \details Used by nuSQUIDSATM implementation to provide the - /// class with the energy list required for retrieval of interpolation - void Set_EmissionEnergies(const marray& ESpace_){ - ESpace = ESpace_; - init_emit_energies = true; - } - /// \brief Returns true if neutrino flux emission from bodies is considered - bool Get_EnergyInit(){ - return init_emit_energies; - } /// \brief Function that sets the height of the atmosphere /// \details Used by nuSQUIDSATM to set the height for every /// instance of the body @@ -861,17 +826,17 @@ class EmittingEarthAtm: public EarthAtm{ /// \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 - void Set_ProducedFlavors(int E_nu_index_, int Mu_nu_index_){ - E_nu_index = E_nu_index_; - Mu_nu_index = Mu_nu_index_; + /// 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; - - //Helpful function for data parsing - void unique_sort(std::vector& vec, double& min_val, double& max_val, long unsigned int& number_val); - /// \brief Returns the body identifier. static unsigned int GetId() {return 8;} From fe93c557a74d2b90370e34937d4465f5dc087911 Mon Sep 17 00:00:00 2001 From: ConnorSponsler <33155897+ConnorSponsler@users.noreply.github.com> Date: Sun, 19 Nov 2023 16:59:51 -0500 Subject: [PATCH 12/12] Add files via upload --- include/nuSQuIDS/nuSQuIDS.h | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/include/nuSQuIDS/nuSQuIDS.h b/include/nuSQuIDS/nuSQuIDS.h index 04303548..74d08731 100644 --- a/include/nuSQuIDS/nuSQuIDS.h +++ b/include/nuSQuIDS/nuSQuIDS.h @@ -1281,15 +1281,9 @@ class nuSQUIDSAtm { iinistate = true; } - /// \brief Passes the log scale energy array into the user constructed body - void Set_AtmEmissionEnergies(){ - earth_atm->Set_EmissionEnergies(enu_array); - } - - /// \brief Sets the index of the produced flavors. defaults to 0,1 for nuE and nuMu - void Set_AtmProducedFlavors(int E_nu_index_, int Mu_nu_index_){ - earth_atm->Set_ProducedFlavors(E_nu_index_, Mu_nu_index_); + 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