diff --git a/source/ArraySlice.h b/source/ArraySlice.h new file mode 100644 index 0000000000000000000000000000000000000000..9aa0ac0be82ef7fef44af86e7b24c5bbf241c044 --- /dev/null +++ b/source/ArraySlice.h @@ -0,0 +1,244 @@ +#ifndef ARRAYSLICE_H +#define ARRAYSLICE_H +/** \class ArraySlice + * \brief A class to slice a normal C array + * + * Includes some nice helper functions for its size and mathematical + * operations. + */ + +#include +#include +#include + +template +class ArraySlice { + public: + /** \brief Constructor taking a pointer and a size + * + * The pointer is supposed to be a part of a c array that extends at + * least to pointer+size. + */ + ArraySlice ( T *pointer, size_t size ) : fbegin(pointer), fsize(size), fbeginConst(pointer) {} + ArraySlice ( const T *pointer, size_t size ) :fbegin(NULL), fsize(size), fbeginConst(pointer) {} + + //! Copy constructor, copy the pointer and the size + ArraySlice (const ArraySlice &old) : fsize(old.fsize), fbegin(old.fbegin), fbeginConst(old.fbeginConst) {} + + //! Return the size of the slice + size_t size() const { return fsize; } + + /** \brief Return a reference to element number i in the Slice + * + * That is element pointer+i in the original + */ + T & operator [] (size_t i) { return *(fbegin+i); } + const T & operator [] (size_t i) const { return *(fbeginConst+i); } + + //Assignment operators, work by assigning values to each element + template + ArraySlice & operator = (C number) { + for (size_t i = 0; i < fsize; ++i) { + *(fbegin+i) = number; + } + return (*this); + } + //We need a special assignment operator for same type slices, otherwise + //the useless default is used + ArraySlice & operator = (const ArraySlice & slice) { + if (fsize != slice.fsize) { + std::string errMsg = "Slices need to have same size in assignment"; + std::cerr< + ArraySlice & operator = (const ArraySlice & slice) { + if (fsize != slice.fsize) { + std::string errMsg = "Slices need to have same size in assignment"; + std::cerr< + ArraySlice & operator = (const std::valarray & valarray) { + if (fsize != valarray.size()) { + std::string errMsg = "Valarray and slice need to have same size in assignment"; + std::cerr< + ArraySlice & operator += (C number) { + for (size_t i = 0; i < fsize; ++i){ + *(fbegin+i) += number; + } + return (*this); + } + //Subtraction with a single number + template + ArraySlice & operator -= (C number) { + for (size_t i = 0; i < fsize; ++i){ + *(fbegin+i) -= number; + } + return (*this); + } + //Multiplication with a single number + template + ArraySlice & operator *= (C number) { + for (size_t i = 0; i < fsize; ++i){ + *(fbegin+i) *= number; + } + return (*this); + } + //Division with a single number + template + ArraySlice & operator /= (C number) { + for (size_t i = 0; i < fsize; ++i){ + *(fbegin+i) /= number; + } + return (*this); + } + + //Addition with another slice + template + ArraySlice & operator += (const ArraySlice & slice) { + if (fsize != slice.fsize) { + std::string errMsg = "Slices need to have same size in addition"; + std::cerr< + ArraySlice & operator -= (const ArraySlice & slice) { + if (fsize != slice.fsize) { + std::string errMsg = "Slices need to have same size in subtraction"; + std::cerr< + ArraySlice & operator *= (const ArraySlice & slice) { + if (fsize != slice.fsize) { + std::string errMsg = "Slices need to have same size in multiplication"; + std::cerr< + ArraySlice & operator /= (const ArraySlice & slice) { + if (fsize != slice.fsize) { + std::string errMsg = "Slices need to have same size in division"; + std::cerr< + ArraySlice & operator += (const std::valarray & valarray) { + if (fsize != valarray.size()) { + std::string errMsg = "Valarray and slice need to have same size in addition"; + std::cerr< + ArraySlice & operator -= (const std::valarray & valarray) { + if (fsize != valarray.size()) { + std::string errMsg = "Valarray and slice need to have same size in subtraction"; + std::cerr< + ArraySlice & operator *= (const std::valarray & valarray) { + if (fsize != valarray.size()) { + std::string errMsg = "Valarry and slice need to have same size in multiplication"; + std::cerr< + ArraySlice & operator /= (const std::valarray & valarray) { + if (fsize != valarray.size()) { + std::string errMsg = "Valarray and slice need to have same size in division"; + std::cerr< +std::ostream & operator << (std::ostream &os, const ArraySlice & arr){ + os<<"["; + for (size_t i = 0; i < arr.size()-1; ++i){ + os<//AWS20050624 +#include //AWS20050624 +#include //IMOS20020112 + +class Configure +{ + public: + + char *galdef_directory; + char *fits_directory; + char *adjunct_directory; + + int directory_length; + +//interface function prototype + int init(); + +//AWS20100409: from galprop for read_isrf.cc + std::string fGaldefDirectory; + std::string fFITSDataDirectory; + std::string fOutputDirectory; + std::string fOutputPrefix; + std::string fVersion; + std::string fGlobalDataPath; + +}; + + + +#endif + + diff --git a/source/Convolved.h b/source/Convolved.h new file mode 100644 index 0000000000000000000000000000000000000000..95e6a491de80c462dc9217fb36c1dc9efad526ac --- /dev/null +++ b/source/Convolved.h @@ -0,0 +1,76 @@ +//energy integrated, instrument convolved model skymaps + +#include"Distribution.h" +#include"HIH2IC.h" //AWS20080526 +#include"Sources.h" //AWS20080812 +#include"Skymap.h" //AWS20080604 + +class Convolved +{ + public: + + Distribution EGRET_bremss_skymap; + Distribution *EGRET_IC_iso_skymap; // 3 ISRF components + Distribution *EGRET_IC_aniso_skymap; // 3 ISRF components //AWS20060907 + Distribution EGRET_pi0_decay_skymap; + Distribution EGRET_isotropic_skymap; + Distribution EGRET_total_skymap; + + Distribution EGRET_HIR_skymap; //AWS20041223 + Distribution EGRET_H2R_skymap; //AWS20041223 + + Distribution EGRET_HIT_skymap; //AWS20050103 + Distribution EGRET_H2T_skymap; //AWS20050103 + + + Distribution GLAST_bremss_skymap; //AWS20080514 + Distribution *GLAST_IC_iso_skymap; // 3 ISRF components + Distribution *GLAST_IC_aniso_skymap; // 3 ISRF components + Distribution GLAST_pi0_decay_skymap; + Distribution GLAST_isotropic_skymap; + Distribution GLAST_total_skymap; + + Distribution GLAST_HIR_skymap; + Distribution GLAST_H2R_skymap; + + Distribution GLAST_HIT_skymap; + Distribution GLAST_H2T_skymap; + + +// GLAST models from gardian +// has to be a pointer apparently (otherwise galplot.cc gives compile error: "no matching function for Galplot()" ) +// something to do with how inheritance from BaseModel works ? +// JoinModel.cc also uses (a vector of) pointers to models +// HIH2IC GLAST_model_HIH2IC; //AWS20080526 NO + HIH2IC *GLAST_model_HIH2IC; //AWS20080526 YES + + + + + Skymap GLAST_convolved_counts; //AWS20080604 + Skymap GLAST_convolved_counts_over_exposure;//AWS20080604 + + Skymap GLAST_convolved_counts_HI; //AWS20080606 + Skymap GLAST_convolved_counts_H2; //AWS20080606 + Skymap GLAST_convolved_counts_IC; //AWS20080606 + Skymap GLAST_convolved_counts_HII; //AWS20090810 + + + Skymap GLAST_convolved_counts_sources; //AWS20081209 + Skymap GLAST_convolved_intensity_sources; //AWS20081212 + + Skymap GLAST_convolved_counts_sourcepop_sublimit; //AWS20090109 + Skymap GLAST_convolved_counts_sourcepop_soplimit; //AWS20090109 + + + Skymap GLAST_convolved_intensity_sourcepop_sublimit; //AWS20090109 + Skymap GLAST_convolved_intensity_sourcepop_soplimit; //AWS20090109 + + + Skymap GLAST_convolved_counts_solar_IC; //AWS20100215 + Skymap GLAST_convolved_intensity_solar_IC; //AWS20100215 + + Skymap GLAST_convolved_counts_solar_disk; //AWS20100215 + Skymap GLAST_convolved_intensity_solar_disk; //AWS20100215 + +}; diff --git a/source/Coordinate.h b/source/Coordinate.h new file mode 100644 index 0000000000000000000000000000000000000000..82352aed3007341f4811a11937fed3abe55d0652 --- /dev/null +++ b/source/Coordinate.h @@ -0,0 +1,47 @@ +#ifndef Coordinate_h +#define Coordinate_h +#include "PhysicalConstants.h" +#include "pointing.h" + +namespace SM { + /** \brief Coordinates for easy access of skymaps + * + * Handles conversion between conventional l and b with the equator at b = 0 to healpix coordinates + */ + class Coordinate { + private: + double m_l, m_b; //!< The coordinates + public: + /**\brief Constructor that takes in galactic coordinates (l,b) in + * degrees. + * + * \param l the longitude in degrees, GC at 0 + * \param b the latitude in degrees, GC at 0 + */ + Coordinate(const double l, const double b) : m_l(l), m_b(b){} + /**\brief Default constructor initializes to 0 */ + Coordinate() : m_l(0), m_b(0) {} + /**\brief Constructor that takes in healpix coordinate pointing */ + Coordinate(const pointing & point); + + /** \brief Return the value of l */ + const double & l() const {return m_l;} + /** \brief Return the value of b */ + const double & b() const {return m_b;} + + /**\brief Get the angular coordinate for healpix + * + * \return a pointing object to use with healpix + */ + pointing healpixAng(void) const; + + /**\brief Output operator + * + * The format is (l,b) + */ + friend std::ostream & operator << (std::ostream &os, const Coordinate & coord); + + }; +} + +#endif diff --git a/source/Counts.h b/source/Counts.h new file mode 100644 index 0000000000000000000000000000000000000000..666e8e9cafb4276270b26104df9b0e08ed6cd0da --- /dev/null +++ b/source/Counts.h @@ -0,0 +1,169 @@ +#ifndef COUNTSMAP_H +#define COUNTSMAP_H + +#include +#include +#include "Skymap.h" +#include +#include "Parameters.h" + +/** \brief Counts map serves as a storage for binned counts. + * + * It has the ability to read counts from both pre-binned healpix based maps, + * or automatically bin counts from FT1 files. + */ +class CountsMap { + public: + /** \brief Standard constructor, basically does nothing */ + CountsMap(){} + /** \brief Construct a counts map from a Parameters object. + * + * \param pars should contain the following parameters + * - countsFile: name of the input file + * - countOrder: the order of the resulting map + * - energyBinning: type of energy binning (log, linear or list) + * - nEnergyBins: number of energy bins + * - eMin: minimum energy in MeV (a space separated list of values if + * energyBinning is list) + * - eMax: maximum energy in MeV (a space separated list of values if + * energyBinning is list) + * - FT1filter: a filter used when opening the FT1 fits file (optional) + * - countsOutFile: name of an output file for the binned FT1 files + * (optional) + * + * If countOrder is not given, the file is assumed to be of healpix type + * and all other parameters are unrelevant. If it is not given, file + * should either be a FT1 fits file or a text file containing a list of FT1 + * file names, which will all be included in the binning. + */ + CountsMap(const Parameters &pars); + /** \brief Construct an empty counts map of given order with given energy + * boundaries. + * + * \param order gives the order of the healpix map + * \param eMin and \param eMax are valarrays containing the lower and upper + * boundaries, respecitvely, of the energy bins in MeV. + */ + CountsMap(int order, std::valarray & eMin, std::valarray & eMax); + /** \brief Construct a counts map from a healpix file. + * + * \param fileName should give the location of a healpix fits file. The + * data should be in an extension name counts, containing a table with a + * single vector column, where each row is a single healpix pixel and the + * vector contains the count as a function of energy. There should be + * another extension, EBINS, holding a table with information about the + * energy bins. That table should have three columns, Channel for the + * channel number, E_Min for the minimum energy and E_Max for the maximum + * energy. + * + * This is the output format when counts map is used to bin FT1 files. + */ + CountsMap(const std::string &fileName); + /** \brief Construct a counts map from a FT1 fits file. + * + * \param fileName is the name of the fits file, or a name of a text file * containing a list of FT1 filenames, one per line. + * \param order controls * the order of the healpix map and + * \param eMin lower energy boundaries + * \param eMax upper energy boundaries + * \param filter an optional parameter that filters the table on opening. This string should follow the cfitsio + * convention for extended file syntax, without the brackets. Example: \code filter = "MC_SRC_ID == 13000" \endcode + */ + CountsMap(const std::string &fileName, int order, std::valarray & eMin, std::valarray & eMax, const std::string &filter = ""); + /** \brief Copy constructor */ + CountsMap(const CountsMap & oldCountsMap); + /** \brief Add photons to the countmap from an valarray of coordinates, + * energy pairs. + * + * \param photonArray is a valarray of coordinate, energy pairs, where the + * energy is given in MeV. + */ + void addPhotons(const std::valarray< std::pair< SM::Coordinate, double > > & photonArray); + /** \brief Add the contents of an FT1 file, or a text file with a list of + * FT1 files to the counts map. + * + * \param fileName should point to the FT1 file or the text file with a + * list of FT1 files, one per line. An optional parameter \param filter + * will be applied when opening the FT1 file(s). It should follow the extended + * file name syntax of cfitsio. + */ + void addFile(const std::string &fileName, const std::string & filter = ""); + /** \brief Reset the counts map to different order and energy bins. + * + * Destroys all data. \param order gives the healpix order while \param + * eMin and \param eMax give energy boundaries. + */ + void reset(int order, std::valarray &eMin, std::valarray & eMax); + /** \brief Return a constant reference to the skymap. + * + * The energy stored in the skymap object refers to the center value of the + * energy bin. + */ + // const Skymap & getCountsMap() const; + const Skymap & getCountsMap() const; //AWS20091002 + /** \brief Return a const reference to the lower energy boundaries */ + const std::valarray & getEMin() const; + /** \brief Return a const reference to the upper energy boundaries */ + const std::valarray & getEMax() const; + /** \brief Return the number of photons lost in binning. + * + * When binning the data, the energy of the photon can be out of range for + * the energy bins. This returns the numbers of photon that where in the + * FT1 files, but did not end in the count map. This does not take into + * account any filter applied, since that is done prior to binning. + */ + int lost() const; + /** \brief Write the counts map to a fits file. + * + * \param fileName is the name of the output file. It will be overwritten + * without an error message. The format is the same as described in the + * healpix fits file constructor. + */ + void write(const std::string &fileName) const; + /** \brief A general exception class for the object to throw */ + class CountsMapException : public std::exception { + private: + std::string eMessage; + public: + /** \brief Construct an error with the default message "Error in CountsMap class" */ + CountsMapException() : eMessage("Error in CountsMap class"){} + /** \brief Construct an error with a specialized error message + * + * \param message is the specialized message + */ + CountsMapException(const std::string & message) : eMessage(message){} + ~CountsMapException() throw(){} + /** \brief Return the error message */ + virtual const char* what () const throw(){ + return(eMessage.c_str()); + } + }; + + /** \brief The assignment operator */ + CountsMap & operator = (const CountsMap & oldCountsMap); + + private: + /** \brief Read an FT1 file. + * + * \param fileName should give the name of file to read and it should be an + * FT1 fits file. If an error occurs in reading the file, the routine + * returns a value different from 0. \param filter is an optional + * parameter for filtering the FT1 file. It should follow the cfitsio + * extended file name syntax. + */ + void ReadFT1Fits(const std::string &fileName, const std::string & filter = ""); + /** \brief Read an healpix fits file. + * + * \param fileName should give the name of the file to read and it should + * be in a compatible format, \see CountsMap::CountsMap. Uses CCFITS which + * throws errors if something goes wrong. + */ + void ReadHealpixFits(const std::string &fileName); + /** \brief Add a coordinate,energy pair to the counts map */ + void Add(const std::pair & photon); + // Skymap fDataSkymap; //!< Stores the counts data + Skymap fDataSkymap; //!< Stores the counts data AWS20091002 + std::valarray fEMin, fEMax; //!< The skymap only stores single energy value, so we need to implement the bins + int fLost; //!< Number of lost photons in the binning, not counting filters + std::string fEnergyCol; //!< The name of the energy column +}; +#endif diff --git a/source/Data.h b/source/Data.h new file mode 100644 index 0000000000000000000000000000000000000000..053fa91408ea0e15d140f04792d87668077c999d --- /dev/null +++ b/source/Data.h @@ -0,0 +1,211 @@ + + + + +#include"Distribution.h" + +// GLAST classes from gardian +#include"Counts.h" //AWS20080519 +#include"Exposure.h" //AWS20080519 +#include"Psf.h" //AWS20080519 +#include"Skymap.h" //AWS20080523 + +#include "healpix_map.h" //AWS20110413 + +class Data +{ + public: + + int n_E_EGRET; + double *E_EGRET; + + double EGRET_psf_dtheta; + + Distribution EGRET_exposure; + Distribution EGRET_counts; + Distribution EGRET_psf; + + + int n_E_GLAST; //AWS20080508 + double *E_GLAST; //AWS20080508 + Distribution GLAST_exposure; //AWS20080508 + Distribution GLAST_counts; //AWS20080508 + Distribution GLAST_psf; //AWS20080508 + double GLAST_psf_dtheta; //AWS20080508 + + // GLAST data in healpix + + CountsMap GLAST_counts_healpix; //AWS20080519 + Exposure GLAST_exposure_healpix; //AWS20080519 + Psf GLAST_psf_healpix; //AWS20080519 + SkymapGLAST_exposure_integrated; //AWS20080523 + int GLAST_exposure_integrated_init;//AWS20080523 + + int n_E_COMPTEL; //AWS20030910 + double *E_COMPTEL; //AWS20030910 + Distribution COMPTEL_intensity; //AWS20030910 + + // original survey data + + Distribution synchrotron____10MHz;//AWS20050728 + Distribution synchrotron____22MHz;//AWS20050728 + Distribution synchrotron____45MHz;//AWS20050728 + Distribution synchrotron___150MHz;//AWS20070315 + Distribution synchrotron___408MHz;//AWS20050728 + Distribution synchrotron___820MHz;//AWS20070209 + Distribution synchrotron__1420MHz;//AWS20050728 + Distribution synchrotron__2326MHz;//AWS20070219 + Distribution synchrotron_22800MHz;//AWS20050728 + Distribution synchrotron_33000MHz;//AWS20070226 + Distribution synchrotron_41000MHz;//AWS20070226 + Distribution synchrotron_61000MHz;//AWS20070226 + Distribution synchrotron_94000MHz;//AWS20070226 + + // rebinned on galprop grid + Distribution synchrotron_skymap____10MHz;//AWS20070112 + Distribution synchrotron_skymap____22MHz;//AWS20070112 + Distribution synchrotron_skymap____45MHz;//AWS20070112 + Distribution synchrotron_skymap___150MHz;//AWS20070315 + Distribution synchrotron_skymap___408MHz;//AWS20070112 + Distribution synchrotron_skymap___820MHz;//AWS20070209 + Distribution synchrotron_skymap__1420MHz;//AWS20070112 + Distribution synchrotron_skymap__2326MHz;//AWS20070219 + Distribution synchrotron_skymap_22800MHz;//AWS20070112 + Distribution synchrotron_skymap_33000MHz;//AWS20070226 + Distribution synchrotron_skymap_41000MHz;//AWS20070226 + Distribution synchrotron_skymap_61000MHz;//AWS20070226 + Distribution synchrotron_skymap_94000MHz;//AWS20070226 + + + // non-WMAP surveys as Healpix_Map's + Healpix_Map synchrotron_Healpix__408MHz;//AWS20130121 from lambda website + + // WMAP 7-year polarization data + Healpix_Map synchrotron_WMAP_I_22800MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_I_33000MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_I_41000MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_I_61000MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_I_94000MHz;//AWS20110413 + + + Healpix_Map synchrotron_WMAP_Q_22800MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_Q_33000MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_Q_41000MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_Q_61000MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_Q_94000MHz;//AWS20110413 + + Healpix_Map synchrotron_WMAP_U_22800MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_U_33000MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_U_41000MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_U_61000MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_U_94000MHz;//AWS20110413 + + Healpix_Map synchrotron_WMAP_P_22800MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_P_33000MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_P_41000MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_P_61000MHz;//AWS20110413 + Healpix_Map synchrotron_WMAP_P_94000MHz;//AWS20110413 + + // rebinned on galprop grid + Distribution synchrotron_I_skymap_22800MHz;//AWS20110426 + Distribution synchrotron_I_skymap_33000MHz;//AWS20110426 + Distribution synchrotron_I_skymap_41000MHz;//AWS20110426 + Distribution synchrotron_I_skymap_61000MHz;//AWS20110426 + Distribution synchrotron_I_skymap_94000MHz;//AWS20110426 + + + Distribution synchrotron_Q_skymap_22800MHz;//AWS20110426 + Distribution synchrotron_Q_skymap_33000MHz;//AWS20110426 + Distribution synchrotron_Q_skymap_41000MHz;//AWS20110426 + Distribution synchrotron_Q_skymap_61000MHz;//AWS20110426 + Distribution synchrotron_Q_skymap_94000MHz;//AWS20110426 + + Distribution synchrotron_U_skymap_22800MHz;//AWS20110426 + Distribution synchrotron_U_skymap_33000MHz;//AWS20110426 + Distribution synchrotron_U_skymap_41000MHz;//AWS20110426 + Distribution synchrotron_U_skymap_61000MHz;//AWS20110426 + Distribution synchrotron_U_skymap_94000MHz;//AWS20110426 + + Distribution synchrotron_P_skymap_22800MHz;//AWS20110426 + Distribution synchrotron_P_skymap_33000MHz;//AWS20110426 + Distribution synchrotron_P_skymap_41000MHz;//AWS20110426 + Distribution synchrotron_P_skymap_61000MHz;//AWS20110426 + Distribution synchrotron_P_skymap_94000MHz;//AWS20110426 + + // WMAP 7-year MEM templates + + // WMAP 7-year free-free MEM templates + Healpix_Map free_free_WMAP_MEM_22800MHz;//AWS20110621 + Healpix_Map free_free_WMAP_MEM_33000MHz;//AWS20110621 + Healpix_Map free_free_WMAP_MEM_41000MHz;//AWS20110621 + Healpix_Map free_free_WMAP_MEM_61000MHz;//AWS20110621 + Healpix_Map free_free_WMAP_MEM_94000MHz;//AWS20110621 + + // skymaps for all frequencies derived from WMAP 7-year free-free MEM templates + // keep consistent nomenclature: append _hp_skymap if Skymap, _skymap if Distribution + Skymap free_free_WMAP_MEM_hp_skymap; //AWS20110621 + Distribution free_free_WMAP_MEM_skymap; //AWS20110621 + + + + // WMAP 7-year dust MEM templates + Healpix_Map dust_WMAP_MEM_22800MHz;//AWS20111221 + Healpix_Map dust_WMAP_MEM_33000MHz;//AWS20111222 + Healpix_Map dust_WMAP_MEM_41000MHz;//AWS20111222 + Healpix_Map dust_WMAP_MEM_61000MHz;//AWS20111222 + Healpix_Map dust_WMAP_MEM_94000MHz;//AWS20111222 + + // skymaps for all frequencies derived from WMAP 7-year dust MEM templates + + Skymap dust_WMAP_MEM_hp_skymap; //AWS20111222 + Distribution dust_WMAP_MEM_skymap; //AWS20111222 + + + + // WMAP 7-year synchrotron MEM templates + Healpix_Map synchrotron_WMAP_MEM_22800MHz;//AWS20111221 + Healpix_Map synchrotron_WMAP_MEM_33000MHz;//AWS20111222 + Healpix_Map synchrotron_WMAP_MEM_41000MHz;//AWS20111222 + Healpix_Map synchrotron_WMAP_MEM_61000MHz;//AWS20111222 + Healpix_Map synchrotron_WMAP_MEM_94000MHz;//AWS20111222 + + // skymaps for all frequencies derived from WMAP 7-year synchrotron MEM templates + + Skymap synchrotron_WMAP_MEM_hp_skymap; //AWS20111222 + Distribution synchrotron_WMAP_MEM_skymap; //AWS20111222 + + + + + + // WMAP 7-year MCMC spinning dust model templates + // Gold etal. 2011 ApJSupp 192,15. + // these are provided at one frequency, with power-law fixed at +2.0 for dust + // and beta_d=+2.0, nu_sd=4.2 GHz in equation(10) for spinning dust. + + Healpix_Map dust_WMAP_MCMC_SD_94000MHz; //AWS20120116 + Healpix_Map dust_Q_WMAP_MCMC_SD_94000MHz; //AWS20120117 + Healpix_Map dust_U_WMAP_MCMC_SD_94000MHz; //AWS20120117 + + Healpix_Map spin_dust_WMAP_MCMC_SD_22800MHz; //AWS20120116 + + // skymaps for all frequencies derived from WMAP 7-year MCMC spinning dust model templates + + Skymap dust_WMAP_MCMC_SD_hp_skymap; //AWS20120116 + Distribution dust_WMAP_MCMC_SD_skymap; //AWS20120116 + + Skymap dust_Q_WMAP_MCMC_SD_hp_skymap; //AWS20120117 + Distribution dust_Q_WMAP_MCMC_SD_skymap; //AWS20120117 + + Skymap dust_U_WMAP_MCMC_SD_hp_skymap; //AWS20120117 + Distribution dust_U_WMAP_MCMC_SD_skymap; //AWS20120117 + + Skymap dust_P_WMAP_MCMC_SD_hp_skymap; //AWS20120117 + Distribution dust_P_WMAP_MCMC_SD_skymap; //AWS20120117 + + + Skymap spin_dust_WMAP_MCMC_SD_hp_skymap; //AWS20120116 + Distribution spin_dust_WMAP_MCMC_SD_skymap; //AWS20120116 + + + }; diff --git a/source/Distribution.h b/source/Distribution.h new file mode 100644 index 0000000000000000000000000000000000000000..129cfb857f37cc4fb6b0ec54eafbe6fe8678cae3 --- /dev/null +++ b/source/Distribution.h @@ -0,0 +1,66 @@ + +//**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****| +// * Distribution.h * galprop package * 4/14/2000 +//**"****!****"****!****"****!****"****!****"****!****"****!****"****!****"****| + +#ifndef Distribution_h +#define Distribution_h +using namespace std; //AWS20050919 +#include //AWS20050919 +#include //AWS20050919 +#include //IMOS20020112 //AWS20050919 +#include //IMOS20020112 //AWS20050919 +#include"Spectrum.h" + +class Distribution +{ + public: + + int n_pgrid; // number of points in momentum + int n_rgrid; // number of points in radius (2D) + int n_zgrid; // number of points in z (1D,2D) + int n_xgrid; // number of points in x (3D) + int n_ygrid; // number of points in y (3D) + + int n_spatial_dimensions;// 1,2,3D + + Spectrum *d1; // 1D array + Spectrum **d2; // 2D array + Spectrum ***d3; // 3D array + + Distribution(){}; // needed for initial construction +// ~Distribution(); // needed for destruction + Distribution(int n_rgrid_,int n_zgrid_,int n_pgrid_); + Distribution(int n_xgrid_,int n_ygrid_,int n_zgrid_,int n_pgrid_); + void delete_array(); + void init(int n_rgrid_,int n_zgrid_,int n_pgrid_); + void init(int n_xgrid_,int n_ygrid_,int n_zgrid_,int n_pgrid_); + void print(); + void print(double units); + float max(); + Distribution operator=(double); + Distribution operator+(double); + Distribution operator+=(double); + Distribution operator*(double); + Distribution operator*=(double); + Distribution operator/=(double); + Distribution operator =(Distribution); + Distribution operator+=(Distribution); + Distribution operator*=(Distribution); + Distribution operator/=(Distribution); +}; + +// following cases have to be non-member functions since they have 2 arguments; +// therefore their prototypes are defined outside the class + +Distribution operator+(double,Distribution); +Distribution operator*(double,Distribution); +Distribution operator+(Distribution,Distribution); +Distribution operator*(Distribution,Distribution); +Distribution operator/(Distribution,Distribution); + +#endif + + + + diff --git a/source/EGRET_catalogue.h b/source/EGRET_catalogue.h new file mode 100644 index 0000000000000000000000000000000000000000..bbfa07a7a9317979757333912581c0a42ffe8eaf --- /dev/null +++ b/source/EGRET_catalogue.h @@ -0,0 +1,18 @@ + + + + +class EGRET_catalogue +{ +public: + int n_source; + double *ra,*dec,*l,*b,*theta,*F,*dF,*index,*counts,*sqrtTS; + char **name,**vp,**ID; + + int *selected; + + int read(int debug); + int print(ofstream &txt_stream); + + int select(int options); +}; diff --git a/source/EnergyDispersionFermiLAT.h b/source/EnergyDispersionFermiLAT.h new file mode 100644 index 0000000000000000000000000000000000000000..b7641352034bc4cf030a7033dc9624c34ee80a7e --- /dev/null +++ b/source/EnergyDispersionFermiLAT.h @@ -0,0 +1,82 @@ +#include"Fermi_Aeff.h" + +class EnergyDispersionFermiLAT +{ + public: + + int nE_true,nE_meas; + + + // for thibaut method +//valarray F,S1,K1,BIAS,S2,K2,PINDEX1,PINDEX2; //AWS20150220 + valarray BIAS; //AWS20150220 for compatibility, to be removed later + + valarray F,S1,K1,BIAS1,S2,K2,BIAS2,PINDEX1,PINDEX2; //AWS20150220 + + valarray parameters_E_low, parameters_E_high; + int n_parameter_sets; + double parameters_E_min,parameters_E_factor; + + valarray Sd_parameters; //AWS20150220 + + double E_true_threshold; //AWS20150223 + double E_meas_threshold; //AWS20150301 + + int parameters_initialized; + int EnergyDispersionMatrix2_initialized; + string method; + int use_Aeff; + + string parameter_file_name; + string parameter_file_type; //AWS20150220 + string exposure_file_name; + string exposure_file_type; + + vector > EnergyDispersionMatrix2; + Fermi_Aeff fermi_Aeff; + int use_matrix; + + // for rsp method + vector > EnergyDispersionMatrix; + + + valarray E_true, E_meas; + double log_E_true_min, log_E_meas_min; + valarray dE_true , dE_meas ; + double dlog_E_true, dlog_E_meas; + void read(int debug); + + int read_parameters(int debug); + int read_exposure (int debug); + int set_method(string method_); + int set_use_Aeff(int use_Aeff_); + int set_E_true_threshold(double E_true_threshold_); //AWS20150223 + int set_E_meas_threshold(double E_meas_threshold_); //AWS20150301 + + int set_parameter_file (string parameter_file_name_); + int set_parameter_file_type(string parameter_file_type_); //AWS20150220 + int set_exposure_file (string exposure_file_name_); + int set_exposure_file_type (string exposure_file_type_); + + int set_use_matrix(int use_matrix_); + int reset(); + + double value(double E_true,double E_meas, int debug); + + double value_thibaut(double E_true,double E_meas, int debug); + double g(double x, double sigma, double k, double b, double p,int debug); + double Sd(double E,int debug); + + void ApplyEnergyDispersion(valarray E, valarray &spectrum, int debug); + void ApplyEnergyDispersion(valarray E, valarray &spectrum, double E_interp_factor,int debug); + void ApplyEnergyDispersion(valarray E_true_, valarray spectrum_true,valarray E_meas_, valarray &spectrum_meas, int debug); + + void test(); + int initialized; // for rsp method + + + + + + +}; diff --git a/source/ErrorLogger.h b/source/ErrorLogger.h new file mode 100644 index 0000000000000000000000000000000000000000..42811d3e2e491a46e54bdf0e13beede70c70b98a --- /dev/null +++ b/source/ErrorLogger.h @@ -0,0 +1,143 @@ +/** \file + Interface to the Auger Error logger + \author Lukas Nellen + \version $Id: ErrorLogger.h 6242 2007-08-16 13:28:12Z darko $ + \date 24-Jan-2003 +*/ + +#ifndef _utl_ErrorLogger_h_ +#define _utl_ErrorLogger_h_ + +#include +#include +#include + +#include + +namespace utl { + + class ErrorLogger : public LeakingSingleton { + public: + /// Message severity levels + enum SeverityLevel { + eDebug = -1, ///< Debugging message + eInfo = 0, ///< General (informational) message + eWarning = 1, ///< Warning message + eError = 2, ///< Error message + eFatal = 3 ///< Fatal error message + }; + + /// Message verbosity + enum VerbosityLevel { + eDefaultVerbosity = -1, ///< Use the default verbosity + eTerse, ///< Terse error messages + eVerbose, ///< Verbose error messages + eSilent ///< Silent messages (turned off) + }; + + /// General interface for logging a message + void Log(const SeverityLevel severity, + const std::string functionName, + const std::string fileName, + const int lineNumber, + const std::string message, + VerbosityLevel verbosity = eDefaultVerbosity) const; + + /// \brief Dump message from \c ostringstream instead of of a \c string + void Log(const SeverityLevel severity, + const std::string functionName, + const std::string fileName, + const int lineNumber, + const std::ostringstream& message, + const VerbosityLevel verbosity = eDefaultVerbosity) const + { Log(severity, functionName, fileName, lineNumber, message.str(), verbosity); } + + /** + \brief Set the error logging stream + + \note + This method will be replaced by a more general mechanism for the + routing of error messages. + */ + void SetStream(std::ostream& stream) { fOStream = &stream; } + + /// Set the verbosity level + void SetVerbosity(const VerbosityLevel verbosity) + { if (verbosity != eDefaultVerbosity) fVerbosity = verbosity; } + + private: + ErrorLogger() : fOStream(0), fVerbosity(eTerse) { } + ErrorLogger(const ErrorLogger& er); + ErrorLogger& operator=(const ErrorLogger& rhs); + + std::ostream* fOStream; ///< Current stream for logging messages + VerbosityLevel fVerbosity; ////< Verbosity level + + friend class LeakingSingleton; + + }; + +} // namespace utl + + +/** + \brief Standard message logging macro + + This macro is used by the convenience macros defined below to write + a log message. It automatically sets the function name, the file + name and the line number. +*/ + +#define LOG_MESSAGE_(severity, message) \ +utl::ErrorLogger::GetInstance().Log(severity, __func__, __FILE__, __LINE__, \ + message) + +/// Brief message logging macro - always do \e terse logging +#define LOG_TERSE_MESSAGE_(severity, message) \ +utl::ErrorLogger::GetInstance().Log(severity, __func__, __FILE__, __LINE__, \ + message, utl::ErrorLogger::eTerse) + +/** + \brief Macro for logging debugging messages. + + This macro is only active if the \c DEBUG macro is + defined. Otherwise, debug messages are suppressed at compile time. + + \remark This macro is not named \c DEBUG since \c DEBUG is usually + used to activate conditionally compiled debugging code. +*/ +#ifdef DEBUG +# define DEBUGLOG(message) LOG_MESSAGE_(utl::ErrorLogger::eDebug, message) +#else +# define DEBUGLOG(message) +#endif + +/// Macro for logging informational messages +#define INFO(message) LOG_MESSAGE_(utl::ErrorLogger::eInfo, message) +/// Macro for logging warning messages +#define WARNING(message) LOG_MESSAGE_(utl::ErrorLogger::eWarning, message) +/// Macro for logging error messages +#define ERROR(message) LOG_MESSAGE_(utl::ErrorLogger::eError, message) +/// Macro for logging fatal messages +#define FATAL(message) LOG_MESSAGE_(utl::ErrorLogger::eFatal, message) + +/// Macro for logging informational messages +#define INFO_TERSE(message) LOG_TERSE_MESSAGE_(utl::ErrorLogger::eInfo, \ + message) +/// Macro for logging warning messages +#define WARNING_TERSE(message) LOG_TERSE_MESSAGE_(utl::ErrorLogger::eWarning, \ + message) +/// Macro for logging error messages +#define ERROR_TERSE(message) LOG_TERSE_MESSAGE_(utl::ErrorLogger::eError, \ + message) +/// Macro for logging fatal messages +#define FATAL_TERSE(message) LOG_TERSE_MESSAGE_(utl::ErrorLogger::eFatal, \ + message) + +#endif // _utl_ErrorLogger_h_ + +// Configure (x)emacs for this file ... +// Local Variables: +// mode: c++ +// compile-command: "make -k -C .. ErrorLogger/ErrorLogger.o" +// End: diff --git a/source/Exposure.h b/source/Exposure.h new file mode 100644 index 0000000000000000000000000000000000000000..7161adfa46d1364e3ac3596e842d4fc78a81c3b2 --- /dev/null +++ b/source/Exposure.h @@ -0,0 +1,151 @@ +#ifndef EXPOSURE_H +#define EXPOSURE_H + +#include "Skymap.h" +#include "ArraySlice.h" + +#include +#include +#include + +/** \brief Exposure stores the exposure maps in appropriate format. + * + * Initializes the healpix exposure map from either a healpix fits file or a mapcube fits file. + * Also does Lagrange polynomial approximation for the exposure between energy points. Does not + * handle exposure correction, that is done by the BaseModel class. + */ +class Exposure { + public: + /** \brief Standard constructor, does nothing */ + Exposure(){} + /** \brief Construct an exposure object from fits mapcube. + * + * \param fileName is the name of the fits file to be opened. This should + * be either a Skymap formatted file or a mapcube in galactic coordinates using the CAR projection. A full sky + * map is needed, but the program will continue with a warning if it finds + * something wrong. The energy dependance of the mapcube should be given + * in a separate extension called ENERGIES, containing a table with a + * single column of energy in MeV. + * + * If the exposure is given as a mapcube, it will be rebinned to healpix + */ + Exposure(const std::string &fileName); + /** \brief Read a file to create the exposure. + * + * \param fileName is the name of file to be read. Its format can be + * either a mapcube as described in the constructor or a Skymap fits file. + */ + void readFile(const std::string &fileName); + /** \brief Write the exposure out in Skymap format + * + * \param fileName is the name of the output file, which will be + * overwritten. + */ + void writeFile(const std::string &fileName); + /** \brief Create a spatially flat exposure with energy dependence + * + * \param order is the healpix order + * \param energies is the energy binning required for the exposure + * \param functionPars is a valarray of function parameters, different size + * is expected for different functions + * \param function is a string representing the energy dependance needed. + * Three values are accepted and the names in the parenthesis indicate the + * parameters needed in functionPars array: + * - Linear (slope, constant) + * - Powerlaw (index, prefactor) + * - Array (exposure at every energy) + */ + void createTestMap(int order, const std::valarray & energies, const std::valarray & functionPars, const std::string & function); + /** \brief Return a constant reference to the exposure skymap. */ + const Skymap & getExposureMap() const; + /** \brief Return a power law weighted exposure for a given coordinate and energy bin + * + * \param co is the coordinate for which the exposure is evaluated. + * \param eMin is the minimum energy boundary + * \param eMax is the maximum energy boundary + * \param index is the power law index (note that it is not negated) + * + * \return the exposure weighted with the power law index over the given energy bin. + */ + double getWeightedExposure(SM::Coordinate &co, double eMin, double eMax, double index) const; + /** \brief Convert an intensity map (units of cm^-2 s^-1 sr^-1 + * MeV^-1, given that exposure energy is MeV) to counts. + * + * This method uses a semi-analytical integration, using a + * power law interpolation on the input map. + * + * \param inMap is the intensity map in units of cm^-2 s^-1 * sr^-1 MeV^-1, given that exposure energy is MeV + * \param eMin gives the minimum energy boundary, has to be the + * same size as eMax + * \param eMax gives the maximum energy boundary + * \param order is the output order. Defaults to the order of + * inMap + * + * \return a double precision skymap of the expected counts of + * inMap + */ + Skymap mapToCounts(Skymap inMap, const std::valarray &eMin, const std::valarray &eMax, int order = -1) const; + /** \brief Convert a spectra (units of cm^-2 s^-1 + * MeV^-1, given that exposure energy is MeV) at a given + * coordinate to counts. + * + * This method uses a semi-analytical integration, using a + * power law interpolation on the input spectra. + * + * \param intensities is the intensity spectra in units of cm^-2 s^-1 MeV^-1, given that exposure energy is MeV + * \param energies is the energies corresponding to the + * intensities. Has to be the same size as intensities. + * \param eMin gives the minimum energy boundary, has to be the + * same size as eMax + * \param eMax gives the maximum energy boundary + * + * \return a double precision valarray of the expected counts of + * the spectra + */ + std::valarray spectraToCounts(const std::valarray &intensities, const std::valarray &energied, const std::valarray &eMin, const std::valarray &eMax, const SM::Coordinate &co) const; + /** \brief Convert a power law spectra (units of cm^-2 s^-1 + * MeV^-1, given that exposure energy is MeV) at a given + * coordinate to counts. + * + * \param index is the power law index + * \param prefactor is the intensity of the spectra at the pivot energy. It is in units of cm^-2 s^-1 sr^-1 MeV^-1, given that exposure energy is MeV + * \param energies is the energies corresponding to the + * intensities. Has to be the same size as intensities. + * \param eMin gives the minimum energy boundary, has to be the + * same size as eMax + * \param eMax gives the maximum energy boundary + * \param order is the output order. Defaults to the order of + * inMap + * + * \return a double precision valarray of the expected counts of + * the spectra + */ + std::valarray spectraToCounts(double index, double prefactor, double pivotEnergy, const std::valarray &eMin, const std::valarray &eMax, const SM::Coordinate &co) const; + /** \brief A constant reference to the energy value the exposure is + * evaluated at. + */ + const std::valarray & getEnergies() const; + private: + /** \brief Create an array of indexes suitable for the exposure + * integration. + * + * \return a vector with all the energies per bin + */ + std::vector > indexSort(const std::valarray &eMin, const std::valarray &eMax, const std::valarray &energies, std::vector > &fluxIndex, std::vector > &expIndex) const; + /** \brief Convert a flux density spectra to counts. The input + * should be in units of cm^-2 s^-1 MeV^-1. + * + * \param fluxIndex and \param expIndex should be initialized with + * indexSort + * \param energies is the output from indexSort + * \param flux is the input flux, evaluated at \param enFlux + */ + std::valarray convertSpectra(const std::vector > &energies, const std::valarray &flux, const std::valarray &enFlux, const std::vector > &fluxIndex, std::vector > &expIndex, const SM::Coordinate &co) const; + std::valarray convertSpectra(const std::vector > &energies, const ArraySlice &flux, const std::valarray &enFlux, const std::vector > &fluxIndex, std::vector > &expIndex, const SM::Coordinate &co) const; + /** \brief Calculate Lagrange interpolation */ + void CalculateInterpolation(); + Skymap fDataSkymap; //!< Store the exposure in a skymap + Skymap fPrefactor, fIndex; //!< Store the power-law interpolation indexes +}; +#endif + diff --git a/source/FITS.h b/source/FITS.h new file mode 100644 index 0000000000000000000000000000000000000000..bdc31d5d080372ddc099259d100a6ade47db3408 --- /dev/null +++ b/source/FITS.h @@ -0,0 +1,84 @@ +class FITS { + + public: + + long NAXIS; + long *NAXES; + double *CRVAL,*CDELT,*CRPIX; + + double CHIPOL,PSIPOL,CHIORI,PSIORI; + + + int nelements; + char filename[1000]; + + + FITS() ; + FITS( int NAXIS1); //AWS20061121 + int init(int NAXIS1); //AWS20061121 + FITS( int NAXIS1,int NAXIS2); + int init(int NAXIS1,int NAXIS2); + FITS( int NAXIS1,int NAXIS2,int NAXIS3); + int init(int NAXIS1,int NAXIS2,int NAXIS3); + FITS( int NAXIS1,int NAXIS2,int NAXIS3,int NAXIS4); + int init(int NAXIS1,int NAXIS2,int NAXIS3,int NAXIS4); + + + void setvalues(); + + float *array; + float & operator()(int i); //AWS20061121 + float & operator()(int i,int j); + float & operator()(int i,int j,int k); + float & operator()(int i,int j,int k,int l); + float & operator()(int i,int j,int k,int l,int m); //AWS20130903 + + FITS operator =(double a); + FITS operator+=(double a); + FITS operator-=(double a); //AWS20101223 + FITS operator*=(double a); + FITS operator/=(double a); //AWS20101223 + + FITS operator =(FITS a); //AWS20101223 + FITS operator+=(FITS a); //AWS20101223 + FITS operator-=(FITS a); //AWS20101223 + FITS operator*=(FITS a); //AWS20101223 + FITS operator/=(FITS a); //AWS20101223 + + + void print(); + double sum(); + double max(); + int smooth(int n_smooth); + + int read (char *filename_); + int read (char *filename_,int hdunum);//AWS20040406 + int read (char *filename_,int hdunum1,int hdunum2,int options);//AWS20070928 + + int read(char *filename_,int naxis_min,int naxis_max); + int write(char *filename_); + + int read_keyword_double (char *filename_,int hdunum,char *keyword,double &value);//AWS20050720 + int write_keyword_double(char *filename_,int hdunum,char *keyword,double value,char * comment);//AWS20061122 + int write_keyword_double( int hdunum,char *keyword,double value,char * comment);//AWS20061122 + + private: + float junk; // used to return pointer when bounds violated AWS20061121 + +}; + + + + + + + + + + + + + + + + diff --git a/source/Fermi_Aeff.h b/source/Fermi_Aeff.h new file mode 100644 index 0000000000000000000000000000000000000000..3db6c3679d0ec8b1eeaad3863b194e852242a67b --- /dev/null +++ b/source/Fermi_Aeff.h @@ -0,0 +1,15 @@ +class Fermi_Aeff +{ + public: + int nE_true; + valarray Aeff_average; + valarray E_true; + + void read(string filename, int debug); + void read(string filename, string filetype, int debug); + + double Aeff_average_interpolated(double E, int debug); + + void test(); + int initialized; +}; diff --git a/source/Fermi_systematic_errors.h b/source/Fermi_systematic_errors.h new file mode 100644 index 0000000000000000000000000000000000000000..19b29b1f0fcea2d2c45f79df1b9c553e12a1bd8b --- /dev/null +++ b/source/Fermi_systematic_errors.h @@ -0,0 +1 @@ +double Fermi_systematic_errors(double E, int mode); diff --git a/source/GCR_data.h b/source/GCR_data.h new file mode 100644 index 0000000000000000000000000000000000000000..8bbd22bff56b1c87db2d2ecb93cf2d53aac242bf --- /dev/null +++ b/source/GCR_data.h @@ -0,0 +1,66 @@ +class GCR_data +{ + + public: + + char database_file[400]; + int n; // number of entries + + // as read + char **entry; // complete input string + char *status; // # etc use, dont use etc. + char **reference; + char **experiment; + char **Y_name; // e.g. flux + double *E_low_input; + double *E_high_input; + double *E_mean_input; + double *value_input; + double *err_minus_input; + double *err_plus_input; + + char *err_type; // a=abs r=relative + + int **Z_numerator; + int **A_numerator; + int **Z_denominator; + int **A_denominator; + + char **comment; + char **X_units; + char **Y_units; + + // converted and derived quantities + + char energy_units[4];// "MeV"|"GeV" + char area_units [4]; // "cm2"|"m2" + + double *E_low; // in energy units + double *E_high; + double *E_mean; + + double *value; + double *err_minus; + double *err_plus; + + + int n_ZA_numerator; + int n_ZA_denominator; + + + // information for plotting: user-defined // AWS20060621 + int *color; // AWS20060621 + int *style; // AWS20060621 + double *size ; // AWS20060621 + + int read(char* database_file_, char *area_units, char *energy_units_); + + int read(char *database_file_, char *area_units_, char *energy_units_, + char *Y_name_, + int n_ZA_numerator_select, int *Z_numerator_select, int* A_numerator_select, + int n_ZA_denominator_select,int *Z_denominator_select,int* A_denominator_select); + + int set_plotting(); // AWS20060621 + + int print(); +}; diff --git a/source/GalacticRadiationField.h b/source/GalacticRadiationField.h new file mode 100644 index 0000000000000000000000000000000000000000..a70bb987ed9343b3129fa02ee1d41e3d1a26fc76 --- /dev/null +++ b/source/GalacticRadiationField.h @@ -0,0 +1,74 @@ +#ifndef _rf_GalacticRadiationField_h_ +#define _rf_GalacticRadiationField_h_ + +#include + +#include +//#include +//#include +#include + +namespace rf { + + class GalacticRadiationField { + + public: + + enum Component { STELLAR, SCATTERED, INFRARED, CMB, TOTAL }; + + GalacticRadiationField(const std::string& filename, + bool useIsotropic = false); + ~GalacticRadiationField(); + + // Energy density returned in (eV cm^-3 micron^-1)*micron + + double GetEnergyDensity(const double wl, + const double r, + const double z, + const Component comp) const; + + // Intensity returned in (eV cm^-2 s^-1 sr^-1 micron^-1)*micron + + double GetIntensity(const double wl, + const double r, + const double z, + const double azimuth, + const double cosZenith, + const Component comp) const; + + const std::valarray& GetWavelengthData() const { return fWavelengthData; } + const std::valarray& GetAzimuthData() const { return fAzimuthData; } + const std::valarray& GetCosZenithData() const { return fCosZenithData; } + const std::valarray& GetRData() const { return fRData; } + const std::valarray& GetZData() const { return fZData; } + + private: + + bool fInitialised; + + double fAzimuthDelta, fCosZenithDelta; + + std::valarray fRData, fZData, fWavelengthData, fAzimuthData, fCosZenithData; + + //float* fEnergyDensity; + //float* fAngularDistribution; + + std::valarray< std::valarray< std::valarray > > fEnergyDensityStellar, fEnergyDensityScattered, fEnergyDensityInfrared; + + std::valarray< std::valarray< std::valarray< std::valarray< std::valarray > > > > //fIntensityStellar, fIntensityScattered, fIntensityInfrared, + fIntensityTotal; + + GalacticRadiationField(); + + void ExtractFITSData(const std::string& filename, const bool useIsotropic); + void ExtractFITSEnergyDensity(fitsfile* fptr); + void ExtractFITSAngularDistribution(fitsfile* fptr); + + double BlackBodyEnergyDensity(const double energy, + const double kT) const; + + }; + +} + +#endif diff --git a/source/Galaxy.h b/source/Galaxy.h new file mode 100644 index 0000000000000000000000000000000000000000..5f211c85b86b641e486fb0793b4e285f75eafbe9 --- /dev/null +++ b/source/Galaxy.h @@ -0,0 +1,142 @@ + +//**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****| +// * Galaxy.h * galprop package * 4/14/2000 +//**"****!****"****!****"****!****"****!****"****!****"****!****"****!****"****| + +#ifndef Galaxy_h +#define Galaxy_h + +#include"Distribution.h" +#include +#include "Skymap.h" + +class Galaxy { + + public: + + Galaxy(); + ~Galaxy(); + + double z_min,z_max,dz; // for 1,2,3D + double r_min,r_max,dr; // for 2D + double x_min,x_max,dx,y_min,y_max,dy; // for 3D + + double p_min ,p_max,p_factor; // momentum start, end, factor + + int n_spatial_dimensions;// 1,2,3D + int n_pgrid; // number of points in momentum + int n_rgrid; // number of points in radius (2D) + int n_zgrid; // number of points in z (1D,2D) + int n_xgrid; // number of points in x (3D) + int n_ygrid; // number of points in y (3D) + + double *x; // x grid + double *y; // y grid + double *z; // z grid + double *r; // r grid + + int n_E_gammagrid; // number of points in gamma-ray energy + double E_gamma_min,E_gamma_max,E_gamma_factor;// min,max, factor for gamma-ray energy + double *E_gamma; // gamma ray energy grid + int n_lat,n_long; // dimensions of gamma-ray skymaps + double long_min,long_max; // longitude range of gamma-ray skymaps + double lat_min, lat_max; // latitude range of gamma-ray skymaps + double d_long,d_lat; // longitude, latitude binsize of gamma-ray skymaps + + int n_nu_synchgrid; // number of points in synchrotron frequency + double nu_synch_min,nu_synch_max,nu_synch_factor; // min,max, factor for synchrotron frequency + double *nu_synch; // synchrotron frequency grid + + Distribution n_HI; + Distribution n_H2; + Distribution n_HII; + + Distribution HIR; // HI column density map in Galactocentric rings AWS2001025 + Distribution COR; // CO column density map in Galactocentric rings AWS2001025 + Skymap hpHIR; // HI column density healpix map in Galactocentric rings + Skymap hpCOR; // HI column density healpix map in Galactocentric rings + + Distribution B_field; + Distribution* ISRF; + Distribution* ISRF_energy_density; + int n_ISRF_components; + + std::valarray fISRFFactors;//AWS20100409 + rf::GalacticRadiationField* fISRFModel; + + double* nu_ISRF; // ISRF frequencies + + Distribution SNR_cell_time; // time between SNR for each cell + Distribution SNR_cell_phase; // phase of SNR for each cell + Distribution SNR_electron_dg; // electron injection spectral index delta (Gaussian distributed) AWS20010410 + Distribution SNR_nuc_dg; // nucleus injection spectral index delta (Gaussian distributed) AWS20010410 + + Distribution bremss_emiss; // bremsstrahlung emissivity on neutral gas + Distribution bremss_ionized_emiss; // bremsstrahlung emissivity on ionized gas + Distribution* IC_iso_emiss; // inverse Compton isotropic emissivity for each ISRF component + Distribution* IC_aniso_emiss; + Distribution pi0_decay_emiss; // pi0 decay emissivity + + Distribution bremss_skymap; // bremsstrahlung intensity skymap on neutral gas + Distribution bremss_ionized_skymap; // bremsstrahlung intensity skymap on ionized gas + Distribution* IC_iso_skymap; // inverse Compton isotropic intensity skymap for each ISRF component + Distribution* IC_aniso_skymap; // inverse Compton anisotropic intensity skymap for each ISRF component IMOS20060420 + Distribution pi0_decay_skymap; // pi0 decay intensity skymap + + Distribution bremss_HIR_skymap; // bremsstrahlung intensity skymap on HI in rings AWS20041214 + Distribution bremss_H2R_skymap; // bremsstrahlung intensity skymap on CO in rings AWS20041214 + Distribution pi0_decay_HIR_skymap; // pi0 decay intensity skymap on HI in rings AWS20041214 + Distribution pi0_decay_H2R_skymap; // pi0 decay intensity skymap on H2 in rings AWS20041214 + + Skymap bremss_hp_skymap; // bremsstrahlung intensity skymap on neutral gas + Skymap bremss_ionized_hp_skymap; // bremsstrahlung intensity skymap on ionized gas + Skymap* IC_iso_hp_skymap; // inverse Compton isotropic intensity skymap for each ISRF component + Skymap* IC_aniso_hp_skymap; // inverse Compton anisotropic intensity skymap for each ISRF component IMOS20060420 + Skymap pi0_decay_hp_skymap; // pi0 decay intensity skymap + + Skymap* bremss_HIR_hp_skymap; // bremsstrahlung intensity skymap on HI in rings AWS20041214 + Skymap* bremss_H2R_hp_skymap; // bremsstrahlung intensity skymap on CO in rings AWS20041214 + Skymap* pi0_decay_HIR_hp_skymap; // pi0 decay intensity skymap on HI in rings AWS20041214 + Skymap* pi0_decay_H2R_hp_skymap; // pi0 decay intensity skymap on H2 in rings AWS20041214 + + + Distribution synchrotron_emiss; // synchrotron emissivity + Distribution synchrotron_Q_emiss; // synchrotron Stokes Q emissivity AWS20110407 + Distribution synchrotron_U_emiss; // synchrotron Stokes U emissivity AWS20110407 + + Distribution synchrotron_skymap; // synchrotron Cartesian skymap + Distribution synchrotron_Q_skymap;// synchrotron Stokes Q Cartesian skymap AWS20110407 + Distribution synchrotron_U_skymap;// synchrotron Stokes U Cartesian skymap AWS20110407 + Distribution synchrotron_P_skymap;// synchrotron polarized intensity Cartesian skymap AWS20110407 + + Skymap synchrotron_hp_skymap; + Skymap synchrotron_Q_hp_skymap; // synchrotron Stokes Q healpix skymap AWS20110407 + Skymap synchrotron_U_hp_skymap; // synchrotron Stokes U healpix skymap AWS20110407 + Skymap synchrotron_P_hp_skymap; // synchrotron polarized intensity healpix skymap AWS20110407 + + Distribution free_free_emiss; // free-free emissivity AWS20110906 + Distribution free_free_skymap; // free-free intensity Cartesian skymap AWS20110906 + Skymap free_free_hp_skymap; // free-free intensity healpix skymap AWS20110906 + + + Distribution ionization_rate; // cosmic-ray ionization rate + + Distribution DM_emiss; // DM gamma-ray emission distribution IMOS20050912 + Distribution DM_skymap; // DM gamma-ray emission skymap IMOS20050912 + Skymap DM_hp_skymap; + + //interface functions prototypes + void init( double r_min_, double r_max_, double dr_, + double z_min_, double z_max_, double dz_); + + void init( double x_min_, double x_max_, double dx_, + double y_min_, double y_max_, double dy_, + double z_min_, double z_max_, double dz_); + + void print(); + +}; + +#endif + + diff --git a/source/Galdef.h b/source/Galdef.h new file mode 100644 index 0000000000000000000000000000000000000000..5e287e19f9b853b1eceef135d04e67fe4bafb93d --- /dev/null +++ b/source/Galdef.h @@ -0,0 +1,199 @@ + +//**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****| +// * Galdef.h * galprop package * 10/12/2003 +//**"****!****"****!****"****!****"****!****"****!****"****!****"****!****"****| + +#ifndef Galdef_h +#define Galdef_h + +using namespace std; //AWS20050624 +#include //AWS20050624 +#include //AWS20050624 + +class Galdef { + + public: + + char version[10]; // galprop version number + char run_no [100]; // identifier of galprop run AWS20100326 + char galdef_ID[100]; // full identifier e.g. 01_123456 + + int n_spatial_dimensions; // 1,2 or 3 + double z_min,z_max,dz; // for 1,2,3D + double r_min,r_max,dr; // for 2D + double x_min,x_max,dx,y_min,y_max,dy; // for 3D + double p_min,p_max,p_factor; // momentum start, end, factor + double Ekin_min,Ekin_max,Ekin_factor; // kinetic energy per nucleon start, end, factor + + char p_Ekin_grid[5]; // "p"||"Ekin": construct grid in p or Ekin + + double E_gamma_min,E_gamma_max,E_gamma_factor; // gamma-ray energy (MeV) start, end, factor + double long_min,long_max; // gamma-ray skymap longitude min,max (degrees) + double lat_min, lat_max; // gamma-ray skymap latitude min,max (degrees) + double d_long,d_lat; // gamma-ray skymap longitude,latitude binsize (degrees) + int integration_mode; // integr.over the particle spectrum: =1-old E*logE; !=1-power-law analyt. + int healpix_order; // The order for healpix maps. + + double nu_synch_min,nu_synch_max,nu_synch_factor;// synchrotron frequency min,max (Hz), factor + + double D0_xx; // diffusion coefficient at break rigidity + double D_rigid_br; // break rigidity for diffusion coefficient in MV + double D_g_1; // diffusion coefficient index below break rigidity + double D_g_2; // diffusion coefficient index above break rigidity + + int diff_reacc; // 1,2=incl.diff.reacc.; 11=Kolmogorov+damping; 12=Kraichnan+damping + double v_Alfven; // Alfven speed in km s-1 + double damping_p0; // ~ 1.e5 MV -characteristic rigidity IMOS20030129 + double damping_max_path_L; // ~ 1.*kpc2cm; Lmax~1 kpc, max free path IMOS20030129 + double damping_const_G; // a const derived from fitting B/C IMOS20030129 + + int convection; // 1=include convection AWS20010323 + double v0_conv; // v=v0_conv+dvdz_conv*abs(z) AWS20010323 + double dvdz_conv; // v=v0_conv+dvdz_conv*abs(z) AWS20010323 + + double nuc_rigid_br; // break rigidity for nuclei injection in MV + double nuc_g_1; // spectral index below break + double nuc_g_2; // spectral index above break + + char inj_spectrum_type[9]; // "rigidity"||"beta_rig"||"Etot": choose the nucleon injection spectrum IMOS20000613.1 + + double electron_g_0; // index below electron_rigid_br0 IMOS20031012 + double electron_rigid_br0; // break rigidity0 for electron injection in MV + double electron_g_1; // spectral index between breaks + double electron_rigid_br; // break rigidity1 for electron injection in MV + double electron_g_2; // spectral index above electron_rigid_br + + double He_H_ratio; // He/H of ISM, by number + int n_X_CO; // number of X_CO values for Galactocentric rings AWS20040227 + double *X_CO; // conversion factor CO integrated temperature -> H2 column density AWS20040227 + int fragmentation; // 1=include fragmentation + int momentum_losses; // 1=include momentum losses + int radioactive_decay; // 1=include radioactive_decay + int K_capture; // 1=include K_capture AWS20010731 + int ionization_rate; // 1=compute ionization rate IMOS20060420 + + double start_timestep; // start time step in years + double end_timestep; // end time step in years + double timestep_factor; // factor to multiply timestep + int timestep_repeat; // number of times to repeat for factor + int timestep_repeat2; // number of times to repeat in timestep_mode=2 + int timestep_print ; // number of timesteps between printings + int timestep_diagnostics; // number of timesteps between propel diagnostics + int control_diagnostics; // control details of propel diagnostics + + int network_iterations; // number of iterations of the entire network + + + int prop_r; // for 2D: 1=propagate in r; + int prop_x,prop_y,prop_z; // for 2D: 1=propagate in z;for 3D 1=propagate in x,y,z + int prop_p; // propagate in momentum + int use_symmetry; // xyz symmetry (3D) + + int vectorized; // 0=unvectorized code, 1=vectorized code + + int source_specification; // 0,1,2 + double source_normalization; // 1. IMOS20030129 + double electron_source_normalization; // 1. IMOS20031016 + + int max_Z; // maximum nucleus Z in galdef file + int *use_Z; // 1=use this nucleus Z + + double **isotopic_abundance; // isotopic abundances + + int total_cross_section; // total inel. cross-sections option AWS20010620 + int cross_section_option; // controls which cross-sections used + + double t_half_limit; // lower limit on radioactive half-life for explicit inclusion AWS20010731 + + int primary_electrons; // 1=propagate primary electrons + int secondary_positrons; // 1=propagate secondary positrons + int secondary_electrons; // 1=propagate secondary electrons + int knock_on_electrons; // 1=propagate knock-on electrons IMOS20060504 + int tertiary_antiprotons; // 1=propagate tertiary antiprotons IMOS20000605.13 + int secondary_antiprotons; // 1=propagate secondary antiprotons + int secondary_protons; // 1=propagate secondary protons IMOS20000605.14 + + int gamma_rays; // 1=compute gamma-ray emission + int pi0_decay; // 1 - Dermer 1986 formalism, 2 - Blattnig et al. 2000,PRD 62,094030 IMOS20050216 + int IC_isotropic; // 1=compute isotropic inverse Compton IMOS20060420 + int IC_anisotropic; // 1=compute anisotropic inverse Compton + int bremss; // 1=compute bremsstrahlung IMOS20060420 + int synchrotron; // 1=compute synchrotron emission + + // DM parameters IMOS20050912 + int DM_positrons; // 1=compute positrons from DM + int DM_electrons; // 1=compute electrons from DM + int DM_antiprotons; // 1=compute antiprotons from DM + int DM_gammas; // 1=compute gamma rays from DM + double // user-defined params of DM (double) + DM_double0, DM_double1, DM_double2, DM_double3, DM_double4, + DM_double5, DM_double6, DM_double7, DM_double8, DM_double9; + int // user-defined params of DM (int) + DM_int0, DM_int1, DM_int2, DM_int3, DM_int4, + DM_int5, DM_int6, DM_int7, DM_int8, DM_int9; + + + int source_model; //1= 2= 3= 5= 6=S&Mattox with cutoff + double source_parameters[10]; // for source_model=1 + + int n_cr_sources; // number of pointlike cosmic-ray sources + double *cr_source_x; // source x positions + double *cr_source_y; // source y positions + double *cr_source_z; // source z positions + double *cr_source_w; // source width sigma in kpc + double *cr_source_L; // source luminosity in TBD units + + int SNR_events; // handle stochastic SNR events + double SNR_interval; // time in years between SNRs in 1 kpc^3 + double SNR_livetime; // CR-producing live-time in years of an SNR + double SNR_electron_sdg; // delta electron source index Gaussian sigma + double SNR_nuc_sdg; // delta nucleus source index Gaussian sigma + double SNR_electron_dgpivot; // delta electron source index pivot rigidity (MeV) AWS20010410 + double SNR_nuc_dgpivot; // delta nucleus source index pivot rigidity (MeV) AWS20010410 + + int HI_survey; // HI survey : 8=original 8 rings+high-latitudes, 9= 9 rings all-sky AWS20050913 + int CO_survey; // CO survey : 8=original 8 rings 9= 9 rings all-sky AWS20050913 + + int B_field_model; // >1000=parameterized model + char B_field_name[100]; // 3D B-field model name ("galprop_original" uses B_field_model) AWS20080313 + int n_B_field_parameters; // number of B-field parameters AWS20080313 + double *B_field_parameters; // parameters for 3D B-field models AWS20080313 + + char ISRF_file[100]; // ISRF input file AWS20050301 + int ISRF_filetype; // 0 for old ( //AWS20051124 +#include //AWS20051124 +#include //AWS20051124 + +class Galplotdef +{ + public: + + // from galprop + char version[10]; // galprop version number + char run_no [100]; // identifier of galprop run AWS20100326 + char galdef_ID[100]; // full identifier e.g. 01_123456 + + /* --- not needed, uses Galdef.h : to be removed + int n_spatial_dimensions; // 1,2 or 3 + double z_min,z_max,dz; // for 1,2,3D + double r_min,r_max,dr; // for 2D + double x_min,x_max,dx,y_min,y_max,dy; // for 3D + double p_min,p_max,p_factor; // momentum start, end, factor + double Ekin_min,Ekin_max,Ekin_factor; // kinetic energy per nucleon start, end, factor + + char p_Ekin_grid[5]; // "p"||"Ekin": construct grid in p or Ekin + + ---------- end of not needed */ + + char *psfile_tag; // tag to be included in postscript file name + int screen_output; // 0=no screen output 1=screen output + int output_format; // 0=none 1=ps 2=gif 3=both + + + double E_gamma_min,E_gamma_max,E_gamma_factor; // gamma-ray energy (MeV) start, end, factor + double long_min,long_max; // gamma-ray skymap longitude min,max (degrees) + double lat_min, lat_max; // gamma-ray skymap latitude min,max (degrees) + double d_long,d_lat; // gamma-ray skymap longitude,latitude binsize (degrees) + + double nu_synch_min,nu_synch_max,nu_synch_factor;// synchrotron frequency min,max (Hz), factor + + + int verbose; // verbosity: 0=min 10=max + int test_suite; // run test suit instead of normal run + + + + int gamma_spectrum, gamma_long_profile, gamma_lat_profile; + double long_min1,long_max1; // gamma-ray skymap longitude min,max (degrees) + double long_min2,long_max2; // gamma-ray skymap longitude min,max (degrees) + double long_binsize; // gamma-ray skymap longitude binsize (degrees) + int long_log_scale; // gamma-ray longitude profile scale 0=linear 1=log AWS20050722 + double lat_min1, lat_max1; // gamma-ray skymap latitude min,max (degrees) + double lat_min2, lat_max2; // gamma-ray skymap latitude min,max (degrees) + double lat_binsize; // gamma-ray skymap latitude binsize (degrees) + int lat_log_scale; // gamma-ray latitude profile scale 0=linear 1=log + + double gamma_spectrum_Emin; // minimum energy on spectrum plot + double gamma_spectrum_Emax; // maximum energy on spectrum plot + double gamma_spectrum_Imin; // minimum intensity*E^2 on spectrum plot + double gamma_spectrum_Imax; // maximum intensity*E^2 on spectrum plot + + int gamma_IC_selectcomp; // select IC components 0=no selection, 1=optical, 2=FIR, 3=CMB AWS20060310 + + int sources_EGRET; // 0=sources removed, 1=sources in + + int energies_EGRET; // number of energy ranges for EGRET AWS20031021 + int convolve_EGRET; // 0=no convolution, 1,2,3=convolve, different normalizations + int convolve_EGRET_iEmax; // maximum EGRET energy range to convolve (starting 1) AWS20031017 + double convolve_EGRET_thmax; // maximum angle for convolution (deg) + + double error_bar_EGRET; // EGRET error bar, as fraction + int spectrum_style_EGRET; // 1=segments, 2=at mean energy 3=both + double spectrum_index_EGRET; // spectral index to plot EGRET data: 0.0=automatic, otherwise value e.g. 2.0 -> E^-2.0 AWS20050916 + int spectrum_50GeV_EGRET; // 1= plot EGRET 10-50 GeV spectrum data points + + int spectrum_cut_total; // 1= cut total prediction outside EGRET energy range AWS20040427 + + int isotropic_use; // 1=ignore isotropic, 2=add to prediction, 3=subtract from data + int isotropic_type; // 1=power law, 2=explicit EGRET ranges + + double isotropic_const; // isotropic background = isotropic_const * E^-isotropic_g + double isotropic_g ; // cm^-2 sr^-1 MeV^-1 s^-1 + + int n_E_EGRET; // number of isotropic background values; set in Galplotdef.cc + double *isotropic_EGRET; // isotropic background in EGRET ranges, cm^-2 s^-1 + + int isotropic_sree_plot; //1=plot Sreekumar et al. isotropic background separately + + int convolve_GLAST; // 0=no convolution, 1=convolve AWS20080515 + int energy_disp_FermiLAT; // 0=no energy dispersion, 1=energy dispersion AWS20140825 + double energy_disp_factor; // factor for energy dispersion interpolation AWS20140829 + + int fit_EGRET ; // 1=fit EGRET data + int fit_nrebin_EGRET ; // number of 0.5*0.5 bins for EGRET fitting rebinning + int fit_counts_min_EGRET ; // minimum counts for EGRET fitting after rebinning + int fit_function_EGRET ; // 1=background only, 2=scaling+background 3=both + int fit_options_EGRET ; // other options + + + int data_EGRET ; // plot EGRET data + + int data_GLAST ; // plot GLAST data AWS20080507 + char *GLAST_counts_file; // GLAST counts healpix data file AWS20080819 + char *GLAST_exposure_file;// GLAST exposure healpix data file AWS20080819 + char *GLAST_psf_file; // GLAST psf data file AWS20080819 + + char *Fermi_edisp_file; // Fermi energy dispersion parameter file AWS20150226 + char *Fermi_edisp_file_type; // Fermi energy dispersion parameter file type AWS20150226 + int Fermi_edisp_use_Aeff; // Fermi energy dispersion use Aeff factor AWS20150226 + double Fermi_edisp_E_true_threshold;// Fermi energy dispersion true energy threshold AWS20150226 + double Fermi_edisp_E_meas_threshold;// Fermi energy dispersion measured energy threshold AWS20150226 + + char *isotropic_bgd_file; // isotropic background data file AWS20081119 + char *Fermi_cat_file ; // Fermi catalogue FITS file AWS20081215 + char *Fermi_cat_file2 ; // Fermi catalogue FITS file AWS20140114 + char *Fermi_cat_file3 ; // Fermi catalogue FITS file AWS20140115 + + char *Fermi_sensitivity_file; // Fermi sensitivity skymap file AWS20170112 + + int data_COMPTEL ; // plot COMPTEL data + int data_OSSE ; // plot OSSE data + int data_RXTE ; // plot RXTE data + int data_INTEGRAL ; // plot INTEGRAL/SPI data + int*data_INTEGRAL_comp ; // INTEGRAL/SPI range of fitted components to plot + int data_INTEGRAL_mode; // INTEGRAL/SPI 1=plot data 2=fit 3=both + double data_INTEGRAL_syserr;// INTEGRAL/SPI systematic error AWS20050412 + double data_INTEGRAL_sysemn;// INTEGRAL/SPI systematic error min. energy to apply AWS20050412 + int data_IBIS ; // plot IBIS data AWS20041018 + int data_MILAGRO ; // plot MILAGRO data + int data_HEGRA ; // plot HEGRA data + int data_WHIPPLE ; // plot WHIPPLE data + int data_GINGA ; // plot GINGA data + int data_Chandra ; // plot Chandra data AWS20050909 + + int spiskymax_profile ; // plot SPI spiskymax profile AWS20050722 + char *spiskymax_image_ID; // spiskymax_image.xxxx.fits: identify spiskymax image AWS20050722 + int spiskymax_iteration; // spiskymax iteration selection AWS20050722 + + + int model_ridge ; // plot ridge model + + int sourcepop1_verbose ; // source population parameters AWS20050911 + double sourcepop1_density0 ; + double sourcepop1_oversample;// source oversampling to avoid statistical fluctuations AWS20060109 + double sourcepop1_L_min ; + double sourcepop1_L_max ; + double sourcepop1_alpha_L ; + double sourcepop1_fluxlimit; + double sourcepop1_alpha_R ; + double sourcepop1_beta_R ; + double sourcepop1_zscale ; + double sourcepop1_alpha_z ; + + int sourcepop1_spectral_model; // AWS20100824 + + double sourcepop1_spectrumg ; + double sourcepop1_spect_g_0 ;// AWS20051219 + double sourcepop1_spect_br0 ; + double sourcepop1_spect_g_1 ; + double sourcepop1_spect_br1 ; + double sourcepop1_spect_g_2 ; + double sourcepop1_E_ref_low ;// AWS20060109 + double sourcepop1_E_ref_high;// AWS20060109 + + double sourcepop1_spectrum_g_0_sigma ;// spectrum with two breaks: Gaussian distribution AWS20100824 + double sourcepop1_spectrum_br0_sigma ; + double sourcepop1_spectrum_g_1_sigma ; + double sourcepop1_spectrum_br1_sigma ; + double sourcepop1_spectrum_g_2_sigma ; + + double sourcepop1_specnormE; + + char *sourcepop_SPI_cat; // AWS20060111 + + int sourcepop_total; // plot intensity of sources total + int sourcepop_sublimit; // plot intensity of sources below limit + int sourcepop_soplimit; // plot intensity of sources above limit + + int sourcepop_ext_model; // external model + + + + int gcr_spectra ; + double gcr_spectra_r_min ; + double gcr_spectra_r_max ; + double gcr_spectra_z_min ; + double gcr_spectra_z_max ; + + double gcr_spectra_Emin; // minimum energy on spectrum plot + double gcr_spectra_Emax; // maximum energy on spectrum plot + double gcr_spectra_Imin; // minimum intensity*E^2 on spectrum plot + double gcr_spectra_Imax; // maximum intensity*E^2 on spectrum plot + + int gcr_spectra_n_ZA ; // implicit in input string, not read + int *gcr_spectra_Z ; + int *gcr_spectra_A ; + + int gcr_spectra_n_ratio ; // number of sec/prim ratios + int *gcr_spectra_n_sec_ZA ; // numbers of Z,A in secondary lists (implicit, not read) + int *gcr_spectra_n_pri_ZA ; // numbers of Z,A in primary lists (implicit, not read) + int **gcr_spectra_sec_Z ; // lists of sec Z for sec/prim ratios + int **gcr_spectra_sec_A ; // lists of sec A for sec/prim ratios + int **gcr_spectra_pri_Z ; // lists of pri Z for sec/prim ratios + int **gcr_spectra_pri_A ; // lists of pri A for sec/prim ratios +double *gcr_spectra_ratio_min ; // minimum value for sec/prim ratios +double *gcr_spectra_ratio_max ; // maximum value for sec/prim ratios + + int gcr_spectra_n_mod_phi ; // number of modulation phi (implicit, not read) +double *gcr_spectra_mod_phi ; // list of modulation phi (MV) + + char *gcr_database_file ; // filename for database file AWS20090220 + + double *sync_index_nu ; // frequencies (MHz) in pairs for synchrotron index AWS20070523 + int sync_index_n_nu ; // number of pairs for synchrotron index AWS20070523 + + int sync_data____10MHz ; // plot synchrotron data at this frequency AWS20071221 + int sync_data____22MHz ; // plot synchrotron data at this frequency AWS20071221 + int sync_data____45MHz ; // plot synchrotron data at this frequency AWS20071221 + int sync_data___150MHz ; // plot synchrotron data at this frequency AWS20071221 + int sync_data___408MHz ; // plot synchrotron data at this frequency AWS20071221 + int sync_data___820MHz ; // plot synchrotron data at this frequency AWS20071221 + int sync_data__1420MHz ; // plot synchrotron data at this frequency AWS20071221 + int sync_data__2326MHz ; // plot synchrotron data at this frequency AWS20071221 + int sync_data_22800MHz ; // plot synchrotron data at this frequency AWS20071221 + int sync_data_33000MHz ; // plot synchrotron data at this frequency AWS20071221 + int sync_data_41000MHz ; // plot synchrotron data at this frequency AWS20071221 + int sync_data_61000MHz ; // plot synchrotron data at this frequency AWS20071221 + int sync_data_94000MHz ; // plot synchrotron data at this frequency AWS20071221 + + int sync_data_WMAP ; // select WMAP data AWS20071221 + int sync_data_options ; // more synch data options AWS20120110 + + int free_free_options ; // 1=model 2=WMAP MEM dust template AWS20111222 + int spin_dust_options ; // 1=WMAP MCMC dust+spinning dust 2=WMAP MEM dust template AWS20120117 + + int galdef_series ; // to plot several runs in plot_luminosity_multiple_galdef AWS20100329 + + int luminosity ; // 1=compute luminosity AWS20100823 + + double ISRF_luminosity_R ; // Galactocentric radius for ISRF luminosity calculation AWS20100421 + double ISRF_luminosity_z ; // Galactocentric height for ISRF luminosity calculation AWS20100421 + + int ISRF_read ; // read ISRF AWS20100818 + + +//interface functions prototypes + int read(char *version_,char *run_no_,char *galdef_directory); + int read_galdef_parameter(char *filename,char *parstring,double *value); + int read_galdef_parameter(char *filename,char *parstring,int *value); + int read_galdef_parameter(char *filename,char *parstring,char *value); + void print(); +}; + +#endif + + + + + + + + + + diff --git a/source/Goodness.h b/source/Goodness.h new file mode 100644 index 0000000000000000000000000000000000000000..6c04bc11f1c8bfedc69fd6490fcec1efbcd9def5 --- /dev/null +++ b/source/Goodness.h @@ -0,0 +1,80 @@ +#ifndef Goodness_h +#define Goodness_h + +#include "Counts.h" +#include "Skymap.h" +#include "Model.h" +#include "Variables.h" +#include +#include + +/**\brief Base class for all goodness of fit evaluator. + * + * This is a virtual class which has to be extended. Defines basic routines + * which have to be implemented by extensions. + */ +class BaseGoodness { + public: + /** \brief Construct a goodness object with counts, a model and a filter + * map. + * + * \param model a reference to the model to fit + * \param filter is a skymap defining areas to exclude from the fit. Every + * pixel with a value of 0 in the filter is excluded. It has to have the + * same size as the counts map used when constructing the model. + */ + BaseGoodness(BaseModel &model, const Skymap & filter) : fcounts(model.getCounts()), fmodel(model), ffilter(filter){} + /** \brief The default operator evaluates the goodness of fit. + * + * \param variables is a Variables object defining the values of the + * model variables + * \return the goodness of fit + */ + virtual double operator () (const Variables &variables) const = 0; + /** \brief Calculate the gradient + * + * \param variables is a Variables object defining the values of the model + * parameters + * \return a map of string double pairs, where the string contains the name + * of the variable while the double holds the value of the gradient for + * that variable. + */ + virtual std::map gradient(const Variables &variables) const = 0; + /** \brief Calculate the Hessian matrix. + * + * \param variables is a Variables object defining the values of the model + * parameters + * \return a double map where the first and second index is the two + * variables the second derivative was evaluated for. Their order does not + * matter, the value is assumed to be the same, independent of order. + */ + virtual std::map > hessian(const Variables &variables) const = 0; + /** \brief Return a Skymap with the goodness of each pixel. + * + * \param vars is a Variables object defining the model variable values. + * \return a double precision Skymap of the goodness. + */ + virtual Skymap getMap (const Variables & vars) const = 0; + /** \brief Return a constant reference to the model */ + const BaseModel & getModel() const {return fmodel;} + protected: + const CountsMap &fcounts; //!< Constant reference to the counts + const Skymap ffilter; //!< The filter Skymap, 0 means exclude, everything else include + BaseModel &fmodel; //!< The reference to the model is not constant, since we want the models to be able to cache their results for faster computations +}; + +/**\brief Uses log likelihood to evaluate the goodness of fit + */ +class LogLikelihood : public BaseGoodness { + public: + /** Constructor to initialize the base goodness */ + LogLikelihood(BaseModel & model, const Skymap & filter); + double operator () (const Variables & variables) const; + std::map gradient(const Variables & variables) const; + std::map > hessian(const Variables &variables) const; + Skymap getMap (const Variables & vars) const; + private: + std::vector flogSumInt; //!< Pre-calculate the log of sum of factorials to get the proper likelihood without much CPU cost +}; + +#endif diff --git a/source/HIH2IC.h b/source/HIH2IC.h new file mode 100644 index 0000000000000000000000000000000000000000..42c36fbe4ae28bdb53e5d8695dd8db92a63625e1 --- /dev/null +++ b/source/HIH2IC.h @@ -0,0 +1,44 @@ +#ifndef HIH2IC_h +#define HIH2IC_h + +#include "Model.h" +#include "Exposure.h" +#include "Counts.h" +#include "Psf.h" +#include "Sources.h" +#include "Parameters.h" +#include "Variables.h" + +/**\brief A sample model class for people to modify. + * + * This model is not intended to be used unmodified. + */ +class HIH2IC : public BaseModel { + public: + /** \brief The constructor takes the same parameters as the BaseModel constructor. + * + * The constructor should load and prepare the data used in the model. It is advised to cache + * the results between computation to speed up model fitting, but please keep memory consumption minimal. + * + * This model loads a skymap from a fits file and scales it. It also adds a Gaussian skymap profile with + * a power law energy dependance. + */ + + // use bool to avoid ambiguity between cases + + HIH2IC(const Exposure &exposure, const CountsMap &counts, const Psf &psf, const Sources &sources, const Parameters &pars, const Skymap &filter, unsigned int configure = 3); //AWS20080717 + HIH2IC(const Exposure &exposure, const CountsMap &counts, const Psf &psf, const Sources &sources, const Parameters &pars, const Skymap &filter, int nrings_in, unsigned int configure = 3); //AWS20080623 AWS20080717 + + void getMap(const Variables &vars, Skymap &map); + gradMap getComponents(const Variables & vars, const std::string &prefix); + void getGrads(const Variables &vars, const std::string & varName, Skymap &map); + gradMap getSecDer(const Variables & vars, const std::string & varName1, const std::string & varName2); + private: + + Skymap fpi0Map,fbremssMap,ficsMap; //AWS20080417 + Skymap fHIMap, fH2Map; //AWS20080421 + Skymap fHIIMap; //AWS20091217 + double flastPrefactor, flastIndex; +}; + +#endif diff --git a/source/HealpixBaseExtended.h b/source/HealpixBaseExtended.h new file mode 100644 index 0000000000000000000000000000000000000000..6ff5ae2537cca004aad5c89f6d7f89491a7e26e9 --- /dev/null +++ b/source/HealpixBaseExtended.h @@ -0,0 +1,78 @@ +/** \class HealpixBaseExtended + * \brief An extension for the Healpix_Base class with local modifications + * + * A place to put local additions to the Healpix_Base class that are necessary + * for the Skymap but not dependent on the template type + */ +#ifndef HealpixBaseExtended_h +#define HealpixBaseExtended_h + +#include "healpix_base.h" +#include "Coordinate.h" +#include "PhysicalConstants.h" +#include "Region.h" +#include +#include +#include + +class HealpixBaseExtended : public Healpix_Base { + private: + /**\brief A generator to insert simple ranges into containers */ + class counter { + public: + counter(int init=0) : n(init) {} + int operator () () { return ++n; } + private: + int n; + }; + + public: + /**\brief Return the resolution in degrees. + * + * This is the square root of the solid angle. + */ + double resolution() const{return sqrt(3/utl::kPi)*60/nside_;} + /**\brief Solid angle of each pixel in radians */ + double solidAngle() const{return utl::kPi/(3*nside_*nside_);} + //! Coordinate for a given pixel + SM::Coordinate pix2coord(int pix) const { return SM::Coordinate(pix2ang(pix)); } + //! Pixel for a given coordinate + int coord2pix(const SM::Coordinate & coord) const { return ang2pix(coord.healpixAng()); } + + + /** \brief Return the pixels for a given region of the sky + * + * \parameters region is the selected region + * + * \returns a set of pixels + * + * All pixels that have their centers within the region are selected. + */ + std::set regionToPixels(const SkyRegion & region) const; + + /**\brief Return the pixel numbers within a rectangular region + * + * \param pointingll is the lower left corner of the rectangle + * \param pointingur is the upper right corner of the rectangle + * \param listpix is the list of pixels who's center lies within the region + * + * This routine assumes a doughnut shaped world, so if the lower limit + * of the rectangle is above (north, theta=0, is up here) the upper limit, + * the pixels within the region above the lower limit and below the upper + * limit is selected. This also applies to the left and right region + * which is not as unusual. + */ + void query_rectangle(const pointing &pointingll, const pointing &pointingur, std::vector &listpix) const; + + /** \brief Return the pixels under a given pixel in a different + * HealpixBase + * + * \parameter hp is the other HEALPix base + * \parameter pixel is the pixel number in the other base + * \param listpix is the list of pixels who's center lies within the region + * + */ + void query_pixel(const Healpix_Base &hp, int pixel, std::vector &listpix) const; +}; + +#endif diff --git a/source/IRF.h b/source/IRF.h new file mode 100644 index 0000000000000000000000000000000000000000..438140c45d0b4608d76f3bf8298af1f0fcf843a7 --- /dev/null +++ b/source/IRF.h @@ -0,0 +1,20 @@ + +//#include"FITS.h" +#include"galprop_classes.h" // since FITS already here and included in galplot + +class IRF +{ + public: + + FITS matrix; + + long n_E; + double *E_lo,*E_hi,*dE; // not yet assigned, taken from RMF + + double *sum,*sum1,*sum2,*sum3;// sums for 3 IRF types as function of energy; + + + + int read (char *filename,int idet1,int idet2); + void print(); +}; diff --git a/source/Isotropic.h b/source/Isotropic.h new file mode 100644 index 0000000000000000000000000000000000000000..77b24c52793ddb24bf4c1c7947482cd5303878a0 --- /dev/null +++ b/source/Isotropic.h @@ -0,0 +1,21 @@ + + + +class Isotropic +{ + public: + + int n_energy; + + double *energy, *spectrum; + + + int read (char *directory, char*filename); + void print(); + + double interpolated(double energy_,int debug); + + double integrated (double energy_1, double energy_2, int debug); //AWS20081202 + + +}; diff --git a/source/Model.h b/source/Model.h new file mode 100644 index 0000000000000000000000000000000000000000..34097016024185f7039fb7bc8321c2ca7f2c9d81 --- /dev/null +++ b/source/Model.h @@ -0,0 +1,433 @@ +#ifndef Model_h +#define Model_h + +#include "Exposure.h" +#include "Counts.h" +#include "Psf.h" +#include "Sources.h" +#include "Variables.h" +#include "Skymap.h" +#include "SparseSkymap.h" +#include "Parameters.h" +#include +#include +#include + +/** \brief A virtual base class for all the models. + * + * This class defines the basic interface for all the models, as well as + * provide some general purpose functions for exposure correction and + * convolution. + */ +class BaseModel { + public: + /** \brief Defines the type returned by the gradient and second derivative + * calculations of the models + */ + typedef std::map > gradMap; + + /** \brief The only constructor available. + * + * \param counts is a reference to a counts map, needed as a reference map + * when doing exposure corrections. + * \param exposure is a reference to an exposure object, needed for + * exposure corrections. + * \param psf is a point spread function, needed for convolution + * \param sources is a reference for source objects. Only + * really used by the SourceModel, but necessary to store here + * for the implementation of source retrieval in the ArrayModel + * \param pars is a reference to a Parameters object, so models can easily be + * configured. + * \param filter is the filter applied to the fit. Every pixel + * where filter == 0 is excluded. + * \param configure is a bitmask controlling the model + * behavior. There are two bits defined in the base model that + * all model should obey + * - EXPOSURE_CORRECT + * - CONVOLVE + * By default, both of them are set. These modify the behavior + * of the built in convolution and exposure correction routine + * in BaseModel. The value of these two are 1 and 2 + * respectively. Other models can use this variable, but + * cannot use these two values for the bitmask. + */ + BaseModel(const CountsMap & counts, const Exposure &exposure, const Psf &psf, const Sources &sources, const Parameters &pars, const Skymap &filter, unsigned int configure = 3); + + const unsigned int EXPOSURE_CORRECT, CONVOLVE; + + /** \brief A virtual function returning the model count skymap given variable values. + * + * \param vars defines the variable values of the model. + * \param map is a double precision Skymap which is of the same + * size as counts map. The results should be added to this + * map. + * + * It is advised that models cache their last values, since the parameters + * need not change between function calls. This will speed up minimization + * considerably. + * + * This method must be implemented in a child object. + */ + virtual void getMap(const Variables & vars, Skymap &map) = 0; + + /** \brief A virtual function returning the different components count skymaps given variable values. + * + * This function can also be used to output files required by + * the model, but please use the prefix for all filenames. + * + * \param vars defines the variable values of the model. + * \return a map of strings and skymaps, where the + * strings are the name of the component and the skymap + * is its counts map + * \param prefix is a string that should be prefixed to all + * output filenames by the model. If no files are written in + * this method, it can be safely ignored. + * + * This method must be implemented in a child object. + */ + virtual gradMap getComponents(const Variables & vars, const std::string &prefix) = 0; + + /** \brief A virtual function returning a skymap of the first derivatives + * of the model with respect to the given variable name + * + * \param vars defines the variable values of the model. + * \param varName is the name of the variable which derivative map should be + * returned + * \param map is a double precision Skymap which is of the same + * size as counts map. The results should be added to this + * map. + * + * It is advised that models cache their last values, since the parameters + * need not change between function calls. This will speed up minimization + * considerably. + * + * This method must be implemented in a child object. + */ + virtual void getGrads(const Variables & vars, const std::string & varName, Skymap &map) = 0; + + /** \brief A virtual function returning a skymap of the second derivatives + * of the model with respect to the given variable names + * + * \param vars defines the variable values of the model. + * \param varName1 and + * \param varName2 are the names of the variables which derivative map should be + * returned + * \return a map with a single varName1 and Skymap pair. If the derivative is 0 in all + * pixels, return an empty map. + * + * It is advised that models cache their last values, since the parameters + * need not change between function calls. This will speed up minimization + * considerably. + * + * This method must be implemented in a child object. + */ + virtual gradMap getSecDer(const Variables & vars, const std::string & varName1, const std::string & varName2) = 0; + + /** \brief Return a reference to the variables */ + Variables & getVariables(); + + /** \brief Return a constant reference to the variables */ + const Variables & getVariables() const; + + /** \brief Set the values of the models of the variables from a variables object + * + * Only changes the value and error of the variables currently available in + * the model, does not set new variables. Throws an error if one of the + * variables in the model is not defined in the new Variables object. + */ + void setVariables(const Variables &vars); + + /** \brief Return a constant reference to the CountsMap object */ + const CountsMap & getCounts() const; + + //! Return a constant reference to the filter skymap + const Skymap & getFilter() const; + + //! Return a copy of the source object with model modifications + Sources getSources(const Variables & vars) const; + + /** \brief Modify the sources object before return. + * + * Mechanism needed for models doing spectral fitting to get the + * information back into the source catalog. Definately not optimal, but + * works. + * + * This function needs to be public for the JoinModel to work properly. This method + * depends on the sources to be the same as the model was initialized with. getSources takes + * care of that. + */ + virtual void modifySources(const Variables &, Sources &) const{}; + + /** \brief Models must be able to throw errors if an exception occurs */ + class ModelError : public std::exception { + private: + std::string eMessage; + public: + /** \brief Construct an error with the default message "Error in Model class" */ + ModelError() : eMessage("Error in Model class"){} + /** \brief Construct an error with a specialized error message + * + * \param message is the specialized message + */ + ModelError(const std::string & message) : eMessage(message){} + ~ModelError() throw(){} + /** \brief Return the error message */ + virtual const char* what () const throw(){ + return(eMessage.c_str()); + } + }; + protected: + /** \brief Converts flux maps to counts. + * + * \param fluxMap is the skymap to be converted to counts. This should be + * in units of 1/Energy/cm^2/s/sr, where energy is the units of energy in the + * skymap's spectra. Note that the energy unit of the flux map and counts + * map have to conform. The energy dependence of the flux is approximated + * with a power-law between spectral points. Finer energy binning in the + * input flux map should be used if this is a bad approximation. + * + * \param psfRebin is a boolean that determines if the return + * map should be binned more finely to get a better accuracy + * for the convolution. + * + * \return a double precision Skymap with the same size and energy + * structure as the counts map unless psfRebin is given, where + * the energy bin structure follows that of the psf. + * + * \note This does not take into account the solid angle of the pixels. + */ + Skymap fluxToCounts(const Skymap &fluxMap, bool psfRebin=false) const; + /** \brief Convolve a skymap with the stored psf + * + * \param skyMap is the input skymap to be convolved. It should be in RING + * scheme to make the computation faster. Output from fluxMap are in that + * scheme. + * + * \param psfRebin is a boolean that determines if the input + * map is more finely binned to get a better accuracy + * for the convolution. The output map will then be rebinned + * to the same binning as the counts map. + * + * \return the convolved skymap, which has the same structure as the + * original one. + * + * Uses the fast spherical harmonics decomposition of healpix to do the + * convolution. Assumes a homogenous spherically symmetric psf, since that + * is the only thing that can be handled by this method. In the case of + * GLAST, this should not matter much. + */ + Skymap convolve(const Skymap & skyMap, bool psfRebin=false) const; + /** \brief Exposure correct a powerlaw spectra, given a prefactor, index + * and a pivot point. + * + * The formula used for the flux is \f$ F(E) = P_0 (E/E_0)^{\\gamma} \f$, + * where \f$ P_0 \f$ is the prefactor, \f$ E_0 \f$ is the pivot point and + * \f$ \\gamma \f$ is the index. + * + * \param prefactor \f$ P_0 \f$ + * \param index \f$ \\gamma \f$ + * \param pivot \f$ E_0 \f$ + * \param co the coordinate we wish to evaluate the exposure + * + * This function does not take the PSF into account, i.e. we assume the + * exposure is uniform enough that it does not matter. This is not true at + * low energies. + */ + std::valarray powerLawCounts(double prefactor, double index, double pivot, const SM::Coordinate &co, bool psfRebin=false, std::valarray *eMin=0, std::valarray *eMax=0) const; + /** \brief Convert a spectra (units of cm^-2 s^-1 + * MeV^-1, given that exposure energy is MeV) at a given + * coordinate to counts. + * + * This method uses a semi-analytical integration, using a + * power law interpolation on the input spectra. + * + * \param intensities is the intensity spectra in units of cm^-2 s^-1 MeV^-1, given that exposure energy is MeV + * \param energies is the energies corresponding to the + * intensities. Has to be the same size as intensities. + * \param co is the coordinate of the flux location + * \param eMinArr gives the minimum energy boundary on return + * \param eMaxArr gives the maximum energy boundary on return + * + * \return a double precision valarray of the expected counts of + * the spectra + */ + std::valarray spectraToCounts(std::valarray &intensities, const std::valarray &energies, const SM::Coordinate &co, bool psfRebin=false, std::valarray *eMinArr=0, std::valarray *eMaxArr=0) const; + const Exposure &fexposure; //!< constant reference to the Exposure object + const CountsMap &fcounts; //!< constant reference to the CountsMap object + const Psf &fpsf; //!< constant reference to the Psf object + const Sources &fsources; //!< constant reference to the Sources object + const Parameters &fparameters; //!< constant reference to the Parameters object + const Skymap &ffilter; //!< constant reference to a filter object. 0 means exclude, 1 include. Models should take that into account if needed. + Variables fvariables; //!< The variables of the model + + unsigned int fconfigure; //!< a bitmask controlling the models behavior + + //Return energy boundaries, compatible with the counts map, but + //as finely binned as the psf. + void getPsfEnergyBoundaries(std::valarray &eMin, std::valarray &eMax) const; + //Rebin a map generated by psf boundaries given above to count + //map boundaries. Does not check spatial binning. + void rebinPsfBoundariesToCounts(Skymap & inMap) const; +}; +/**\brief a model to join other models in a combined one + * + * Stores a vector of pointers to models. Sums up all the skymaps from the models + * as a total model map. Handles all the functions from basemodel gracefully + * so all models are included. + */ +class ArrayModel : public BaseModel{ + public: + ArrayModel(const Exposure &exposure, const CountsMap &counts, const Psf &psf, const Sources &sources, const Parameters &pars, const Skymap &filter, unsigned int configure = 3); + void getMap(const Variables &vars, Skymap&map); + gradMap getComponents(const Variables & vars, const std::string &prefix); + void getGrads(const Variables &vars, const std::string & varName, Skymap&map); + gradMap getSecDer(const Variables & vars, const std::string & varName1, const std::string & varName2); + void setModels(const std::vector &); + protected: + void modifySources(const Variables &vars, Sources &sources) const; + std::vector fmodels; //!< A vector of model pointers to store the used models +}; + + +/**\brief A model to join other models in a combined one. + * + * The constructor has a list of models it accepts. To extend this object, + * just add your model to the list in the constructor. + * + * The models are selected with the parameter + * - modelNames: a space separated list of model names. + */ +class JoinModel : public ArrayModel { + public: + /**\brief The constructor takes the same parameters as the BaseModel constructor. + * + * Initializes the BaseModel as well as all the models given in the modelNames parameter. + */ + JoinModel(const Exposure &exposure, const CountsMap &counts, const Psf &psf, const Sources &sources, const Parameters &pars, const Skymap &filter, unsigned int configure = 3); + ~JoinModel(); +}; + +/**\brief A spatially constant power law model for the Extragalactic + * background. + * + * This model has two variables + * - EGBPrefactor + * - EGBSlope + * + * Uses the BaseModel::powerLawCounts function, so no convolution done. + */ +class EGBModel : public BaseModel { + public: + /** \brief The constructor takes the same parameters as the BaseModel constructor. + * + * Initializes the BaseModel and calculates the EGB skymap from the initial values of the model variables. + * Subsequent calls to the model do not perform exposure correction, but rather just scale the cached + * map to fit the new variables. There is a slight error involved, but it should be far less than the systematic + * errors of the effective area. Refitting with better initial values minimizes this error. + */ + EGBModel(const Exposure &exposure, const CountsMap &counts, const Psf &psf, const Sources &sources, const Parameters &pars, const Skymap &filter, unsigned int configure = 3); + void getMap(const Variables &vars, Skymap &map); + gradMap getComponents(const Variables & vars, const std::string &prefix); + void getGrads(const Variables &vars, const std::string & varName, Skymap &map); + gradMap getSecDer(const Variables & vars, const std::string & varName1, const std::string & varName2); + private: + double flastPrefactor, flastSlope, fcompareSlope; //!< Cache the previous variable values. If the final values are very different from the initial values, do the exposure correction again. + Skymap fcacheMap; //!< The cached skymap + void exposureCorrectCacheMap(const Variables &variables); //!< A helper routine to calculate the skymap using accurate exposure correction + const std::string fprefName, findName; //!< The names of the variables of the model +}; + +/**\brief A model which handles the point sources from the catalog. + * + * Currently it sorts the sources by their flux >100MeV. + * The only available spectral model for the sources is a power law. + * This model is controlled by 3 parameters: + * - numberOfSources: maximum number of sources to include in the map + * - numberOfSourcesPrefactor: number of sources where the prefactor is fitted + * - numberOfSourcesIndex: number of sources where the index is fitted + * + * Number of sources never exceed the value given by numberOfSources. If the + * number of fitted sources in the other cases exceeds that, then they are + * truncated. In all cases, only the brightest sources are selected, so it is + * not yet possible to select which sources to fit. + */ +class SourceModel : public BaseModel { + public: + /** \brief The constructor takes the same parameters as the BaseModel constructor. + * + * Creates the fixed source maps, as well as the necessary sparse skymaps for sources that need to be fitted. + */ + SourceModel(const Exposure &exposure, const CountsMap &counts, const Psf &psf, const Sources &sources, const Parameters &pars, const Skymap &filter, unsigned int configure = 3); + void getMap(const Variables &vars, Skymap &map); + gradMap getComponents(const Variables & vars, const std::string &prefix); + void getGrads(const Variables &vars, const std::string & varName, Skymap &map); + gradMap getSecDer(const Variables & vars, const std::string & varName1, const std::string & varName2); + /** \brief Enum for defining the type of variable. + * + * - Pre: pre factor correction + * - Ind: index correction + * - Unk: unknown parameter type + */ + enum varType { Pre, Ind, Unk }; + private: + /** \brief Add the number i as a string to the prefix. + * + * \param i the number to append + * \param prefix is the prefix + * \return the joined string + */ + std::string prependNumber (int i, const std::string &prefix) const; + /** \brief Given a source number i, find the source corrections. + * + * \param i the source number + * \param prefCorr is the prefix correction + * \param indCorr is the index correction + */ + void getCorr(int i, double &prefCorr, double &indCorr, const Variables & variables) const; + /** \brief Find the type of variables. + * + * \param varName is the name of the variable + * \return Ind or Pre for index and prefactor, respectively, or Unk if the parameter is unknown. + */ + varType getType(const std::string & varName) const; + /** Find the index of a variable name. + * + * \param varName is the name of the variable, and it is assumed it is of a known type + * \return the index for the parameter or -1 if an error occurs. + */ + int getIndex(const std::string & varName) const; + Skymap ffixedSkymap; // > fvariableSkymaps; // fprevIndCorr, fprevPrefCorr; //!< We cache the last known values for the corrections to speed up calculations when they don't change +}; + +/**\brief Basic fitting class that varies the electron and proton + * normalization. + * + * Uses healpix GALPROP output and scales the electron and proton normalization parameters. The variables of the model are + * + * - ElectronNormalization + * - ProtonNormalization + * + * Additionally, one also has to point the model to the GALPROP output files with the parameters + * - galdefID: The galdef ID + * - galpropFitsDirectory: The directory which holds the GALPROP output + */ +class ElectronProton : public BaseModel { + public: + /** \brief The constructor takes the same parameters as the BaseModel constructor. + * + * Loads the skymaps from GALPROP, exposure corrects them, convolves them and adds the inverse Compton and bremsstrahlung + * components together for electron normalization. + */ + ElectronProton(const Exposure &exposure, const CountsMap &counts, const Psf &psf, const Sources &sources, const Parameters &pars, const Skymap &filter, unsigned int configure = 3); + void getMap(const Variables &vars, Skymap &map); + gradMap getComponents(const Variables & vars, const std::string &prefix); + void getGrads(const Variables &vars, const std::string & varName, Skymap &map); + gradMap getSecDer(const Variables & vars, const std::string & varName1, const std::string & varName2); + private: + Skymap felectronSkymap, fpionSkymap; +}; + +#endif diff --git a/source/Parameters.h b/source/Parameters.h new file mode 100644 index 0000000000000000000000000000000000000000..df1e39b2519b45af9346dfb74146f591dd302f05 --- /dev/null +++ b/source/Parameters.h @@ -0,0 +1,241 @@ +#ifndef Parameters_h +#define Parameters_h + +#include +#include +#include +#include +#include +#include + +/** \brief Trim white space from end of string */ +std::string TrimWhitespace(std::string & str); + +/** \brief Class to handle parsing of simple parameter = value pairs. + * + * Uses templated functions and the >> input operator to handle different types + * of values. + * + * For now, parameter can not contain whitespace and there must be a whitespace + * between the parameter and the = sign. + */ +class Parameters { + private: + typedef std::map maptype; //!< Alias typedef for a map of string pairs, used to store the parameter name and value + typedef maptype::const_iterator cmapiterator; //!< constant interator for the map + maptype fParameters; //!< Storage for the parmeter name/value pairs + std::string fCommentString; //!< The comment string is defined in the constructor. + /** \brief This class holds the neccessary input/output operations for each + * line of parameter = value pair. + * + * It is assumed that the parameter does not contain spaces and there is a + * space between the parameter and the equal sign. + */ + class Pair { + /** \brief Output operator for a pair + * + * Outputs a parameter = value string + */ + friend std::ostream & operator << (std::ostream & os, const Pair &pair); + /** \brief Input operator for a pair + * + * Input must be of the form + * \code + * ParameterName = value can be a vector with whitespace + * \endcode + * \note The space between the parameter name and the =-sign is + * mandatory. + */ + friend std::istream & operator >> (std::istream & is, Pair &pair); + private: + std::string fKey, fValue; //!< Store the name and value + public: + /** \brief Create the pair from an input stream. + * + * \param is is the input stream to be parsed. + */ + Pair(std::istream & is) { is >> *this;} + /** \brief Create a pair from two strings. + * + * \param key is the parameter name + * \param value is it value in string form + */ + Pair(const std::string & key, const std::string & value) : fKey(key), fValue(value) {} + /** \brief Create a pair from a std::pair of two strings. + * + * \param pair is the std::pair where the first value is the parameter + * name and the second value is the parameter value. + */ + Pair(const std::pair &pair) : fKey(pair.first), fValue(pair.second) {} + /** \brief Return a std::pair with the Pair info. + * + * \return a std::pair of two string, where the first is the parameter + * name and the second is the parameter value. + */ + std::pair getPair() const {return std::pair(fKey,fValue);} + }; + /** \brief Helper class to print all the parameters. + * + * Used in conjuction with the for_each algorithm to print the pairs. + */ + class Print { + private: + std::ostream & fos; + Print(); + public: + /** \brief The constructor takes an output stream as an argument */ + Print(std::ostream & os) : fos(os) {} + /** \brief Default operator works on std::pairs of strings */ + void operator () (const std::pair & pair){ + fos << Pair(pair); + } + }; + /** \brief Strip comment string and everything after it from the input str. + * + * \param str is the string from which to strip comments. + * \return a copy of the original string after comments have been removed. + * + * This utilises the comment string set in the constructor + */ + std::string StripComments(std::string & str); + + public: + /** \brief Error class for the Parameters class. */ + class ParameterError : public std::exception { + private: + std::string eMessage; //!< store the error message + public: + /** \brief Construct an error with the default message "Error in + * Parmaeters class" + */ + ParameterError() : eMessage("Error in Parameters class"){} + /** \brief Construct an error with a customized error message. + * + * \param message is the error message to display. + */ + ParameterError(const std::string & message) : eMessage(message){} + ~ParameterError() throw(){} + /** \brief Return the error message of this error */ + virtual const char* what () const throw(){ + return(eMessage.c_str()); + } + }; + /** \brief Construct an empty Parameters object. + * + * \param commentString is the comment string; all characters on a line, including and after this string + * in the parsed input streams will be discarded. Defaults to "#". + */ + Parameters(std::string commentString="#") : fCommentString(commentString) {} + /** \brief Construct a Parameters object from an input stream. + * + * \param is is the input stream to be parsed. It will be parsed line by + * line and all comments stripped off. + * \param commentString is the comment string; all characters on a line, including and after this string + * in the parsed input streams will be discarded. Defaults to "#". + */ + Parameters(std::istream & is, std::string commentString="#") : fCommentString(commentString) {ParseIstream(is);} + /** \brief Parse a parameter stream into the Parameters object. + * + * \param is is the input stream to be parsed. It will be parsed line by + * line and all comments stripped off. + */ + void ParseIstream(std::istream & is); + /** \brief Print all the parameters to a output stream + * + * \param os is the output stream to which the parameters will be printed. + */ + void PrintOstream(std::ostream & os) const; + /** \brief Return a parameter value given its name. + * + * This is a templated function and the results depend on the type of the + * input parameter. + * \param parameter is the name of the parameter + * \param value is a reference to a variable which will contain the + * parameter value on return. + * + * The basic input mechanism of c++ is used to parse the value string, so + * only the first non-whitespace containing string will be parsed for the + * value. + */ + template + void getParameter(const std::string & parameter, T & value) const{ + const cmapiterator it = fParameters.find(parameter); + if (it == fParameters.end()) throw(ParameterError("Parameter \""+parameter+"\" not found")); + //Create a stream to read the parameter + std::istringstream iss((*it).second); + iss >> value; + if (iss.fail() || iss.bad()) throw(ParameterError("Error reading parameter \""+parameter+"\"")); + } + /** \brief Return a parameter value as a string given its name. + * + * \param parameter is the name of the parameter + * \param value is a reference to a string which will contain the + * parameter value on return. + * + * The unparsed value string is returned. + */ + void getParameter(const std::string & parameter, std::string & value) const { + const cmapiterator it = fParameters.find(parameter); + if (it == fParameters.end()) throw(ParameterError("Parameter \""+parameter+"\" not found")); + value = (*it).second; + } + /** \brief Return a parameter value as a vector given its name. + * + * This is a templated function and the results depend on the type of the + * input parameter. + * \param parameter is the name of the parameter + * \param vec is a reference to a vector which will contain the + * parameter values on return. + * + * The basic input mechanism of c++ is used to parse the value string, so + * the value string will be parsed to type T for each white space seperated + * object in the value string. + */ + template + void getParameter(const std::string & parameter, std::vector & vec) const{ + const cmapiterator it = fParameters.find(parameter); + if (it == fParameters.end()) throw(ParameterError("Parameter \""+parameter+"\" not found")); + //Clear the vector + vec.clear(); + //Create a stream to read the parameter + std::istringstream iss((*it).second); + while (iss.good()){ + T value; + iss >> value; + if (iss.fail() || iss.bad()) throw(ParameterError("Error reading parameter \""+parameter+"\"")); + vec.push_back(value); + } + } + + /** \brief Set the value of a parameter given its name. + * + * This is a templated function and the results depend on the type of the + * input parameter. + * \param parameter is the name of the parameter + * \param value is a reference to a variable which holds the parameter + * value + * + * The basic output mechanism of c++ is used to convert the value to string. + */ + template + void setParameter(const std::string & parameter, const T & value){ + std::istringstream iss; + iss << value; + fParameters[parameter] = iss.str(); + } + /** \brief Set the value of a parameter given its name. + * + * \param parameter is the name of the parameter + * \param value is a reference to a string which holds the parameter + * value + * + * The value of the parameter is set to that of the string. + */ + void setParameter(const std::string & parameter, const std::string & value){ + fParameters[parameter] = value; + } + friend std::ostream & operator << (std::ostream & os, const Pair &pair); + friend std::istream & operator >> (std::istream & is, Pair &pair); +}; + +#endif diff --git a/source/Particle.h b/source/Particle.h new file mode 100644 index 0000000000000000000000000000000000000000..56fb491d0d8cd18441e47662ff595d2e19bc8401 --- /dev/null +++ b/source/Particle.h @@ -0,0 +1,88 @@ + +//**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****| +// * Particle.h * galprop package * 4/14/2000 +//**"****!****"****!****"****!****"****!****"****!****"****!****"****!****"****| + +#ifndef Particle_h +#define Particle_h + +using namespace std; //AWS20050624 +#include //AWS20050624 +#include //AWS20050624 +#include"Distribution.h" + +class Particle +{ + public: + char name[100]; + int Z,A; + int K_electron; // number of K-electrons (0 or 1) AWS20010731 + double mass; + double t_half; // half-life in years + double primary_abundance; // primary isotopic abundance + + double z_min,z_max,dz; // for 1,2,3D + double r_min,r_max,dr; // for 2D + double x_min,x_max,dx,y_min,y_max,dy; // for 3D + + double p_min, p_max, p_factor; // momentum start, end, factor + double Ekin_min, Ekin_max, Ekin_factor; // kinetic energy/nucleon start,end,factor + + int n_spatial_dimensions;// 1,2,3D + int n_pgrid; // number of points in momentum + int n_rgrid; // number of points in radius (2D) + int n_zgrid; // number of points in z (1D,2D) + int n_xgrid; // number of points in x (3D) + int n_ygrid; // number of points in y (3D) + + double *x; // x grid + double *y; // y grid + double *z; // z grid + double *r; // r grid + + double *p; // total momentum of particle in MV + double *Etot; // total energy of particle in MeV + double *Ekin; // kinetic energy per nucleon in MeV + double *beta; // velocity/c + double *gamma; // Etot/mass + double *rigidity; // rigidity p/Z in MV + + char species[20]; + + int arrays_assigned; + + double normalization_factor;// to normalize protons or electrons AWS20010121 + + Distribution cr_density; + Distribution primary_source_function; + Distribution secondary_source_function; + + Distribution fragment;// fragmentation destruction rate + Distribution decay ;// radioactive decay rate + Distribution dpdt; // momentum change rate + Distribution Dxx ; // spatial diffusion coefficient + Distribution Dpp ; // momentum diffusion coefficient + Distribution v_conv; // convection velocity in z-direction + +//interface functions prototypes + void init(char *name_,int Z_,int A_,int K_electron_,double t_half_, //AWS20010731 + double r_min_, double r_max_, double dr_, + double z_min_, double z_max_, double dz_, + double p_min_,double p_max_,double p_factor_, + double Ekin_min_,double Ekin_max_,double Ekin_factor_, + char *p_Ekin_grid); + void init(char *name_,int Z_,int A_,int K_electron_,double t_half_, //AWS20010731 + double x_min_, double x_max_, double dx_, + double y_min_, double y_max_, double dy_, + double z_min_, double z_max_, double dz_, + double p_min_,double p_max_,double p_factor_, + double Ekin_min_,double Ekin_max_,double Ekin_factor_, + char *p_Ekin_grid); + void init(); // flags object as having arrays not yet assigned + int delete_transport_arrays(); + int create_transport_arrays(); + void print(); + Particle operator=(Particle particle); +}; + +#endif diff --git a/source/PhysicalConstants.h b/source/PhysicalConstants.h new file mode 100644 index 0000000000000000000000000000000000000000..dca28b4c1a7842b5c65cb78a43760ff5f7e0b8c8 --- /dev/null +++ b/source/PhysicalConstants.h @@ -0,0 +1,199 @@ +#ifndef _utl_PhysicalConstants_h_ +#define _utl_PhysicalConstants_h_ + +#include + +#include + +namespace utl { + + const double kPi = M_PI; + const double kExp = exp(1.0); + const double kSqrtTwo = sqrt(2.0); + const double kSqrtThree = sqrt(3.0); + const double kOneOnThree = 1.0/3.0; + const double kFourPiOnThree = 4.0*kPi/3.0; + const double kPiOnTwo = kPi/2.0; + const double kTwoPi = 2.0*kPi; + const double kFourPi = 2.0*kTwoPi; + const double kOneOnPi = 1.0/kPi; + const double kOneOnTwoPi = 1.0/kTwoPi; + const double kOneOnFourPi = 1.0/kFourPi; + const double kConvertDegreesToRadians = kPi/180.0; + + // All taken from PDG data tables (2002) + + // Physical constants + + const double kSpeedOfLight_SI = 299792458.0; + const double kSpeedOfLight = kSpeedOfLight_SI*m/s; + const double kPlanck_SI = 6.62606876e-34; + const double kPlanckReduced_SI = kPlanck_SI*kOneOnTwoPi; + const double kPlanck = kPlanck_SI*joule*s; + const double kPlanckReduced = kPlanckReduced_SI*joule*s; + const double kMuZero_SI = 4.0*kPi*1.0e-7; + const double kMuZero = kMuZero_SI*newton/(ampere*ampere); + const double kBoltzmann_SI = 1.3806503e-23; + const double kBoltzmann = kBoltzmann_SI*joule/kelvin; + + const double kPlanckTimesSpeedOfLight_SI = kPlanck_SI*kSpeedOfLight_SI; + const double kPlanckTimesSpeedOfLightSquared_SI = + kPlanckTimesSpeedOfLight_SI*kPlanckTimesSpeedOfLight_SI; + + // Particle and other masses + + const double kMassConversion_SI = e_SI/(kSpeedOfLight_SI*kSpeedOfLight_SI); + + const double kHydrogenMass_SI = 1.6735e-27; + const double kHydrogenMass = kHydrogenMass_SI*kg; + const double kSolarMass_SI = 1.98892e30; + const double kSolarMass = kSolarMass_SI*kg; + + const double kElectronMass = 0.510998902*MeV; + const double kElectronMass_SI = kElectronMass*kMassConversion_SI; + const double kMuonMass = 105.658357*MeV; + const double kMuonMass_SI = kMuonMass*kMassConversion_SI; + const double kTauMass = 1776.99*MeV; + const double kTauMass_SI = kTauMass*kMassConversion_SI; + + const double kProtonMass = 938.271998*MeV; + const double kProtonMass_SI = kProtonMass*kMassConversion_SI; + const double kNeutronMass = 939.56533*MeV; + const double kNeutronMass_SI = kNeutronMass*kMassConversion_SI; + const double kDeuteronMass = 1875.612762*MeV; + const double kDeuteronMass_SI = kDeuteronMass*kMassConversion_SI; + + const double kLambdaMass = 1115.683*MeV; + const double kLambdaMass_SI = kLambdaMass*kMassConversion_SI; + const double kSigmaZeroMass = 1192.642*MeV; + const double kSigmaZeroMass_SI = kSigmaZeroMass*kMassConversion_SI; + const double kSigmaPlusMass = 1189.37*MeV; + const double kSigmaPlusMass_SI = kSigmaPlusMass*kMassConversion_SI; + const double kSigmaMinusMass = 1197.449*MeV; + const double kSigmaMinusMass_SI = kSigmaMinusMass*kMassConversion_SI; + const double kXiZeroMass = 1314.83*MeV; + const double kXiZeroMass_SI = kXiZeroMass*kMassConversion_SI; + const double kXiMinusMass = 1321.31*MeV; + const double kXiMinusMass_SI = kXiMinusMass*kMassConversion_SI; + const double kOmegaMinusMass = 1672.45*MeV; + const double kOmegaMinusMass_SI = kOmegaMinusMass*kMassConversion_SI; + + const double kPiZeroMass = 134.9766*MeV; + const double kPiZeroMass_SI = kPiZeroMass*kMassConversion_SI; + const double kPiChargedMass = 139.57018*MeV; + const double kPiChargedMass_SI = kPiChargedMass*kMassConversion_SI; + const double kKaonChargedMass = 493.677*MeV; + const double kKaonChargedMass_SI = kKaonChargedMass*kMassConversion_SI; + + const double kAtomicMassUnit_SI = 1.660538e-27; + + const double kCarbonMass_SI = 12.0107*kAtomicMassUnit_SI; + const double kOxygenMass_SI = 15.9994*kAtomicMassUnit_SI; + const double kMagnesiumMass_SI = 24.3050*kAtomicMassUnit_SI; + const double kSiliconMass_SI = 28.0855*kAtomicMassUnit_SI; + const double kIronMass_SI = 55.845*kAtomicMassUnit_SI; + + const double kSilicateMass_SI = kMagnesiumMass_SI + kIronMass_SI + + kSiliconMass_SI + 4.0*kOxygenMass_SI; // Silicate is MgFeSi0_4 + + const double kGraphiteDensity_SI = 2.24*(gram/kilogram)/pow(cm/m, 3.0); + const double kSilicateDensity_SI = 3.50*(gram/kilogram)/pow(cm/m, 3.0); + + // Particle lifetimes + + const double kMuonLifetime = 2.19703e-6*s; + + const double kNeutronLifetime = 885.7*s; + + const double kLambdaLifetime = 2.632e-10*s; + const double kSigmaZeroLifetime = 7.4e-20*s; + const double kSigmaPlusLifetime = 0.8018e-10*s; + const double kSigmaMinusLifetime = 1.479e-10*s; + const double kXiZeroLifetime = 2.9e-10*s; + const double kXiMinusLifetime = 1.639e-10*s; + const double kOmegaMinusLifetime = 0.821-10*s; + + const double kPiZeroLifetime = 8.4e-17*s; + const double kPiChargedLifetime = 2.6033e-8*s; + const double kKaonChargedLifetime = 1.2384e-8*s; + + // Derived constants + + const double kEpsilonZero_SI = + 1.0/(kMuZero_SI*kSpeedOfLight_SI*kSpeedOfLight_SI); + const double kAlpha = (e_SI*e_SI)/ + (4.0*kPi*kEpsilonZero_SI*kPlanckReduced_SI*kSpeedOfLight_SI); + const double kElectronRadius_SI = (e_SI*e_SI)/ + (4.0*kPi*kEpsilonZero_SI*kElectronMass_SI* + kSpeedOfLight_SI*kSpeedOfLight_SI); + const double kThomsonCrossSection_SI = + 8.0*kPi*kElectronRadius_SI*kElectronRadius_SI/3.0; + + // Distance conversions + + const double kParsec = 3.0856775807e+16*m; + const double kKiloParsec = kParsec*1.0e+3; + const double kMegaParsec = kKiloParsec*1.0e+3; + + const double pc = kParsec; + const double kpc = kKiloParsec; + const double Mpc = kMegaParsec; + + // Some other conversions and constants + + const double kGalactocentricRadiusSun = 8.5*kpc; + + const double kYearToSec = 365.25*24.0*60.0*60.0; + const double kSolarLuminosity_SI = 3.846e26; + const double kSolarLuminosity = kSolarLuminosity_SI*watt; + + // Conversion from LSun/kpc^3 to eV cm^-3 sr^-1 kpc^-1 + // -> LSun (J s^-1) / (kpc/cm)^3 * kpc/cm / (4*Pi*c) + const double kWainscoatConversion = + (kSolarLuminosity_SI/e_SI)/pow(kpc/cm, 3.0)* + (kpc/kSpeedOfLight_SI)*(1.0/(4.0*kPi)); + + // Conversion from LSun/pc^3 to eV cm^-3 sr^-1 kpc^-1 + // -> LSun (J s^-1) / (pc/cm)^3 * kpc/cm / (4*Pi*c) + const double kMathisConversion = + (kSolarLuminosity_SI/e_SI)/pow(pc/cm, 3.0)* + (kpc/kSpeedOfLight_SI)*(1.0/(4.0*kPi)); + + // Conversion from MJy sr^-1 kpc^-1 -> eV cm^-3 sr^-1 kpc^-1. + // Take M == 10^6, Jy == 10^-26 W m^-2 Hz^-1 -> 10^-26/e_SI eV m^-2 Hz^-1 + // 1/c converts cm^-2 s^-1 to cm^-3 if in kpc s^-1 + const double kFreudenreichConversion = + 1.0e6*(1.0e-26/e_SI)*(cm/kSpeedOfLight_SI)*(cm*cm); + + // Conversion of W sr^-1 H-atom^-1 H-atom cm^-3 -> eV cm^-3 sr^-1 kpc^-1 + // W -> J s^-1/e_SI -> eV s^-1 + // 1/c converts cm^-2 s^-1 -> cm^-2 kpc^-1 if c in kpc s^-1 + // So get W sr^-1 cm^-3 -> eV cm^-3 sr^-1 kpc^-1 + // times factor 10^-32 associated with emissivity from table + + const double kSodroskiConversion = (1.0e-32/e_SI)*(kpc/kSpeedOfLight_SI); + + // Temperature, absolute magnitude at V, bolometric correction and magnitude + // for Sun + + const double kTemperatureSun = 5770; + const double kMVSun = 4.82; + const double kBCSun = -0.05; + const double kMBolometricSun = kMVSun + kBCSun; + + const double kFrequencyConstant = + 2.0*kPlanck_SI/pow(kSpeedOfLight_SI, 2.0); + const double kFrequencyExponentConstant = kPlanck_SI/kBoltzmann_SI; + const double kWavelengthConstant = + 2.0*kPlanck_SI*pow(kSpeedOfLight_SI, 2.0); + const double kWavelengthExponentConstant = + kPlanck_SI*kSpeedOfLight_SI/kBoltzmann_SI; + const double kPlanckIntegral = 2.0*pow(kBoltzmann_SI*kPi, 4.0) + /(15.0*pow(kPlanck_SI*kSpeedOfLight_SI, 2.0)*kPlanck_SI); + const double kStefanBoltzmann = kPlanckIntegral*kPi; + +} + +#endif + + diff --git a/source/Psf.h b/source/Psf.h new file mode 100644 index 0000000000000000000000000000000000000000..355b1b937a5e74f5a3c889b5e5e5cfaae5dbbfcd --- /dev/null +++ b/source/Psf.h @@ -0,0 +1,158 @@ +#ifndef PSF_H +#define PSF_H + +#include +#include +#include +#include "healpix_base.h" +#include "healpix_map.h" +#include "Exposure.h" + +/** \brief Handles the conversion of the PSF from a table to a healpix map. + * + * Reads the point spread function from a fits table, having the same format as + * the output of gtpsf. + */ +class Psf{ + public: + /** \brief Standard constructor */ + Psf() {} + /** \brief Construct the object from a fits file. + * + * \param fileName gives the name of the file to open. This is a fits + * file, where the point spread function is defined in an extension named + * PSF. This extension has a table with 3 columns, energy, exposure and + * psf, where the last column is a vector containing the function values + * for each value of theta. The values of theta are stored in another + * extension, THETA in a table with a single column containing the theta + * values in degrees. The energy in the PSF table is supposed to be in MeV. + */ + Psf(const std::string &fileName); + /** \brief Construct the object from a 2 dimensional vector. + * + * \param psf is a 2 dimensional vector, where the first index corresponds + * to theta and the second index to energy + * \param energies a vector of energy values in MeV + * \param theta a vector of angles in degrees + */ + Psf(const std::vector< std::vector > & psf, const std::vector & energies, const std::vector &theta); + /** \brief Return a Healpix_Map of the point spread function for a given + * energy. + * + * \param energy is the energy for which to evaluate the psf. The closest + * available energy is actually used, no interpolation performed. + * \param order gives the order of the map. If it is 0 (the default) the + * order is calculated from the resolution of the psf. + * \param center gives the direction of the psf. This is useful when + * calculating the maps of point sources. The method used to calculate the + * map when the pointing is (0,0) (the default) is more accurate than if it + * is rotated. + * \return A Healpix_Map with the point spread function in the correct + * place. + */ + Healpix_Map getPsfMap(double energy, int order = 0, const pointing ¢er = pointing(0,0), double psfFraction=1) const; + /** \brief Returns an effective PSF map given an energy range and a power + * law index for the energy dependance of the source flux. + * + * \param eMin the minimum energy + * \param eMax the maximum energy + * \param index the power law index of the source flux + * \param order the order of the resulting map. If it is 0 (the default) + * then the order is approximated from the psf resolution. + * \param center the position of the psf center. Calculation is more + * accurate when it is (0,0) (the default). + * \return The power law weighted map of the psf in the correct position. + */ + Healpix_Map getEffectiveMap(double eMin, double eMax, double index, int order = 0, const pointing ¢er = pointing(0,0), double psfFraction=1) const; + /** \brief Convolve a skymap with the psf + * + * \param skyMap is the input skymap to be convolved. It should be in RING + * scheme to make the computation faster. + * + * \param deconvolve can be set to true for deconvolution. + * Does not work very well. + * + * \return the convolved skymap, which has the same structure as the + * original one. + * + * Uses the fast spherical harmonics decomposition of healpix to do the + * convolution. Assumes a homogeneous spherically symmetric psf, since that + * is the only thing that can be handled by this method. In the case of + * GLAST, this should not matter much. + */ + Skymap convolve(const Skymap & skyMap, bool deconvolve=false) const; + /** \brief Convolve a skymap with the psf + * + * \param skyMap is the input skymap to be convolved. It should be in RING + * scheme to make the computation faster. + * + * \param fraction is set to the fraction of psf to take into + * account. Will not speed up the convolution but combined + * with SparseSkymap it can make for less memory consumption. + * + * \return the convolved skymap, which has the same structure as the + * original one. + * + * Uses brute force to do the convolution, still assuming a + * homogeneous spherically symmetric psf. + */ + Skymap convolveBruteForce(const Skymap & skyMap, double fraction = 1.0) const; + /** \brief Given an energy range and a power law index, calculate an effective psf vector, weighted with the power law. + * + * \param eMin the minimum energy + * \param eMax the maximum energy + * \param index the power law index of the source flux + * \return a vector of weighted psf values, corresponding to the theta + * values of the actual underlying psf. + */ + std::vector getEffectivePsf(double eMin, double eMax, double index) const; + /** \brief Given an energy value, return the approximate width of the psf + * at the energy closest to the expected value. + * + * \param energy is approximated to the nearest energy point of the psf + * object. + * \param fraction is the psf containment limit + * \return the value of theta in degrees + */ + double getPsfWidth(double energy, double fraction) const; + /** \brief Return a constant reference to the raw psf values */ + const std::vector< std::vector > &getPsfVector() const; + /** \brief Return a constant reference to the theta values */ + const std::vector &getTheta() const; + /** \brief Return a constant reference to the energy values */ + const std::vector &getEnergies() const; + private: + /** \brief Read the psf table from fits file. + * + * \param fileName is the name of fits file to read. + */ + void ReadPsfFits(const std::string &fileName); + /** \brief Calculate Lagrange polynomyal interpolation of the energy + * dependance of the Psf. + * + * Used when calculating the effective psf. + */ + void CalculateInterpolation(); + /** \brief Given a psf vector values, return a map. + * + * \param center the position of the psf + * \param psf a vector containing the psf values as a function of theta + * \param order is the order of the output map + * \return The psf in a Healpix_Map. + */ + Healpix_Map createMap(const pointing & center, const std::vector & psf, int order, double psfFraction=1) const; + /** \brief Given a psf, return the approximate width of the psf + * + * \param psf is a vector containing the psf values as a + * function of theta + * \param energy is approximated to the nearest energy point of the psf + * object. + * \param fraction is the psf containment limit + * \return the value of theta in degrees + */ + double getPsfWidth(const std::vector & psf, double fraction) const; + std::vector fEnergies, fTheta; // > fPsf; // > > fLagrange; // +#include +#include +#include + +#include + +#include + +namespace rf { + + // Explicit assumption made is that the supplied radiation field + // is generated on a regular spatial grid. That is, there are n*m entries + // in the grid. For each entry in the grid, there is a range of validity + // (given by the cell size). Outside this range the radiation field is + // determined by interpolation within the grid. For points outside the grid + // the radiation field is considered zero + + class RadiationField { + + public: + + enum STELLARCOMPONENT { TOTAL, DIRECT, SCATTERED, TRANSIENT, THERMAL }; + + typedef CLHEP::Hep3Vector ThreeVector; + + RadiationField(); + + RadiationField(const std::string& filename, + const std::valarray& freq, + int rebinnedSkymapOrder = 0); + + ~RadiationField(); + + // Note: The method below is *very* expensive if the desired + // healpix order is different from the rebinned order supplied + // in the constructor above. Use with caution! + + const Skymap GetSkymap(const ThreeVector& pos, + const STELLARCOMPONENT component, + const int healpixOrder); + + private: + + const Skymap GetSkymap(const ThreeVector& pos, + const STELLARCOMPONENT component); + + std::valarray fCacheBuilt; + + string fPrefix; + + int fRebinnedSkymapOrder; + + unsigned int fNumberOfComponents; + + double fLuminosity, fDustMass; + + std::valarray fFrequency, fWavelength, fStellarComponentLuminosity; + + std::vector fStellarComponentName; + + std::map > fFilenameData; + + std::vector< std::vector > fPositionData, fRangeData; + + std::valarray fXData, fYData, fZData, fXRangeData, fYRangeData, fZRangeData; + + std::vector< std::vector< std::vector* > > fFilenameOrderedData; + + std::vector< std::vector< std::vector* > > > fSkymapOrderedData; + + void ReadRadiationField(const std::string& filename); + + void BuildSkymapCache(const unsigned long component); + void FlushSkymapCache(const unsigned long component); + + void FlushSkymapCache(); + + void ClearData(); + + }; + +} + +#endif diff --git a/source/Region.h b/source/Region.h new file mode 100644 index 0000000000000000000000000000000000000000..0bef690e67b1d2146e2af456e807402052b51e86 --- /dev/null +++ b/source/Region.h @@ -0,0 +1,59 @@ +#ifndef REGION_H +#define REGION_H + +#include "Coordinate.h" +#include + +enum RegionType { + DISC, //!< Region selected by query_disk method of Healpix + RECTANGLE //!< A square region on a CAR projection +}; + +class SkyRegion { + public: + SkyRegion(); + /** \brief Create a region with the given parameters + * + * \parameter type is either DISC or RECTANGLE + * \parameter center is the center of the region + * \parameter radius is the radius of the DISC, ignored if the region is + * a RECTANGLE + * \parameter deltal is the longitude range of the RECTANGLE, ignored if the + * region is a DISC + * \parameter deltab is the latitude range of the RECTANGLE, ignored if the + * region is a DISC + * + * The RECTANGLE region takes into account that we are selecting on a + * sphere and does not loop from the north pole to the south pole. A + * rectangle that traverses either pole results in two triangles, both with + * a vertex at the pole (unless the latitude range covers the whole + * sphere). + */ + SkyRegion(RegionType type, const SM::Coordinate ¢er, double radius, double deltal, double deltab); + //! Methods to get information on the region + const RegionType & type() const; + const SM::Coordinate & center() const; + const double & radius() const; + const double & deltal() const; + const double & deltab() const; + /** \brief Read a region from a stream + * + * The format of the region is + * + * where is either RECTANGLE or DISC. DISC requires 3 arguments + * (center and radius) while RECTANGLE requires 4 (center and width and + * height). All dimensions are in degrees. + */ + friend std::istream & operator >> (std::istream &is, SkyRegion ®ion); + /** \brief Write a region to a stream + * + * The format is the same as for reading + */ + friend std::ostream & operator << (std::ostream &os, const SkyRegion ®ion); + private: + RegionType ftype; + //Parameters of the region + SM::Coordinate fcenter; + double fradius, fdeltal, fdeltab; +}; +#endif diff --git a/source/SPIERSP.h b/source/SPIERSP.h new file mode 100644 index 0000000000000000000000000000000000000000..bbdc44d87cc551a6009b18c665acd01b1742f456 --- /dev/null +++ b/source/SPIERSP.h @@ -0,0 +1,28 @@ +// SPI energy response + +#include"RMF.h" +#include"IRF.h" + +class SPIERSP +{ + public: + + RMF rmf; + IRF irf; + + + + + double *E_in, *E_out_lo,*E_out_hi; // interpolated response energies + int n_E_in,n_E_out; // interpolated response energies + double *response; // interpolated response grid + + int read (char *filename,int idet1,int idet2, + double *E_in_, int n_E_in_, + double *E_out_lo_,double *E_out_hi_,int n_E_out_); + + void print(); + + + double R(double E_in, double E_out); +}; diff --git a/source/Singleton.h b/source/Singleton.h new file mode 100644 index 0000000000000000000000000000000000000000..9b5bb015f945db180f3a97f8b4bd075628a4a376 --- /dev/null +++ b/source/Singleton.h @@ -0,0 +1,125 @@ +#ifndef _utl_Singleton_h_ +#define _utl_Singleton_h_ + +namespace utl { + + /** + * \class Singleton Singleton.h utl/Singleton.h + * + * \brief Curiously Recurring Template Pattern (CRTP) for Meyers singleton + * + * The singleton class is implemented as follows + * \code + * #include + * + * class SomeClass : public utl::Singleton { + * ... + * private: + * // prevent creation, destruction + * SomeClass() { } + * ~SomeClass() { } + * + * friend class utl::Singleton; + * }; + * \endcode + * Singleton automatically prevents copying of the derived class. + * + * \author Darko Veberic + * \date 9 Aug 2006 + * \version $Id: Singleton.h 7249 2008-05-07 22:03:06Z darko $ + * \ingroup stl + */ + + template + class Singleton { + public: + static + T& + GetInstance() +#ifdef __MAKECINT__ + ; +#else + { + static T instance; + return instance; + } +#endif + + protected: + // derived class can call ctor and dtor + Singleton() { } + ~Singleton() { } + + private: + // no one should do copies + Singleton(const Singleton&); + Singleton& operator=(const Singleton&); + + }; + + + /** + * \class LeakingSingleton Singleton.h utl/Singleton.h + * + * \brief CRTP for leaking singleton + * + * This type of creation (Gamma singleton) leaks the object at + * the end of the run, i.e. class destructor does not get called + * in at_exit(). + * + * This singleton can be implemented as follows + * \code + * #include + * + * class SomeClass : public utl::LeakingSingleton { + * ... + * private: + * // prevent creation, destruction + * SomeClass() { } + * ~SomeClass() { } + * + * friend class utl::LeakingSingleton; + * }; + * \endcode + * LeakingSingleton automatically prevents copying of the derived + * class. + * + * \author Darko Veberic + * \date 9 Aug 2006 + * \version $Id: Singleton.h 7249 2008-05-07 22:03:06Z darko $ + * \ingroup stl + */ + + template + class LeakingSingleton { + public: + static + T& + GetInstance() + { + if (!fgInstance) + fgInstance = new T; + return *fgInstance; + } + + protected: + // derived class can call ctor and dtor + LeakingSingleton() { } + // will never get called anyway + ~LeakingSingleton() { } + + private: + // no one should do copies + LeakingSingleton(const LeakingSingleton&); + LeakingSingleton& operator=(const LeakingSingleton&); + + static T* fgInstance; + }; + + template + T* LeakingSingleton::fgInstance = 0; + +} + + +#endif diff --git a/source/Skymap.h b/source/Skymap.h new file mode 100644 index 0000000000000000000000000000000000000000..cec6f9637a60399933665feccc2954141841e0d6 --- /dev/null +++ b/source/Skymap.h @@ -0,0 +1,1772 @@ +/**\class Skymap + * \brief Skymaps for gamma ray calculation in galprop + * + * The skymap class is used to store skymaps calculated in gardian + * The information is stored in a healpix format. + * Internally we use a std::valarray to store the spectra for each pixel + * \author Gudlaugur Johannesson + */ + +#ifndef SKYMAP_H +#define SKYMAP_H + +#include "HealpixBaseExtended.h" +#include "Coordinate.h" +#include "PhysicalConstants.h" +#include "ArraySlice.h" + +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +//Coordinate transformation if we find libastro +#ifdef HAVE_ASTRO +extern "C"{ +#include +} +//Do to stupid convention in astro, we must create wrapper functions +void eq_gal2(double mj, double DEC, double RA, double *b, double *l); +void eq_ecl2(double mj, double DEC, double RA, double *lt, double *lg); +void gal_eq2(double mj, double b, double l, double *DEC, double *RA); +void ecl_eq2(double mj, double lt, double lg, double *DEC, double *RA); +//Functions to convert from gal to ecl and wise versa +void gal_ecl2(double mj, double b, double l, double *lt, double *lg); +void ecl_gal2(double mj, double lt, double lg, double *b, double *l); +void empty(double mj, double x, double y, double *z, double *w); +#endif + + + +/**\brief Store a full sky map in a healpix pixelation. + * + * This class is templated to allow for a flexible storage. Inherits from Healpix_Base so all its methods + * are available. Can only be initialized with a whole number order, i.e. nside = 2^order. + */ +template +class Skymap : public HealpixBaseExtended { + public: + enum CoordSys {EQ, GAL, ECL}; + private: + std::valarray fMap; //!< A valarray to store the values, we use slices to get single spectra + std::valarray fSpectra; //!< The values at which the spectra is evaluated, usually energy for photon maps. + std::valarray fSpecMin, fSpecMax; //!< If the skymap is binned, we store maximum and minimum as well as mid value for compatibility. + bool fBinned; //!< Tells if the skymap is spectrally binned or instantaneous. + + /** \brief Returns true if the maps have the same size and same fSpectra */ + bool equiv(const Skymap & otherMap) const{ + //First we check for binning equivalance + bool eq = (fBinned == otherMap.fBinned); + bool spectra = (fSpectra.size() == otherMap.fSpectra.size()); + if (eq && spectra) { + //Check for equivalence of spectra + if (fBinned) { + for (int i = 0; i < fSpectra.size(); ++i){ + spectra &= (fabs((fSpecMin[i] - otherMap.fSpecMin[i])/fSpecMin[i]) < 1e-6); + spectra &= (fabs((fSpecMax[i] - otherMap.fSpecMax[i])/fSpecMax[i]) < 1e-6); + } + }else{ + for (int i = 0; i < fSpectra.size(); ++i){ + spectra &= (fabs((fSpectra[i] - otherMap.fSpectra[i])/fSpectra[i]) < 1e-6); + } + } + } + return (scheme_ == otherMap.scheme_ && order_ == otherMap.order_ && spectra && eq); + } + + /**\brief Convert a CAR map to healpix format + * + * Given arrays for CRVAL, CDELT and CRPIX (at the center of the pixel), + * turn the mapcube into pixels. The skymap has to be resized properly + * before entering this method. The \a image is assumed to be three + * dimensional in the standard FITS way, longitude evolving fastest, then + * latitude and last spectra. + * + * \param image is the valarray of image data, as from ccfits + * \param crval is the CRVAL array from the fits header + * \param crpix is the CRPIX array from the fits header + * \param cdelt is the CDELT array from the fits header + */ + void fillMapcube(const std::valarray &image, const std::valarray &axis, const std::valarray &crval, const std::valarray &cdelt, const std::valarray &crpix, bool correctSA = false){ + //Divide each HEALPix pixel into 256 smaller ones and fill those with the + //value from the map pixel directly beneeth its center. The HEALPix + //pixel value is the average of those 256 smaller ones. + const int dOrder = std::min(5,13-order_); //Order can not exceed 13 + const int nPix = (1< totSpectra(T(0), fSpectra.size()); + //Loop over all of the subpixels + pointing pnt1 = pix2ang(p); + for (int sp = pp*nPix; sp<(pp+1)*nPix; ++sp){ + const pointing pnt = finer.pix2ang(sp); + + //Find the correct index in the array, first the longitude which can + //loop + double dil = (pnt.phi*180/utl::kPi-crval[0])/cdelt[0] + crpix[0] + 0.5; + //If the pixel value is negative, we must loop it over the 360 + //degrees + if (dil < 1) { + dil += fabs(360./cdelt[0]); + //Pixel values above the boundaries could also indicate looping + } else if (dil > axis[0]) { + dil -= fabs(360./cdelt[0]); + } + int il = int(dil - 1); //Since the pixels are 1 based and the array is 0 based + + //Then the latitude, which can't loop + int ib = int( ( (90 - pnt.theta*180/utl::kPi)-crval[1])/cdelt[1] + crpix[1] + 0.5) - 1; + + + //They must both be within bounds to do something useful + if (ib >= 0 && ib < axis[1] && il >= 0 && il < axis[0]) { + const int ind1 = il + ib*axis[0]; + //We must divide with the solid angle if requested, since the + //method does not work unless the value of the pixel is independent + //of the solid angle of the pixel + double lbSolidAngle = 1; + if (correctSA) { + double bMiddle = crval[1] + (ib+1 - crpix[1])*cdelt[1]; + //Note that sin(x-pi/2) = -cos(x) + lbSolidAngle = cdelt[0]*utl::kPi/180.*(cos((bMiddle-cdelt[1]/2.)*utl::kPi/180.) - cos((bMiddle+cdelt[1]/2.)*utl::kPi/180.)); + } + for (int is=0; is spectra, const Healpix_Ordering_Scheme scheme=RING, const T & default_value = T(0)){ + Resize(order,spectra,scheme,default_value); + } + + /**\brief Constructor that takes an order of the skymap, the + * boundaries of the spectral bins, an ordering scheme and a default value for the map. + * + * \param order is the order of the healpix coordinate system. The number + * of sides per base pixel, nside, is 2^order. + * \param specMin are the lower boundaries of the spectral bins + * \param specMax are the upper boundaries of the spectral bins + * \param scheme is the ordering scheme of the healpix array + * \param default_value is the default value for the skymap + */ + Skymap(const int order, const std::valarray specMin, const std::valarray specMax, const Healpix_Ordering_Scheme scheme=RING, const T & default_value = T(0)){ + Resize(order,specMin,specMax,scheme,default_value); + } + + /**\brief Constructor that takes an order of the skymap, the size of the spectra, + * an ordering scheme and a default value for the map. + * + * \param order is the order of the healpix coordinate system. The number + * of sides per base pixel, nside, is 2^order. + * \param nSpectra is the size of the spectra + * \param scheme is the ordering scheme of the healpix array (defaults to RING) + * \param default_value is the default value for the skymap (defaults to T(0)) + * + * The values at which the spectra is located is set to 1 in all cases. This is done so we can take the + * logarithm of the values, but it is expected that the values at which the spectra is evaluated is in most + * cases logarithmically distributed. + */ + Skymap(const int order, const int nSpectra, const Healpix_Ordering_Scheme scheme=RING, const T & default_value = T(0)){ + Resize(order,nSpectra,scheme,default_value); + } + + /**\brief Construct a skymap from file, either a healpix fits file or CAR + * projected fits image. + * + * \param fileName is the name of the file to be opened + */ + Skymap(const std::string & fileName, int order=-1) { + load(fileName, order); + } + + /**\brief Copy constructor */ + Skymap(const Skymap & oldMap){ + if (oldMap.fBinned){ + Resize(oldMap.order_, oldMap.fSpecMin, oldMap.fSpecMax, oldMap.scheme_); + }else{ + Resize(oldMap.order_, oldMap.fSpectra, oldMap.scheme_); + } + fMap = oldMap.fMap; + } + + /**\brief Resize the map to order and size of spectra + * + * \param order is the order of the healpix coordinate system + * \param nSpectra is the size of the spectral array + * \param scheme is the ordering scheme of the healpix map (defaults to RING) + * \param defaultValue is the default value for the map (defaults to T(0)) + * + * This method destroys the data in the skymap. The value at which the spectra is evaluated is set to 1. + */ + void Resize(const int order, const int nSpectra, const Healpix_Ordering_Scheme scheme=RING, const T & defaultValue = T(0)){ + std::valarray spectra(1.0,nSpectra); + Resize(order, spectra, scheme, defaultValue); + } + + /**\brief Resize the map to order and size of spectra + * + * \param order is the order of the healpix coordinate system + * \param spectra is the values at which the spectra is evaluated + * \param scheme is the ordering scheme of the healpix map (defaults to RING) + * \param defaultValue is the default value for the map (defaults to T(0)) + * + * This method destroys the data in the skymap. + */ + void Resize(const int order, const std::valarray &spectra, const Healpix_Ordering_Scheme scheme=RING, const T & defaultValue = T(0)){ + Set(order, scheme); // Call the construcor for healpix with correct arguments + fSpectra.resize(spectra.size()); + fSpectra = spectra; + fBinned = false; + fMap.resize(Npix()*fSpectra.size(), defaultValue); + } + + /**\brief Resize the map to order and size of spectra in a binned fashion + * + * \param order is the order of the healpix coordinate system + * \param specMin are the lower boundaries of the spectral bins + * \param specMax are the upper boundaries of the spectral bins + * \param scheme is the ordering scheme of the healpix map (defaults to RING) + * \param defaultValue is the default value for the map (defaults to T(0)) + * + * This method destroys the data in the skymap. + */ + void Resize(const int order, const std::valarray &specMin, const std::valarray &specMax, const Healpix_Ordering_Scheme scheme=RING, const T & defaultValue = T(0)){ + //The spectral sizes must be the same + if (specMin.size() != specMax.size()){ + std::cerr<<"Spectral sizes not equal for boundary arrays"< values(1,""); + keywords[0] = "PIXTYPE"; + values[0] = "HEALPIX"; + CCfits::FITS fits(fileName, CCfits::Read, keywords, values); + + //A reference to the table containing the skymap data + CCfits::ExtHDU &skymapTable = fits.currentExtension(); + + //Read the keywords to set up the skymap object + int nRows, nSpectra, nSide,TFIELDS; + std::string ordering; + skymapTable.readKey("NAXIS2", nRows); + + skymapTable.readKey("TFIELDS", TFIELDS); // number of columns in table + nSpectra=TFIELDS; // equal to number of energies if column format + + if(TFIELDS==1) // vector format has only one column, and uses NBRBINS for number of energies + skymapTable.readKey("NBRBINS", nSpectra); + + skymapTable.readKey("NSIDE", nSide ); + skymapTable.readKey("ORDERING", ordering); + //Calculate the order + int hporder = int(log(double(nSide))/log(2.0)+0.1);// original: int order = int(log(nSide)/log(2)+0.1); AWS20080519 + + //Try to find the EBOUNDS or ENERGIES extensions, either must exist + bool foundEnExt = false; + //First try the energies table + try { + fits.read(std::vector (1,"ENERGIES")); + CCfits::ExtHDU &energyTable = fits.extension("ENERGIES"); + std::valarray energies; + energyTable.column(1).read(energies, 1, nSpectra); + + if (ordering == "NEST") { + if (Order() != hporder || nSpectra != fSpectra.size() || Scheme() != NEST ) { + Resize(hporder, energies, NEST); + } else { + setSpectra(energies); + } + } else { + if (Order() != hporder || nSpectra != fSpectra.size() || Scheme() != RING ) { + Resize(hporder, energies, RING); + } else { + setSpectra(energies); + } + } + foundEnExt=true; + } catch (CCfits::FITS::NoSuchHDU) { + try { //Then EBOUNDS table + fits.read(std::vector (1,"EBOUNDS")); + CCfits::ExtHDU &energyTable = fits.extension("EBOUNDS"); + std::valarray eMin, eMax; + energyTable.column(2).read(eMin, 1, nSpectra); + energyTable.column(3).read(eMax, 1, nSpectra); + + // Fermi-LAT Science Tools counts format has EBOUNDS in keV: convert to MeV galprop/galplot standard + string TUNIT2; + energyTable.readKey("TUNIT2",TUNIT2); + if(TUNIT2=="keV") {eMin*=1.0e-3; eMax*=1.0e-3;} + + if (ordering == "NEST") { + if (Order() != hporder || nSpectra != fSpectra.size() || Scheme() != NEST ) { + Resize(hporder, eMin, eMax, NEST); + } else { + setSpectra(eMin, eMax); + } + } else { + if (Order() != hporder || nSpectra != fSpectra.size() || Scheme() != RING ) { + Resize(hporder, eMin, eMax, RING); + } else { + setSpectra(eMin, eMax); + } + } + foundEnExt=true; + } catch (CCfits::FITS::NoSuchHDU) { } + } + + //If we found an energy extension, read in the data + if (foundEnExt) { + //Read in the data from the table + //Skip the unnecessary overhead of using CCfits this + //time and use cfitsio routines + //Create a map of typeid's to cfitsio datatype + std::map formatMap; + formatMap[typeid(char).name()] = TSBYTE; + formatMap[typeid(short).name()] = TSHORT; + formatMap[typeid(int).name()] = TINT; + formatMap[typeid(long).name()] = TLONG; + formatMap[typeid(float).name()] = TFLOAT; + formatMap[typeid(double).name()] = TDOUBLE; + formatMap[typeid(unsigned char).name()] = TBYTE; + formatMap[typeid(unsigned short).name()] = TUSHORT; + formatMap[typeid(unsigned int).name()] = TUINT; + formatMap[typeid(unsigned long).name()] = TULONG; + //Select the appropriate datatype + int dataType = formatMap[typeid(T).name()]; + + + //Check for skymap format with column vectors + + if (TFIELDS > 1) //AWS20131217 + { + //column format + for (int icol(1); icol <= nSpectra; ++icol) { + std::valarray binMap(Npix()); + skymapTable.column(icol).read(binMap,1,Npix()); + for (size_t j(0); j < binMap.size(); ++j) + fMap[j*nSpectra + (icol-1)] = binMap[j]; + } + } else { + //vector format + //Point the fits file pointer to the correct extension + skymapTable.makeThisCurrent(); + //Get the fits pointer + fitsfile* ff = fits.fitsPointer(); + //Read the data + int status(0), anynul(0); + T null(0); + fits_read_col(ff, dataType, 1, 1, 1, fMap.size(), &null, &fMap[0], &anynul, &status); + if (status != 0) { + fits_report_error(stderr, status); + throw(1); + } + } //AWS20131217 + } else { + std::cerr<<"Not a compatible fits file, did not find an EBOUNDS or ENERGIES extension"< crval(0.0, 3), crpix(1.0, 3), cdelt(1.0, 3); + std::valarray axis(3); + for (int i = 0; i < axes; ++i) { + axis[i] = mapCube.axis(i); + //Seek to the beginning of the stringstream to overwrite old values + ss.str(""); + ss << "CRPIX" << i+1; + try { + mapCube.readKey(ss.str(), crpix[i]); + } catch (CCfits::HDU::NoSuchKeyword) {} //Assume the value is 1 if undefined + //Seek to the beginning of the stringstream to overwrite old values + ss.str(""); + ss << "CDELT" << i+1; + try { + mapCube.readKey(ss.str(), cdelt[i]); + } catch (CCfits::HDU::NoSuchKeyword) { + //Assume whole sky maps and 1 for all + //other axis + if (i == 0) { + cdelt[i] = 360./axis[i]; + } else if (i == 1) { + cdelt[i] = 180./axis[i]; + } + } + //Seek to the beginning of the stringstream to overwrite old values + ss.str(""); + ss << "CRVAL" << i+1; + try { + mapCube.readKey(ss.str(), crval[i]); + } catch (CCfits::HDU::NoSuchKeyword) { + //Assume full sky maps and 0 for everything else + if (i == 0) { + crval[i] = 0 + cdelt[i]/2.; + } else if (i == 1) { + crval[i] = -90 + cdelt[i]/2.; + } + } + } + + //Read the data and resize the skymap. Let the skymap be of RING + //structure, most favorable for convolution. + //The resolution determines the order for the map, which is always + //greater or equal to the resolution of the map + double res = std::min(fabs(cdelt[0]),fabs(cdelt[1])); + if (order == -1) + order = int(log(sqrt(3./utl::kPi)*60/res)/log(2.0)) + 1; //log(nside)/log(2) + 1 log(2)->log(2.0) AWS20080519 + + //If this is a proper mapcube, read in the energies value + std::valarray energies, emin, emax; + try { + //Try energies extension first + CCfits::ExtHDU & energyTable = fits.extension("ENERGIES"); + int nSpectra; + energyTable.readKey("NAXIS2", nSpectra); + energyTable.column(1).read(energies, 1, nSpectra); + } catch (CCfits::FITS::NoSuchHDU) { + try{ + //Then ebounds + CCfits::ExtHDU & energyTable = fits.extension("EBOUNDS"); + int nSpectra; + energyTable.readKey("NAXIS2", nSpectra); + energyTable.column(2).read(emin, 1, nSpectra); + energyTable.column(3).read(emax, 1, nSpectra); + energies.resize(emin.size()); + } catch (CCfits::FITS::NoSuchHDU) {} + } + + //If no energies extension is found, create the values from CRVAL and + //CDELT values + if (energies.size() == 0) { + if (axes < 3) { + energies.resize(1); + energies[0] = 1; + axis[2] = 1; + }else{ + energies.resize(axis[2]); + //Assume CDELT represents logarithmic values + for (int i = 0; i < axis[2]; ++i) { + energies[i] = pow(10,crval[2] + i*cdelt[2]); + } + } + } + + //Now we can set up the map + if (emin.size() == 0){ + if (Order() != order || energies.size() != fSpectra.size() || Scheme() != RING ) { + Resize(order, energies, RING); + }else{ + setSpectra(energies); + } + } else { + if (Order() != order || emin.size() != fSpectra.size() || Scheme() != RING ) { + Resize(order, emin, emax, RING); + }else{ + setSpectra(emin, emax); + } + } + + //Read the image data and fill the skymap + std::valarray image; + long npix(1); + for (int i = 0; i < axes; ++i){ + npix *= axis[i]; + } + mapCube.read(image,1,npix); + fillMapcube(image, axis, crval, cdelt, crpix); + std::cout<<"All done"< formatMap; + formatMap[typeid(char).name()] = "S"; + formatMap[typeid(short).name()] = "I"; + formatMap[typeid(int).name()] = "J"; + formatMap[typeid(long).name()] = "K"; + formatMap[typeid(float).name()] = "E"; + formatMap[typeid(double).name()] = "D"; + formatMap[typeid(unsigned char).name()] = "B"; + formatMap[typeid(unsigned short).name()] = "U"; + formatMap[typeid(unsigned int).name()] = "V"; + formatMap[typeid(unsigned long).name()] = "K"; + + + //Create the data table to store the actual image. The extension name is + //SKYMAP + + // new format is column per energy + + string format="column"; //AWS20140127 + // format="vector"; //AWS20140127 + + + + std::ostringstream strForm; + strForm << fSpectra.size() << formatMap[typeid(T).name()]; + std::vector colNames(1,"Spectra"), colForm(1,strForm.str()), colUnit(1,""); + + if(format=="column") //AWS20140128 + { + + std::cout<<"Skymap.h column format: formatting output file "<column("Spectra").write(fMap,1); + }else{ + imgTable->column("Spectra").write(fMap,Npix(),1); + } + + } //format=="vector" + + + if(format=="column") //AWS20140129 + { + std::cout<<"Skymap.h column format: writing data to file "<columndata; + columndata.resize(Npix()); + + for (int icol=0;icoladdKey("PIXTYPE", "HEALPIX", "Healpix pixel scheme"); + std::string keyvalue = "RING"; + if (NEST == Order()) keyvalue = "NEST"; + imgTable->addKey("ORDERING", keyvalue, "Ring or nested ordering of pixels"); + imgTable->addKey("NSIDE", Nside(), "Number of sides in a base pixel"); + imgTable->addKey("FIRSTPIX", 0, "Number of first pixel"); + imgTable->addKey("LASTPIX", Npix()-1, "Number of last pixel"); + imgTable->addKey("NBRBINS", nSpectra(), "Number of energy bins in spectra"); + imgTable->addKey("EMIN", fSpectra[0], "Minimum energy of spectra (MeV)"); + double keyv = 0; + if (fSpectra.size() > 1){ + keyv = log(fSpectra[1]/fSpectra[0]); + } + imgTable->addKey("EBIN", keyv, "Energy bin size (logarithmic)"); + imgTable->addKey("COORDTYPE", "GAL", ""); + + //Create another table to store the energy + if (fBinned){ + colNames.resize(3); colUnit.resize(3); colForm.resize(3); + colNames[0] = "CHANNEL"; colNames[1] = "E_MIN"; colNames[2] = "E_MAX"; + colUnit[0] = ""; colUnit[1] = "MeV"; colUnit[2] = "MeV"; //AWS20140130 + colForm[0] = "I"; colForm[1] = "D"; colForm[2] = "D"; + CCfits::Table *eTable = fits.addTable("EBOUNDS", nSpectra(), colNames, colForm, colUnit); + //Create the channel array + std::valarray channel(fSpectra.size()); + for (int i = 0; i < int(channel.size()); ++i){ + channel[i] = i+1; + } + eTable->column("CHANNEL").write(channel,1); + eTable->column("E_MIN").write(fSpecMin,1); + eTable->column("E_MAX").write(fSpecMax,1); + }else{ + colNames.resize(1); colUnit.resize(1); colForm.resize(1); //AWS20140130 + colNames[0] = "Energy"; + colUnit[0] = "MeV"; //AWS20140130 + colForm[0] = "D"; + CCfits::Table *eTable = fits.addTable("ENERGIES", nSpectra(), colNames, colForm, colUnit); + + eTable->column("Energy").write(fSpectra,1); + } + } + + /**\brief Convert all the pixels to a different type and returns the new copy of the map + * + * \param dummy controls the output type + */ + template + Skymap convert(C dummy) const{ + //Create an output, the same size and same spectra + Skymap output; + if (fBinned) { + output.Resize(order_, fSpecMin, fSpecMax, scheme_); + } else { + output.Resize(order_, fSpectra, scheme_); + } + for (size_t i = 0; i < Npix(); ++i){ + for (size_t j = 0; j < nSpectra(); ++j){ + output[i][j] = C((*this)[i][j]); + } + } + return output; + } + + /** \brief Interpolate the map to a different order. + * + * \param order specifies the new order. It makes little sense + * to interpolate to lower order, but it is not forbidden. + * + * Internally we use the get_interpol method of Healpix_Base to + * calculate the interpolation value at the centers of the new + * pixels + */ + Skymap interpolate(int order) const{ + //If the order is the same, just return + if (order == order_) return *this; + + //Create a new skymap with the new order, keeping the spectra and the + //scheme_ + Skymap newMap; + if (fBinned) { + newMap.Resize(order, fSpecMin, fSpecMax, scheme_); + }else{ + newMap.Resize(order, fSpectra, scheme_); + } + + //Loop the new map and calculate the interpolation + fix_arr pixels; + fix_arr weight; + for (int i = 0; i < newMap.Npix(); ++i) { + get_interpol(newMap.pix2ang(i),pixels,weight); + for (int j = 0; j < 4; ++j){ + for (int k = 0; k < fSpectra.size(); ++k){ + newMap[i][k] += weight[j]*(*this)[pixels[j]][k]; + } + } + } + return newMap; + } + + + + /** \brief Rebin the map to a different order. + * + * \param order specifies the new order. If it is greater than the old + * order, the map size is increased and the new map contains multiple + * pixels with the same values. If it is less than the old order, the new + * pixel values will be the average of the old pixels contained in the new + * set. + * + * It is not wise to rebin the map to a higher order, rather use the + * coordinate selector to get the pixel values. That is what this + * rebinning does anyway. + */ + Skymap rebin(int order, bool SAcorrect = true) const{ + //If the order is the same, just return + if (order == order_) return *this; + + //Create a new skymap with the new order, keeping the spectra and the + //scheme_ + Skymap newMap; + if (fBinned) { + newMap.Resize(order, fSpecMin, fSpecMax, scheme_); + }else{ + newMap.Resize(order, fSpectra, scheme_); + } + + //What we do now depends if the order, if it is greater, loop the new map + //and set the pixel from their coordinate + if (order > order_){ + int dOrder = order - order_; + const int nPix = (1< avStore(0.0, fSpectra.size()); + ArraySlice average(&avStore[0], fSpectra.size()); +#pragma omp for schedule(static) + for (long p = 0; p < newMap.Npix(); ++p){ + average = 0; + //Use pixel numbering from the nested scheme + const int pp = (newMap.Scheme() == NEST) ? p : newMap.ring2nest(p); + for (long po = pp*nPix; po < (pp+1)*nPix; ++po){ + //Convert pixel back to correct scheme if needed + const int pop = (Scheme() == NEST) ? po : nest2ring(po); + average += (*this)[pop]; + } + //Divide by the number of pixels to get the average + if (SAcorrect) average /= double(nPix); + newMap[p] = average; + } + }//End parallel + } + //Assign the newMap to us + return newMap; + } + + /**\brief The size of the spectra */ + int nSpectra() const{ return fSpectra.size();} + + /**\brief Set the values at which the spectra is evaluated + * + * \param spectra are the new values, stored in a valarray + * + * \return if the size of the input values do not confirm with the spectral size, a false value is returned. + * + * This method changes the skymap to unbinned mode + */ + bool setSpectra(const std::valarray &spectra){ + if (spectra.size() == fSpectra.size()){ + fSpectra = spectra; + fBinned = false; + return true; + } + return false; + } + + /**\brief Set the values at which the spectra is evaluated in + * binned mode + * + * \param specMin are the new lower boundaries + * \param specMax are the new upper boundaries + * + * \return if the size of the input values do not confirm with the spectral size, a false value is returned. + * + * This method changes the skymap to binned mode + */ + bool setSpectra(const std::valarray &specMin, const std::valarray &specMax){ + if (specMin.size() == fSpectra.size() && specMax.size() == fSpectra.size()){ + fSpecMin.resize(fSpectra.size()); + fSpecMax.resize(fSpectra.size()); + fSpectra = 0.5*(specMin+specMax); + fSpecMin = specMin; + fSpecMax = specMax; + fBinned = true; + return true; + } + return false; + } + + /**\brief Set the values at which the spectra is evaluated + * + * \param spectra are the new values in a C array + * \param size is the size of the new array. + * + * \return if the size of the input values do not confirm with the spectral size, a false value is returned. + */ + bool setSpectra(const double spectra[], const int size){ + if (size == fSpectra.size()){ + for (int i=0; i & getSpectra() const{return fSpectra;} + + /**\brief Return the boundaries for binned maps. + * + * \return false if the map is not binned + */ + bool getBoundaries(std::valarray & specMin, std::valarray & specMax) const{ + if (fBinned) { + specMin.resize(fSpecMin.size()); + specMax.resize(fSpecMax.size()); + specMin = fSpecMin; + specMax = fSpecMax; + } + return fBinned; + } + + /**\brief Fill skymap with data from an l and b map + * + * \param map is an one dimensional array formatted in mapcube format (order l,b,spectra) + * \param nl the number of l bins + * \param nb the number of b bins + * \param correctSA is a boolean which controls wether we correct for solid angle or not. + * + * The spectral size is assumed to be the same as the skymap and no checking is performed. + * This method splits the healpix map into finer bins and selects the undelying l and b pixel. + * The input map must be full sky and the first pixel in the input map must have the edge at + * l = 0 and b = -90. It is assumed that the pixel values are given at the center of the pixels, + * so the pixel (ib, il) spans the area l = [il*360/nl, (il+1)*360/nl] and b = [ib*180/nl-90, (ib+1)*180/nl-90]. + * Note that ib and il are numbered from 0. + */ + void filllbCARarray( const T map[], const int nl, const int nb, const bool correctSA=false){ + //Set up the crval, cdelt, crpix and pass it on + std::valarray crval(0.0, 3), cdelt(1.0,3), crpix(1.0,3); + std::valarray axis(3); + + //Full skymaps + cdelt[0] = 360./nl; + cdelt[1] = 180./nb; + crval[0] = cdelt[0]/2.; + crval[1] = -90+cdelt[1]/2.; + axis[0] = nl; + axis[1] = nb; + axis[2] = nSpectra(); + + //The valarray for the image data + std::valarray image(map,axis[0]*axis[1]*axis[2]); + + fillMapcube(image,axis,crval,cdelt,crpix,correctSA); + } + + /**\bried Convert to a mapcube in CAR projection in galactic + * coordinates and write it to a file + * + * Spatial coordinates are given as in the fits standard, with + * vectors of length 2: NAXIS, CRVAL, CDELT, CRPIX. The first + * element corresponds to longitude but the second to latitude. + * + * If NAXIS is 0, full sky is assumed, trying to take CDELT + * and CRVAL into account. If CDELT is also 0, the resolution + * of the skymap is used. We put a pixel at 0,0 by default. + * + * \param fileName is the name of the output file, with full + * path. + * \param NAXIS is the number of pixels per axis (length 2) + * \param CRVAL is the value of pixel CRPIX (length 2) + * \param CDELT is the difference in value between two pixels + * (length 2) + * \param CRPIX is the pixel number that is at CRVAL (note that + * pixels are numbered from 1 like when you count, not from 0 + * like in C/C++) + * \param correctSA should be true if you are working with + * counts, false if you are working with flux. + * + * \return a 1D valarray that contains all the pixels, where l + * loops fastest, then b and last the energy. + * + * Note that if the length of NAXIS et al. is greater than 2, + * the values are simply ignored. There is no way to rebin in + * energy. On output, the arrays will be resized to length 3 + * and the corresponding value for the energy dimension added + * to it. It will be in log energy. + * + * This routine works by splitting the output pixels up into + * 1024 smaller pixels and assigning it the value that is + * directly under the center. Then the final bigger pixels is + * achieved by taking the average. + */ + void writeMapcube(const std::string &fileName, std::vector &NAXIS, std::vector &CRVAL, std::vector &CDELT, std::vector &CRPIX, bool correctSA = false) const{ + //Check for consistency in the input parameters + NAXIS.resize(2,0); + CDELT.resize(2,0.0); + CRVAL.resize(2,0.0); + CRPIX.resize(2,1.0); + double res = resolution(); + if (NAXIS[0] == 0) { + if (CDELT[0] == 0) { + NAXIS[0] = long(360/res); + CDELT[0] = 360./double(NAXIS[0]); + } else { + NAXIS[0] = long(360/CDELT[0]); + } + CRVAL[0] = 0; + CRPIX[0] = (NAXIS[0]+1)/2; + } + if (NAXIS[1] == 0) { + if (CDELT[1] == 0) { + NAXIS[1] = long(180/res); + } else { + NAXIS[1] = long(180/CDELT[1]); + } + if (NAXIS[1] % 2) { + CDELT[1] = 180./double(NAXIS[1]-1); + } else { + CDELT[1] = 180./double(NAXIS[1]); + ++NAXIS[1]; + } + CRVAL[1] = 0; + CRPIX[1] = (NAXIS[1]+1)/2; + } + //Fix CDELT, if it is still 0 + if (CDELT[0] == 0) { + CDELT[0] = 360./double(NAXIS[0]); + } + if (CDELT[1] == 0) { + CDELT[1] = 180./double(NAXIS[1]); + } + + //Do the conversion + std::valarray mapCube = toMapcube(NAXIS,CRVAL,CDELT,CRPIX,correctSA); + + //Create the fits file + CCfits::FITS fits(fileName, DOUBLE_IMG, 3, &NAXIS[0]); + + fits.pHDU().write(1, NAXIS[0]*NAXIS[1]*NAXIS[2], mapCube); + + //Write the keywords + fits.pHDU().addKey("CRVAL1", CRVAL[0], "Value of longitude in pixel CRPIX1"); + fits.pHDU().addKey("CDELT1", CDELT[0], "Step size in longitude"); + fits.pHDU().addKey("CRPIX1", CRPIX[0], "Pixel that has value CRVAL1"); + fits.pHDU().addKey("CTYPE1", "GLON-CAR", "The type of parameter 1 (Galactic longitude in CAR projection)"); + fits.pHDU().addKey("CUNIT1", "deg", "The unit of parameter 1"); + fits.pHDU().addKey("CRVAL2", CRVAL[1], "Value of latitude in pixel CRPIX2"); + fits.pHDU().addKey("CDELT2", CDELT[1], "Step size in latitude"); + fits.pHDU().addKey("CRPIX2", CRPIX[1], "Pixel that has value CRVAL2"); + fits.pHDU().addKey("CTYPE2", "GLAT-CAR", "The type of parameter 2 (Galactic latitude in CAR projection)"); + fits.pHDU().addKey("CUNIT2", "deg", "The unit of parameter 2"); + fits.pHDU().addKey("CRVAL3", CRVAL[2], "Energy of pixel CRPIX3"); + fits.pHDU().addKey("CDELT3", CDELT[2], "log10 of step size in energy (if it is logarithmically distributed, see ENERGIES/EBOUNDS extension)"); + fits.pHDU().addKey("CRPIX3", CRPIX[2], "Pixel that has value CRVAL3"); + fits.pHDU().addKey("CTYPE3", "Energy", "Axis 3 is the spectra"); + fits.pHDU().addKey("CUNIT3", "MeV", "The unit of axis 3"); + + //Write the energy extension + if (fBinned){ + std::vector colNames(3), colForm(3), colUnit(3); + colNames[0] = "CHANNEL"; colNames[1] = "E_MIN"; colNames[2] = "E_MAX"; + colUnit[0] = ""; colUnit[1] = ""; colUnit[2] = ""; + colForm[0] = "I"; colForm[1] = "D"; colForm[2] = "D"; + std::valarray channel(nSpectra()); + for (int i = 0; i < channel.size(); ++i){ + channel[i] = i+1; + } + CCfits::Table *eTable = fits.addTable("EBOUNDS", nSpectra(), colNames, colForm, colUnit); + eTable->column("CHANNEL").write(channel,1); + eTable->column("E_MIN").write(fSpecMin,1); + eTable->column("E_MAX").write(fSpecMax,1); + }else{ + std::vector colNames(1), colForm(1), colUnit(1); + colNames[0] = "Energy"; + colUnit[0] = ""; + colForm[0] = "D"; + CCfits::Table *eTable = fits.addTable("ENERGIES", nSpectra(), colNames, colForm, colUnit); + + eTable->column("Energy").write(fSpectra,1); + } + + } + + /**\bried Convert to a mapcube in CAR projection in galactic + * coordinates + * + * Spatial coordinates are given as in the fits standard, with + * vectors of length 2: NAXIS, CRVAL, CDELT, CRPIX. The first + * element corresponds to longitude but the second to latitude. + * No effort is put into correcting for input errors. + * + * \param NAXIS is the number of pixels per axis (length 2) + * \param CRVAL is the value of pixel CRPIX (length 2) + * \param CDELT is the difference in value between two pixels + * (length 2) + * \param CRPIX is the pixel number that is at CRVAL (note that + * pixels are numbered from 1 like when you count, not from 0 + * like in C/C++) + * \param correctSA should be true if you are working with + * counts, false if you are working with flux. + * + * \return a 1D valarray that contains all the pixels, where l + * loops fastest, then b and last the energy. + * + * Note that if the length of NAXIS et al. is greater than 2, + * the values are simply ignored. There is no way to rebin in + * energy. On output, the arrays will be resized to length 3 + * and the corresponding value for the energy dimension added + * to it. It will be in log energy. + * + * This routine works by splitting the output pixels up into + * 1024 smaller pixels and assigning it the value that is + * directly under the center. Then the final bigger pixels is + * achieved by taking the average. + */ + std::valarray toMapcube(std::vector &NAXIS, std::vector &CRVAL, std::vector &CDELT, std::vector &CRPIX, bool correctSA = false) const{ + //Create the output array + std::valarray output(NAXIS[0]*NAXIS[1]*fSpectra.size()); + + //Calculating the finer resolution, taking the resolution of + //the current HEALPix map into account, not making it too + //small + const double res = resolution(); + const int nfiner = 32; //The nominal number of subdivision, unless the pixel size of the HEALPix map is small enough + const int nlfiner = CDELT[0] < res ? int(CDELT[0]/res*nfiner)+1 : nfiner; + const int nbfiner = CDELT[1] < res ? int(CDELT[1]/res*nfiner)+1 : nfiner; + int npfiner = nlfiner*nbfiner; + + //The finer resolution + double cdeltl = CDELT[0]/nlfiner; + double cdeltb = CDELT[0]/nbfiner; + + //Add the spectral information to the vectors + NAXIS.resize(3); + CRVAL.resize(3); + CDELT.resize(3); + CRPIX.resize(3); + NAXIS[2] = fSpectra.size(); + CRVAL[2] = fSpectra[0]; + CRPIX[2] = 1; + CDELT[2] = fSpectra.size() > 1 ? log10(fSpectra[1]/fSpectra[0]) : 0; + + //Loop over the output pixels, dividing them up into the + //finer grid +#pragma omp parallel for default(shared) schedule(static) private(npfiner) + for (long ib = 0; ib < NAXIS[1]; ++ib){ + //Create storage to calculate the average + std::valarray totStore(T(0),fSpectra.size()); + ArraySlice totSpectra(&totStore[0], fSpectra.size()); + const size_t indb = ib*NAXIS[0]; + //The value of b at the edge of the pixel + const double b = CRVAL[1] + CDELT[1]*(ib+0.5-CRPIX[1]); + //Correct for solid angle if necessary + double lbSolidAngle = 1; + if (correctSA){ + lbSolidAngle = CDELT[0]*utl::kPi/180. * (sin(utl::kPi*(b+CDELT[1])/180.) - sin(utl::kPi*b/180.)); + } + for (long il = 0; il < NAXIS[0]; ++il){ + totSpectra = 0; + const size_t indl = indb + il; + //The value of l at the edge of the pixel + const double l = CRVAL[0] + CDELT[0]*(il+0.5-CRPIX[0]); + //Loop over the finer grid and calculate the sum + npfiner = 0; + for (int ibf = 0; ibf < nbfiner; ++ibf){ + const double bf = b + (ibf+0.5)*cdeltb; + if (fabs(bf) >= 90) continue; + for (int ilf = 0; ilf < nlfiner; ++ilf){ + double lf = l + (ilf+0.5)*cdeltl; + while (lf < 0) lf += 360; + while (lf >= 360) lf -= 360; + totSpectra += (*this)[SM::Coordinate(lf,bf)]; + ++npfiner; + } + } + //Now assign the average to the pixels + for (int j = 0; j < totSpectra.size(); ++j){ + output[j*NAXIS[0]*NAXIS[1]+indl] = totSpectra[j]/double(npfiner)*lbSolidAngle; + } + } + } + + //Correct for solid Angle if needed + if (correctSA) + output/=solidAngle(); + + return output; + } + + /**\brief Convert to a l and b map in CAR projection + * + * \param nl is the number of l values in the output array + * \param nb is the number of b values in the output array + * \param correctSA is a boolean which tells the routine weather to correct for solid angle or not. + * (When converting counts map, correctSA should be true.) + * + * \return a pointer to a dynamically allocated one dimensional c array compatible with mapcubes. The + * l values are looped fastest, then b values and spectra last. Please deallocate the array once done using + * it to prevent memory leaks. + */ + T * lbCARarray( const int nl, const int nb, double CRVAL[3], int CRPIX[3], double CDELT[3], const bool correctSA = false) const{ + T *array; + array = new T[nl*nb*fSpectra.size()]; + const int nfiner = 32; + const double lres = 360./nl; + const double bres = 180./nb; + const double mapres = resolution(); + //The resolution of the conversion has to be 1/10th of the resolution of + //the resulting map. But it only has to have 1/10th of the resolution of + //the healpix map. + const int nlfiner = lres < mapres ? int(lres/mapres*nfiner)+1 : nfiner; + const int nbfiner = bres < mapres ? int(bres/mapres*nfiner)+1 : nfiner; + const int nboxes = nbfiner*nlfiner; + //Populate the fits header arrays + CRVAL[0] = 0; CRVAL[1] = -90+bres/2.; CRVAL[2] = fSpectra[0]; + CRPIX[0] = CRPIX[1] = CRPIX[2] = 1; //fits starts counting at 1 + CDELT[0] = lres; CDELT[1] = bres; CDELT[2] = fSpectra.size() > 1 ? fSpectra[1] - fSpectra[0] : 0; + //Divide by the solid angle if needed to get the correct results + double hpSolidAngle = 1; + if (correctSA){ + hpSolidAngle = solidAngle(); + } + //Create storage + std::valarray totStore(T(0),fSpectra.size()); + ArraySlice totSpectra(&totStore[0], fSpectra.size()); +#pragma omp parallel for default(shared) schedule(static) + for (int j = 0; j & inmap){ + //Check for iSpectra bounds and the order of the map + if (iSpectra < fSpectra.size() && order_ == inmap.Order() && scheme_ == inmap.Scheme()){ + for(int i=0; i toHealpixMap(const int iSpectra) const{ + //Create the output map and fill it with zeros + Healpix_Map outmap(order_, scheme_); + outmap.fill(0.0); + + //iSpectra must be within bounds, otherwise return a zero map + if (iSpectra < fSpectra.size()){ + for(int i=0; i operator [] (const SM::Coordinate & coordinate){ + return (*this)[ang2pix(coordinate.healpixAng())]; + } + /** \overload */ + const ArraySlice operator [] (const SM::Coordinate & coordinate) const{ + return (*this)[ang2pix(coordinate.healpixAng())]; + } + /**\brief Reference to a spectra given a pixel + * + * \param pixel is a pixel number for the skymap + * \return a std::valarray of the spectra in the pixel + */ + ArraySlice operator [] (const int pixel){ + return ArraySlice(&fMap[pixel*fSpectra.size()], fSpectra.size()); + } + /** \overload */ + const ArraySlice operator [] (const int pixel) const{ + return ArraySlice(&fMap[pixel*fSpectra.size()], fSpectra.size()); + } + + + /** \brief Bidirectional iterator for the skymap */ + class Iterator : public std::iterator > { + private: + int m_index; //!< The index for the array + Skymap & m_map; //!< A reference to a skymap + Iterator(){} //!< The default iterator constructor should not be allowed + public: + /**\brief Constructor */ + Iterator(const int index, Skymap & map) : m_index(index), m_map(map) {}; + /**\brief Pre-increment operator */ + Iterator & operator ++ () { ++m_index; return *this; } + /**\brief Post-increment operator */ + Iterator operator ++ (int i) { Iterator tmp(*this); ++m_index; return tmp; } + /**\brief Pre-decrement operator */ + Iterator & operator -- () { --m_index; return *this; } + /**\brief Post-decrement operator */ + Iterator operator -- (int i) { Iterator tmp(*this); --m_index; return tmp; } + /**\brief Dereference operator */ + ArraySlice operator * () { return m_map[m_index];} + /**\brief Unequal operator */ + bool operator != (const Iterator compare) { return m_index != compare.m_index; } + /**\brief Return the corresponding coordinate */ + SM::Coordinate coord () { SM::Coordinate co(m_map.pix2ang(m_index)); return co; } + }; + /** \brief Constant bidirectional iterator for the skymap */ + class constIterator : public std::iterator > { + private: + int m_index; //!< The index for the array + const Skymap & m_map; //!< A reference to a skymap + constIterator(){} //!< The default iterator constructor should not be allowed + public: + /**\brief Constructor*/ + constIterator(const int index, const Skymap & map) : m_index(index), m_map(map) {}; + /**\brief Pre-increment operator */ + constIterator & operator ++ () { ++m_index; return *this; } + /**\brief Post-increment operator */ + constIterator operator ++ (int i) { Iterator tmp(*this); ++m_index; return tmp; } + /**\brief Pre-decrement operator */ + constIterator & operator -- () { --m_index; return *this; } + /**\brief Post-decrement operator */ + constIterator operator -- (int i) { Iterator tmp(*this); --m_index; return tmp; } + /**\brief Dereference operator */ + const ArraySlice operator * () { return m_map[m_index];} + /**\brief Unequal operator */ + bool operator != (const constIterator compare) { return m_index != compare.m_index; } + /**\brief Return the corresponding coordinate */ + SM::Coordinate coord () { SM::Coordinate co(m_map.pix2ang(m_index)); return co; } + }; + + + /**\brief Returns an Iterator for the first pixel */ + Iterator begin() { return Iterator(0, *this);} + /** \overload */ + constIterator begin() const { return constIterator(0, *this);} + /**\brief Returns an Iterator to the last+1 pixel */ + Iterator end() { return Iterator(npix_, *this); } + /** \overload */ + constIterator end() const { return constIterator(npix_, *this); } + +#ifdef HAVE_ASTRO + /**\brief Convert the skymap from one coordinate system to the + * other + * + * \parameter from is the coordinate system the map is + * currently in + * \parameter to is the coordinate system to transfer to + * \parmeter mj is the modified Julian date + * + * \return Skymap converted to the new coordinate system + * + * Uses similar methods as transformation to a flat CAR + * projection, creating a finer grid for the transformation and + * rebinning to the same grid size again. + */ + Skymap CoordTransform(CoordSys from, CoordSys to, double mj){ + //If from and to are the same, just return this + if (from == to) + return *this; + //Create the transformation function. We do it reversed, + //since we take the to coordinate and transfer to from + TransformFunction f(to,from); + //Create the output map and a finer healpix grid to do the + //transformation on + Skymap out(*this); + const int dOrder = std::min(5,13-order_); //Order can not exceed 13 + const int nPix = (1< totSpectra(T(0), fSpectra.size()); + ArraySlice totsl(&totSpectra[0], fSpectra.size()); + //Loop over all of the subpixels + pointing pnt1 = pix2ang(p); + for (int sp = pp*nPix; sp<(pp+1)*nPix; ++sp){ + const pointing pnt = finer.pix2ang(sp); + pointing old; + f(mj, pnt, old); + totsl += (*this)[old]; + } + totsl /= double(nPix); + out[p] = totsl; + } + return out; + } +#endif + + /**\brief Print the skymap + * + * \param os is the output stream to which to print the skymap + * + * First outputs the coordinate and then the spectra. This routine is for debugging of + * small skymaps. The output gets very large very quickly. The SkymapFitsio routines are + * recommended. + */ + void print (std::ostream & os){ + for (int i = 0; i < npix_; ++i){ + SM::Coordinate co(pix2ang(i)); + os << co; + os << std::endl << " "; //!< First output the coordinates + for(int j = 0; j < fSpectra.size(); ++j){ + os << fMap[i][j] << " "; + } + os<= 0 && iSpectra < int(fSpectra.size())){ + m = fMap[iSpectra]; + for (int i = 1; i < npix_; ++i){ + m = std::min(fMap[i*fSpectra.size()+iSpectra], m); + } + } + return m; + } + + /** \brief Returns the maximum value of the map */ + T max () const { + return fMap.max(); + } + /** \brief Returns the maximum value of the map at a given point in the spectra + * + * \param iSpectra is the point in spectra to evaluate the maximum value. + */ + T max(int iSpectra) const { + T m(0); + if (iSpectra >= 0 && iSpectra < int(fSpectra.size())){ + m = fMap[iSpectra]; + for (int i = 1; i < npix_; ++i){ + m = std::max(fMap[i*fSpectra.size()+iSpectra], m); + } + } + return m; + } + + /**\brief Sum of all the pixels for a point in spectra + * + * \param iSpectra is the point in spectra to evaluate the sum + */ + T sum (int iSpectra) const { + T s(0); + //Check if iSpectra is within bounds + if (iSpectra >= 0 && iSpectra < fSpectra.size()){ + for (int i = 0; i < npix_; ++i){ + s += fMap[i*fSpectra.size()+iSpectra]; + } + } + return s; + } + + /**\brief Sum all of the pixels in the map */ + T sum () const { + return fMap.sum(); + } + + /**\brief Assigns a single value to the map */ + template + Skymap & operator = (const C number){ + fMap = T(number); + return(*this); + } + + /**\brief Assignment for maps, basically return the same map */ + Skymap & operator = (const Skymap & oldMap){ + //Avoid self-assignment + if (this != &oldMap){ + if (oldMap.fBinned){ + Resize(oldMap.order_, oldMap.fSpecMin, oldMap.fSpecMax, oldMap.scheme_); + }else{ + Resize(oldMap.order_, oldMap.fSpectra, oldMap.scheme_); + } + fMap = oldMap.fMap; + } + return(*this); + } + + //! Add a number to the current skymap + template + Skymap & operator += (C number) { + fMap += T(number); + return (*this); + } + //! Add a number to a skymap + template + Skymap operator + (C number) const { + Skymap returnMap(*this); + returnMap += number; + return returnMap; + } + + //! Withdraw a number from the current skymap + template + Skymap & operator -= (C number) { + fMap -= number; + return (*this); + } + //! Withdraw a number from a skymap + template + Skymap operator - (C number) const { + Skymap returnMap(*this); + returnMap -= number; + return returnMap; + } + + //! Multiply the current skymap with a number + template + Skymap & operator *= (C number) { + fMap *= number; + return (*this); + } + //! Multiply a skymap with a number + template + Skymap operator * (C number) const { + Skymap returnMap(*this); + returnMap *= number; + return returnMap; + } + + //! Divide the current skymap with a number + template + Skymap & operator /= (C number) { + fMap /= number; + return (*this); + } + //! Divide a skymap with a number + template + Skymap operator / (C number) const { + Skymap returnMap(*this); + returnMap /= number; + return returnMap; + } + + /**\brief Add a skymap to the current one + * + * Only adds the skymap if the values at which the spectra is evaluated are the same and the skymaps are the + * same size. + */ + Skymap & operator += (const Skymap & otherMap){ + if (! equiv(otherMap) ){ + std::cerr<<"Skymaps are not equivalent in addition"< operator + (const Skymap & otherMap) const{ + Skymap returnMap(*this); + returnMap += otherMap; + return (returnMap); + } + + /**\brief Subtract a skymap from the current one + * + * Only subtracts the skymap if the values at which the spectra is evaluated are the same and the skymaps are the + * same size. + */ + Skymap & operator -= (const Skymap & otherMap){ + if (! equiv(otherMap) ){ + std::cerr<<"The skymaps arent equivalent when subtracting them. Check that their spectra and sizes are the same."< operator - (const Skymap & otherMap) const{ + Skymap returnMap(*this); + returnMap -= otherMap; + return (returnMap); + } + + /**\brief Multiply a skymap with the current one + * + * Only multiplies the skymap if the values at which the spectra is evaluated are the same and the skymaps are the + * same size. + */ + Skymap & operator *= (const Skymap & otherMap){ + if (! equiv(otherMap) ){ + std::cerr<<"The skymaps arent equivalent when multiplying them. Check that their spectra and sizes are the same."< operator * (const Skymap & otherMap) const{ + Skymap returnMap(*this); + returnMap *= otherMap; + return (returnMap); + } + + /**\brief Divide a skymap into the current one + * + * Only divides the skymap if the values at which the spectra is evaluated are the same and the skymaps are the + * same size. + */ + Skymap & operator /= (const Skymap & otherMap){ + if (! equiv(otherMap) ){ + std::cerr<<"The skymaps arent equivalent when dividing them. Check that their spectra and sizes are the same."< 1e5) + for (int i = 0; i < fMap.size(); ++i){ + fMap[i] /= otherMap.fMap[i]; + } + return (*this); + } + + /**\brief Divide a skymap with another skymap + * + * Only divides the skymaps if the values at which the spectra is evaluated are the same and the skymaps are the + * same size. + */ + Skymap operator / (const Skymap & otherMap) const{ + Skymap returnMap(*this); + returnMap /= otherMap; + return (returnMap); + } + + + /**\brief Equality operator, returns true if every pixel has the same value and the values at which the spectra + * is evaluated is the same as well. + */ + bool operator == (const Skymap & otherMap) const{ + //Check for equivalence of the map + bool map = equiv(otherMap); + if (! map) return map; + //And the map, if needed + for (int i = 0; i < fMap.size(); ++i){ + map &= (fMap[i] == otherMap.fMap[i]); + if (! map) { + return map; + } + } + return map; + } + + /** \brief Non-equality operator, the inverse of the equality operator. */ + bool operator != (const Skymap & otherMap) const{ + return ( ! ((*this) == otherMap)); + } + +}; + +#endif diff --git a/source/SkymapFitsio.h b/source/SkymapFitsio.h new file mode 100644 index 0000000000000000000000000000000000000000..2b4151b07b88ff630a2f8a10f3bacf046674673f --- /dev/null +++ b/source/SkymapFitsio.h @@ -0,0 +1,37 @@ +#ifndef SkymapFitsIO_h +#define SkymapFitsIO_h + +#include "Skymap.h" +#include + +/** \brief Write skymap to fits file. + * + * \param skymap is the skymap to output + * \param filename is the filename to write to. This routine will overwrite + * files without a warning. + * \param sUnit the units of the spectra to put in the fits file + * \param sType the type of the spectra to put in the fits file + * \return integer value telling the cfitsio return status. A value of 0 means + * all went well + * + * The output format is slightly modified from the standard healpix output. We + * still use a table extension to represent the skymap, but each pixel has its + * own row in the table and there is a single column containing a vector + * element which has the spectra for each pixel. The extension for the healpix + * table is called SKYMAP. In addition, there is a ENERGIES extension, + * containing a table of energies (frequencies) which the spectra is evaluated + * at. + */ +int SkymapToFits (const Skymap & skymap, const std::string & filename, const std::string & sUnit, const std::string & sType); +/** \overload */ +int SkymapToFits (const Skymap & skymap, const std::string & filename, const std::string & sUnit, const std::string & sType); +/** \brief Read a skymap from a fits file. + * + * \param skymap will be replaced with the one from the fits file + * \param filename is the name of the file to read + */ +int FitsToSkymap (Skymap & skymap, const std::string & filename); +/** \overload */ +int FitsToSkymap (Skymap & skymap, const std::string & filename); + +#endif diff --git a/source/SourceModel2.h b/source/SourceModel2.h new file mode 100644 index 0000000000000000000000000000000000000000..da55a171ace61aac3bbe41fc1bfb87b72a3973ab --- /dev/null +++ b/source/SourceModel2.h @@ -0,0 +1,119 @@ +#ifndef SourceModel2_h +#define SourceModel2_h +#include "Model.h" +#include "Exposure.h" +#include "Counts.h" +#include "Psf.h" +#include "Sources.h" +#include "Variables.h" +#include "Skymap.h" +#include "SparseSkymap.h" +#include "Parameters.h" +#include +#include + +/**\brief A model which handles the point sources from the catalog. + * + * The sources to include in the model, and how, is controlled by a special + * column in the source catalog, called FIT_TYPE. It understands the following + * values: + * - \e FIXED Include the source as the power law specified in the catalog. + * Do not fit. + * - \e PREFACTOR Include the source as a power law and vary the prefactor in + * the fit but keep the index at the original value. + * - \e INDEX Include the source as a power law and vary the index in the fit + * but keep the prefactor at the original value. + * - \e POWERLAW Include the source as a power law and vary both the prefactor + * and the index. + * - \e POWERLAWEXPCUTOFF Include the source as a power law with an + * exponential cut off. Vary all three parameters of the model. The initial + * value for the cut off is 10 GeV. + * - \e BROKENPOWERLAW Include the source as a broken power law. Vary all + * four parameters of the model. The initial value for the break energy is + * 10 GeV and the initial value for the second index is the same as the first + * index. The gradient calculation for this model has not been done, so it is + * disabled for now. + * - \e FREE Include the source and vary its spectral profile freely. It uses + * one parameter per energy bin in the fit, so the number of parameters can + * grow quickly if the number of energy bins are not kept minimal. + * + * For every other value in this column, the source is not included in the + * model. + */ +class SourceModel2 : public BaseModel { + public: + /** \brief The constructor takes the same parameters as the BaseModel constructor. + * + * Creates the fixed source maps, as well as the necessary sparse skymaps for sources that need to be fitted. + */ + SourceModel2(const Exposure &exposure, const CountsMap &counts, const Psf &psf, const Sources &sources, const Parameters &pars, const Skymap &filter, const std::string &scalerVariable="FixedSourceScaling", unsigned int configure = 3); + // Destructor needed to delete all the sparseSkymaps + // dynamically allocated + ~SourceModel2(); + void getMap(const Variables &vars, Skymap &map); + gradMap getComponents(const Variables & vars, const std::string &prefix); + void getGrads(const Variables &vars, const std::string & varName, Skymap &map); + gradMap getSecDer(const Variables & vars, const std::string & varName1, const std::string & varName2); + void modifySources(const Variables &vars, Sources &sources) const; + /** \brief Add skymap to fixed sources Skymap + * + * Used when bootstrapping all sources to skip the creation of the fixed + * source map each time. + */ + void addToFixedSources(const Skymap&); + //! Enum for defining the type of the variable + enum varType { + PREF, //!< Prefactor + INDX, //!< Index in a power law + ECUT, //!< Cut off energy + EBRE, //!< Break energy + IND1, //!< First index in a broken power law + IND2, //!< Second index in a broken power law + SCAL, //!< Free scaling per bin + UNKN //!< Unknown variable type + }; + private: + /** \brief Add a count spectra to a SparseSkymap + * + * Adds the convolved count spectra of the source with number + * SourceNumber to counts. + * + * \param SourceNumber is the number of the source to add to + * the map + * \param counts is the sparse skymap which the source is added + * to. + */ + void genCountMap(int SourceNumber, SparseSkymap *counts) const; + /** \brief Add count spectras to a SparseSkymap + */ + void genCountMaps(const std::valarray &sourceNumbers, SparseSkymap *counts) const; + /** \brief Add count spectras to a SparseSkymap + */ + void genCountMaps(int numberOfSources, SparseSkymap *counts) const; + /** \brief Add count spectras to a SparseSkymap + */ + void genCountMaps(double lowFlux, SparseSkymap *counts) const; + /** \brief Find the type of the variable and index of a source given a + * variable name. + * + * \param varName is the name of the variable + * \param index is the source index + * \param type is the varType of the variable + */ + void getIndexAndType(const std::string & varName, int &index, varType &type) const; + std::vector fFixedSources, //!< Store indices of fixed sources + fPrefactorSources, //!< Store indices of sources with prefactor fitting + fIndexSources, //!< Store indices of sources with index fitting + fPowerLawSources, //!< Store indices of sources with power law fitting + fPowerLawExpCutOffSources, //!< Store indices of sources with power law fitting and exponential cutoff + fBrokenPowerLawSources, //!< Store indices of sources with broken power law fitting + fFreeSources; //!< Store indices of sources with one parameter per energy bin + Skymap ffixedSkymap; //!< Cache for the fixed sources. They take a long time to compute + std::map* > fvariableSkymaps; //!< We store the skymaps for the variable sources separately in a sparse skymap + std::map > flastPars; //!< We cache the last known values for the variables to speed up calculations + std::string ffixedScalerVar; //!< Name of the fixed sources skymap scaler variable. Scaler not used if it is an empty string. + double fpsfFraction; //!< The fraction of psf containment to use for the sources. + double fscalerMin; //!< The minimum value for the free source scalers. +}; + +#endif diff --git a/source/SourcePopulation.h b/source/SourcePopulation.h new file mode 100644 index 0000000000000000000000000000000000000000..01c241f4c00104420e7eff37287c158a5aecfe0a --- /dev/null +++ b/source/SourcePopulation.h @@ -0,0 +1,183 @@ + + + +#include"Distribution.h" +#include"Skymap.h" //AWS20090107 +#include //AWS20100817 +#include //AWS20100817 + +class SourcePopulation +{ + public: + // parameters of population + // luminosity in photons s-1 in given energy band + + double pi, rtd, dtr, kpc_to_cm; + + char title[200]; //AWS20140102 was 100 + + int luminosity_function_type;// luminosity-function: 1=delta-function, 2=power-law + double L; // luminosity if delta-function + + // double density0,alpha_R,alpha_z; // parameters of exponential disk + double density0,alpha_R,beta_R,zscale; // parameters of source distribution + double alpha_z; + + double oversample;// source oversampling to avoid statistical fluctuations AWS20060118 + + + double L_min,L_max,alpha_L; // parameters of power-law luminosity function + + int spectral_model; // 1 = 2 breaks fixed, 2 = 2 breaks with dispersion AWS20100824 + + double spectrum_g,spectrum_norm_E; // spectral index and normalization energy of sources + double spectrum_g_0 ;// spectrum with two breaks AWS20051219 + double spectrum_br0 ; + double spectrum_g_1 ; + double spectrum_br1 ; + double spectrum_g_2 ; + + + double spectrum_g_0_sigma ;// spectrum with two breaks: Gaussian distribution AWS20100823 + double spectrum_br0_sigma ; + double spectrum_g_1_sigma ; + double spectrum_br1_sigma ; + double spectrum_g_2_sigma ; + + int n_E_bin; // number of energy bins + double *E_bin_low,*E_bin_high; // boundaries of energy bins + double E_ref_low, E_ref_high; // range for which luminosity is defined + + + int n_E; // number of energies for skymap spectra + double *E; // energies for skymap spectra + + double NL(double L); + double density(double R, double z); // number per kpc^3 + + // derived quantities + long n_sample; + double *sample_R, *sample_x, *sample_y, *sample_z; + double *sample_theta,*sample_d,*sample_l,*sample_b,*sample_flux; + double *sample_flux100,*sample_flux1000,*sample_flux10000; //AWS20120615 integral fluxes + double *sample_flux_spectrum_E_ref_low,*sample_flux_spectrum_E_ref_high; //AWS20111001 + double sample_d_min,sample_d_max,sample_flux_min,sample_flux_max; + double sample_total_flux,sample_total_L; + double sample_selected_total_flux, sample_selected_total_flux_sublimit, sample_selected_total_flux_soplimit; + double sample_selected_sourcecnt , sample_selected_sourcecnt_sublimit, sample_selected_sourcecnt_soplimit; + double sample_selected_d_min, sample_selected_d_max, sample_selected_d_av; + double sample_selected_d_min_sublimit,sample_selected_d_max_sublimit,sample_selected_d_av_sublimit; + double sample_selected_d_min_soplimit,sample_selected_d_max_soplimit,sample_selected_d_av_soplimit; + double *sample_l_min,*sample_b_min, *sample_l_max,*sample_b_max, *sample_l_av,*sample_b_av; + + + double **sample_flux_spectra, **sample_binned_flux_spectra; //AWS20100817 + double *sample_spectrum_g_0,*sample_spectrum_br0,*sample_spectrum_g_1,*sample_spectrum_br1,*sample_spectrum_g_2; //AWS20100823 parameters per source + + double latprof_l_min, latprof_l_max, latprof_dlat ; // for profiles + double longprof_b_min,longprof_b_max,longprof_dlong; + + double long_min1,long_max1; // for logN-logS + double long_min2,long_max2; + double lat_min1, lat_max1; + double lat_min2, lat_max2; + + double *longprof_intensity,*longprof_intensity_sublimit ; + double *longprof_sourcecnt,*longprof_sourcecnt_sublimit ; + int longprof_nlong; + double * latprof_intensity,* latprof_intensity_sublimit ; + double * latprof_sourcecnt,* latprof_sourcecnt_sublimit ; + int latprof_nlat; + + double flux_detection_limit; + long i_flux_min,i_flux_max; + + + + double *dlnN_dlnS; + double *dlnN_dlnS_deg; + + double *dlnN_dlnS_int; + double *dlnN_dlnS_int_deg; + + double *FS,*FS_int; + + double lnS_min,lnS_max,dlnS; + double dlnN_dlnS_total_N, dlnN_dlnS_total_S; + int n_dlnN_dlnS; + + + double *dlnN_dlnS_soplimit; // N(S) above flux limit AWS20170123 + double *dlnN_dlnS_sublimit; // N(S) below flux limit AWS20170123 + + double *dlnN_dlnS_int_soplimit; //AWS20170123 + double *dlnN_dlnS_int_sublimit; //AWS20170123 + + double *FS_soplimit,*FS_int_soplimit; //AWS20170123 + double *FS_sublimit,*FS_int_sublimit; //AWS20170123 + + double R0 ;// Sun-Galactic centre in kpc + +// skymaps of source counts N(S) + Distribution skymap_sourcecnt; + Distribution skymap_sourcecnt_sublimit; + Distribution skymap_sourcecnt_soplimit; + +// skymaps integrated over energy ranges + Distribution skymap_intensity; + Distribution skymap_intensity_sublimit; + Distribution skymap_intensity_soplimit; + +// skymaps for full energy spectrum + Distribution skymap_intensity_spectrum; + Distribution skymap_intensity_spectrum_sublimit; + Distribution skymap_intensity_spectrum_soplimit; + +// --- healpix versions of the same skymaps --- + int healpix_skymap_order; // order of healpix maps + +// skymaps of source counts N(S) + Skymap healpix_skymap_sourcecnt; + Skymap healpix_skymap_sourcecnt_sublimit; + Skymap healpix_skymap_sourcecnt_soplimit; + +// skymaps integrated over energy ranges + Skymap healpix_skymap_intensity; + Skymap healpix_skymap_intensity_sublimit; + Skymap healpix_skymap_intensity_soplimit; + +// skymaps for full energy spectrum + Skymap healpix_skymap_intensity_spectrum; + Skymap healpix_skymap_intensity_spectrum_sublimit; + Skymap healpix_skymap_intensity_spectrum_soplimit; + +// Fermi-LAT source sensitivity maps AWS20170112 + Healpix_Map healpix_skymap_Fermi_sensitivity; //AWS20170112 + char* Fermi_sensitivity_file; //AWS20170112 + + int verbose; + + ofstream *txt_stream; + + int init(); + int print(int options=1); + int gen_sample(unsigned seed); + int analyse_sample(); + int gen_sample_EGRET_catalogue(int options); + int gen_sample_SPI_catalogue (char *directory, char *filename, int options); + int gen_sample_Fermi_catalogue(char *directory, char *filename, int options);//AWS20081215 + + int gen_Fermi_catalogue_from_sample(char *directory, char *filename, int options);//AWS20110829 + + int read_sample_MSPs (char *directory, char *filename, int options);//AWS20120521 + + + SourcePopulation operator=(SourcePopulation); // copy constructor + + double integratepowbr(double A0,double g0,double g1,double g2, double br0,double br1, double E1,double E2); + double shapepowbr(double A0,double g0,double g1,double g2, double br0,double br1, double E); + + double flux_cutoff_powerlaw(double E , double g, double Ecutoff,double Emin_ref, double Emax_ref, double flux_ref); + double flux_cutoff_powerlaw(double Emin, double Emax, double g, double Ecutoff,double Emin_ref, double Emax_ref, double flux_ref); + + }; diff --git a/source/Sources.h b/source/Sources.h new file mode 100644 index 0000000000000000000000000000000000000000..5bc0e3b7e8fd82c74674111a2d2db37deff32d9c --- /dev/null +++ b/source/Sources.h @@ -0,0 +1,240 @@ +#ifndef SOURCES_H +#define SOURCES_H + +#include +#include +#include +#include +#include "Skymap.h" +#include "Psf.h" +#include "Exposure.h" + +//The value for E0 in the definition of the source spectra +#define E0 100.0 + +/** \brief Store relevant information from the source catalog. + * + * Has the ability to generate counts map from the source catalog, given an Exposure and a Psf. + */ +class Sources { + public: + /** \brief Default constructor, does nothing */ + Sources() {} + /** \brief Construct the object from a catalog in a fits file + * + * \param fileName is the file name of the fits file containing the catalog. + * + * The file should contain a table in the first extension with at least the + * following columns (formats) + * - SOURCE_NAME (string) the name of the source + * - L (float) the galactic longitude of the source in degrees + * - B (float) the galactic latitude of the source in degrees + * - FLUX100 (float) integrated flux above 100 MeV. Used for sorting on + * brightness + * - PREFACTOR (float) the power law prefactor associated with the source spectra + * - SPECTRAL_INDEX (float) the power law index of the source spectra. The + * spectra is assumed to have a pivot point around 100 MeV + * + * Using the following columns will also work and is + * automatically selected + * - NickName (string) the name of the source + * - GLON (float) galactic longitude in degrees + * - GLAT (float) galactic latitude in degrees + * - Flux100 (float) integrated flux above 1000 MeV + * - Flux_Density (float) the prefactor for the power law + * - Spectral_Index (float) the index for the power law + * - Pivot_Energy (float) the pivot energy for the power law + * (F = F_0 (E/E_p)^g + * + * Additionally, these columns will be read if they exist + * - SP_TYPE (string) the name of the spectral type. Available types are + * - POWERLAW + * - POWERLAWEXPCUTOFF + * - BROKENPOWERLAW + * - FREE + * - SP_PARS (float array) the parameters for the spectra + * - POWERLAW: prefactor, index (100 MeV pivot point) + * - POWERLAWEXPCUTOFF: prefactor, index, cutoff energy (100 MeV pivot point) + * - BROKENPOWERLAW: prefactor, index1, index2, break energy + * - FREE: energy, scaler pairs, where the scaler is the deviation from + * the power law given in the columns before. Power law interpolation + * is used if the energy ranges are not compatible. + * - SP_PARS_UNC (float array) the uncertainty of the spectral + * parameters. Has the same structure as SP_PARS. For the + * free floating spectra, gives the half-width of the energy + * bin used in the fit. + * - FIT_TYPE (string) how to fit the source in SourceModel2. \see + * FitType structure for available values. + * + * The SP_TYPE, SP_PARS, and SP_PARS_UNC columns are ment to be used as output from fits + * performed with SourceModel2, but are of course read in in subsequent + * calls and can be used manually. + * + * Also, if the spectral type is not defined in the catalog, + * the routine tries to read the following columns + * - Flux100_300 + * - Flux300_1000 + * - Flux1000_3000 + * - Flux3000_10000 + * - Flux10000_100000 + * + * which indicate the integrated flux for those bins. + * Subsequently, the spectral type is set to free and these + * values used for the spectra. It also reads the accompaning + * uncertainty columns. + */ + Sources(const std::string &fileName); + /** \brief Enum values for available spectral shapes. + * + * This defines both the shape of the spectra and how we wan't the source to be + * fit. The information is read from the source catalog. The catalog is + * assumed to contain a index and a prefactor for a power law spectral shape, + * which is used as a initial value for the fit. The exponential cut off is + * put at 10 GeV as a default and the break energy in the broken power law likewise. + * + * Currently, the broken power law is not functioning. + */ + enum FitType { + FIXED, //!< Use the values specified in the catalog + PREFACTOR, //!< Fit only the prefactor of the simple power law + INDEX, //!< Fit only the index of the simple power law + POWERLAW, //!< Simple power law with a prefactor and an index + POWERLAWEXPCUTOFF, //!< Simple power law with a high energy cutoff + BROKENPOWERLAW, //!< Broken power law, with two indices, a prefactor and a break energy. The prefactor is the spectral value at the break energy. + FREE, //!< One parameter per energy bin + NONE //!< Do not include the source + }; + + //! Convert FitType to string. + static std::string fitTypeToStr(const FitType & type); + + /** Convert strings to FitType. + * + * Converts strings containing the enum names for FitType into enum + * objects. This routine is case insensitive. + */ + static FitType strToFitType(const std::string & type); + + /** \brief Structure to store the source information */ + struct Source { + std::string name; //!< Name of the source + double l, b, flux100, prefactor, index, pivot; + double flux1000; //AWS20110806 + double flux10000;//AWS20120530 + double ra,dec; //AWS20120612 + + FitType fitType; //!< How we are supposed to fit + FitType spType; //!< What kind of spectra we have from the extended catalog + /** \brief The spectral parameters + * + * The size of this array depends on the spectral Type + * - \e POWERLAW: prefactor, index + * - \e POWERLAWEXPCUTOFF: prefactor, index, cutoff energy + * - \e BROKENPOWERLAW: prefactor, index1, index2, break energy + * - \e FREE: energy, prefactor pairs + * + * In the case of a free spectra, the prefactor is a correction for the power law given in the catalog for the given energy, which is sqrt(Emin*Emax) + */ + std::valarray spParameters; + /** \brief Uncertainty of the spectral parameters + * + * This is the uncertainty for the parameters given in spParameters. + * It has the same shape as the spParameters array. + * For the free fit, the uncertainty in energy is half the bin width. + */ + std::valarray spParametersUnc; + }; + + /** \brief Read a new catalog from file + * + * \see the constructor for details on the catalog format + */ + void read(const std::string &fileName); + /** \brief Write the current catalog to file, including all spectral + * types. + */ + void write(const std::string &fileName); + /** \brief Generate the spectra for a given source + * + * \param sourceNumber is the source for which the spectra is + * to be generated + * \param energies is an array containing the energies at which + * the spectra is to be evaluated. + * + * \return the spectra in units of cm^-2 s^-1 MeV^-1 at the + * required energies. + */ + std::valarray genSpectra(int sourceNumber, const std::valarray &energies) const; + /** \brief Add fluxes of sources brighter than a limit to a skymap + * + * \param lowFlux is the minimum flux >100MeV a source has to have to be included. It is given in photon/cm^2/s + * \param fluxMap is the skymap to add the sources to in + * units of photons/cm^2/MeV/s + * + * Does not take the point spread function into account. + */ + void genFluxMap(double lowFlux, Skymap &fluxMap) const; + /** \brief Add fluxes of numberOfSources brightest + * sources to a skymap + * + * \param numberOfSources is the number of sources to include + * \param fluxMap is the skymap to add the sources to in + * units of photons/cm^2/MeV/s + * + * Does not take the point spread function into account. + */ + void genFluxMap(int numberOfSources, Skymap &fluxMap) const; + /** \brief Return a filter map for the sources. + * + * Uses the fraction of the psf as the region to filter. + */ + void genFilterMaps(double lowFlux, Skymap &filter, double fraction, const Psf &psf) const; + /** overload */ + void genFilterMaps(int numberOfSources, Skymap &filter, double fraction, const Psf &psf) const; + /** overload */ + void genFilterMaps(const std::valarray &sourceNumbers, Skymap &filter, double fraction, const Psf &psf) const; + /** overload */ + void genFilterMap(int sourceNumber, Skymap &filter, double fraction, const Psf &psf) const; + /** \brief Return a constant reference to a source */ + const Source & operator [] (int i) const { return fSources[i]; } + /** \brief Return a reference to a source */ + Source & operator [] (int i) { return fSources[i]; } + //! The number of sources + int size() const { return fSources.size(); } + //! Given a name, find the source number + int sourceNumber(const std::string & name) const; + //! Assignment operator + Sources & operator = (const Sources &); + /** \brief Function to convert fluxes in energy bins into + * intensities at the geometrical mean energy of the bin + * + * \parameter eMin, eMax are the boundaries of the bins, must + * be in ascending order, but not necessarily adjacent + * \parameter flux are the flux values of the bins + * \return the intensities at the geometrical + * mean energy of the bins + */ + static std::valarray flux2int(const std::valarray &eMin, const std::valarray &eMax, const std::valarray &flux); + + int add_source (Source source); //AWS20110829 + + private: + struct fparams { + std::valarray eMin, eMax, flux; + }; + static int minFunc(const gsl_vector *x, void *params, gsl_vector *f); + /** \brief Sort the sources on their flux > 100MeV value */ + void sort(); + /** \brief Read the source catalog from fits file */ + void readFromFile(const std::string & fileName); + std::vector fSources; //!< Store the sources + std::map fnameIndex; //!< Store the name, index in a map for fast name lookup + bool foldFormat; //!< We must know what format we are using + /** \brief Helper class to compare the sources. Needed in the sort method */ + class SourceCompare { + public: + bool operator () (const Source &a, const Source &b); + }; +}; + +#endif diff --git a/source/SparseSkymap.h b/source/SparseSkymap.h new file mode 100644 index 0000000000000000000000000000000000000000..4b802b7586828e91af90f8624f2a0b83ece7ccef --- /dev/null +++ b/source/SparseSkymap.h @@ -0,0 +1,251 @@ +#ifndef SparseSkymap_h +#define SparseSkymap_h + +#include "Skymap.h" +#include +#include + +typedef std::pair indexType; //!< This is the index for each value in the sparse map, the pixel index and the spectral index + +/** \brief Store only set values of pixels. + * + * This class was built for counts map for sources, where the psf excludes any emission in a large portion of the map. + * Stores pixel and spectra indices as well as the pixel value. + * + * Only possible to initialize from Skymap object, where 0 values are eliminated. + */ +template +class SparseSkymap : public HealpixBaseExtended, public std::map{ + private: + std::valarray fSpectra, fSpecMin, fSpecMax; //!< The values at which the spectra is evaluated + bool fBinned; + typedef std::pair valueType; //!< Typedef for the pixel storage + + /** \brief Helper class to handle copying mapTypes */ + class copyFunctor{ + private: + SparseSkymap *fmap; + public: + copyFunctor( SparseSkymap * map) : fmap(map) {} + void operator () (const valueType & pixel){ + fmap->insert(pixel); + } + }; + /** \brief Helper class to turn mapTypes to Skymaps */ + class fillFunctor{ + private: + Skymap &fmap; //! &map): fmap(map){} + void operator () (const valueType & pixel){ + fmap[pixel.first.first][pixel.first.second] = pixel.second; + } + }; + + + public: + /** \brief Default constructor, sets size to 0 */ + SparseSkymap() : fBinned(false) {Resize(0,0);} + /** \brief Construct a map of given order and spectral size + * + * \param order is the healpix order of the map + * \param nSpectra is the spectral size of the map + * \param scheme is the healpix ordering scheme (default value is RING) + */ + SparseSkymap(int order, int nSpectra, Healpix_Ordering_Scheme scheme=RING){ + Resize(order,nSpectra,scheme); + } + /** \brief Construct a map of given order and values at which spectra is evaluated + * + * \param order is the healpix order of the map + * \param spectra is the values at which the spectra is evaluated + * \param scheme is the healpix ordering scheme (default value is RING) + */ + SparseSkymap(int order, const std::valarray & spectra, Healpix_Ordering_Scheme scheme=RING){ + Resize(order,spectra,scheme); + } + + /** \brief Construct a map of given order and values at which spectra is evaluated + * + * \param order is the healpix order of the map + * \param specMin are the lower boundaries for binned maps + * \param specMax are the lower boundaries for binned maps + * \param scheme is the healpix ordering scheme (default value is RING) + */ + SparseSkymap(int order, const std::valarray & specMin, const std::valarray & specMax, Healpix_Ordering_Scheme scheme=RING){ + Resize(order,specMin,specMax,scheme); + } + + /** \brief Construct a map from a Skymap object. + * + * \param Skymap is the skymap to convert to sparse skymap. All 0 value + * pixels are removed. + */ + SparseSkymap(const Skymap & map){ + fromSkymap(map); + } + + //Copy constructor + /* Not working at the moment, don't know why, some error about too many + * arguments to function from for_each + SparseSkymap(const SparseSkymap & oldMap){ + Resize(oldMap.order_, oldMap.fSpectra, oldMap.scheme_); + //Copy all the pixels from the old map + copyFunctor f(); + std::for_each(oldMap.begin(), oldMap.end(), f); + }*/ + + /** \brief Resize the map + * + * \param order is the new healpix order + * \param nSpectra is the new spectral size for the map + * \param scheme is the new healpix ordering scheme. + * + * \note This does not clear the data in the map and it could be rendered invalid with this method + * if the order is changed. + */ + void Resize(int order, int nSpectra, Healpix_Ordering_Scheme scheme=RING){ + std::valarray spectra(1.0,nSpectra); + Resize(order,spectra,scheme); + } + /** \brief Resize the map + * + * \param order is the new healpix order + * \param spectra is the new values at which the spectra is evaluated + * \param scheme is the new healpix ordering scheme. + * + * \note This does not clear the data in the map and it could be rendered invalid with this method. + */ + void Resize(int order, const std::valarray & spectra, Healpix_Ordering_Scheme scheme=RING){ + Set(order,scheme); + fSpectra.resize(spectra.size()); + fSpectra=spectra; + fBinned=false; + } + /** \brief Resize the map + * + * \param order is the new healpix order + * \param specMin are the lower boundaries for binned maps + * \param specMax are the lower boundaries for binned maps + * \param scheme is the new healpix ordering scheme. + * + * \note This does not clear the data in the map and it could be rendered invalid with this method. + */ + void Resize(int order, const std::valarray & specMin, const std::valarray & specMax, Healpix_Ordering_Scheme scheme=RING){ + Set(order,scheme); + fSpectra.resize(specMin.size()); + fSpecMin.resize(specMin.size()); + fSpecMax.resize(specMin.size()); + fSpectra=0.5*(specMin+specMax); + fSpecMin=specMin; + fSpecMax=specMax; + fBinned=true; + } + + /** \brief Return the number of bins in the spectra */ + int nSpectra() const { return fSpectra.size(); } + + /** \brief Return a constant reference to the values at which the spectra is evaluated */ + const std::valarray & getSpectra() const{ + return fSpectra; + } + /**\brief Return the boundaries for binned maps. + * + * \return false if the map is not binned + */ + bool getBoundaries(std::valarray & specMin, std::valarray & specMax) const{ + if (fBinned) { + specMin.resize(fSpecMin.size()); + specMax.resize(fSpecMax.size()); + specMin = fSpecMin; + specMax = fSpecMax; + } + return fBinned; + } + + /** \brief Convert to a Skymap */ + Skymap toSkymap() const{ + Skymap output; + if (fBinned){ + output.Resize(order_,fSpecMin,fSpecMax,scheme_); + }else{ + output.Resize(order_,fSpectra,scheme_); + } + fillFunctor f(output); + std::for_each(this->begin(), this->end(), f); + return output; + } + + /** \brief Add a healpix map to the sparse map + * + * \param map is the input Healpix map. All non-zero values + * are added to the current map. + * \param eBin is the bin number to add the map to + */ + void addHealpixMap(const Healpix_Map &map, int eBin) { + //Check the order + if (map.Order() != Order() ) { + std::cerr<<"Order does not match when adding Healpix map to SparseSkymap. Nothing added."<= fSpectra.size()){ + std::cerr<<"Energy bin out of bounds when adding Healpix map to SparseSkymap. Nothing added."<::iterator it = this->find(ind); + if ( it == this->end() ) { + it = this->insert(valueType(ind, T(0))).first; + } + it->second += map[i]; + } + } + } + + + /** \brief Convert from a Skymap + * + * \param map is the input skymap. All non-zero values are inserted into the sparse skymap. + */ + void fromSkymap(const Skymap & map){ + std::valarray specMin, specMax; + if (map.getBoundaries(specMin, specMax)){ + Resize(map.Order(),specMin,specMax,map.Scheme()); + }else{ + Resize(map.Order(), map.getSpectra(), map.Scheme()); + } + //Clear the map + this->clear(); + //Insert all values that are not zero; + for (long i = 0; i < map.Npix(); ++i){ + for (int j = 0; j < map.nSpectra(); ++j){ + if (map[i][j] != T(0)){ + //Check for existence of the pixel + indexType index(i,j); + (*this)[index] = map[i][j]; + } + } + } + } + + //Assignment operator + SparseSkymap & operator = (const SparseSkymap & oldMap){ + //Avoid self-assignment + if (this != &oldMap){ + if (oldMap.fBinned) { + Resize(oldMap.order_, oldMap.fSpecMin, oldMap.fSpecMax, oldMap.scheme_); + }else{ + Resize(oldMap.order_, oldMap.fSpectra, oldMap.scheme_); + } + this->clear(); + copyFunctor f(this); + std::for_each(oldMap.begin(), oldMap.end(), f); + } + return(*this); + } +}; + +#endif diff --git a/source/Spectrum.h b/source/Spectrum.h new file mode 100644 index 0000000000000000000000000000000000000000..3daa565d6602ce4f303da6b257547c70b63fc4c3 --- /dev/null +++ b/source/Spectrum.h @@ -0,0 +1,16 @@ + +//**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****| +// * Spectrum.h * galprop package * 4/14/2000 +//**"****!****"****!****"****!****"****!****"****!****"****!****"****!****"****| + +#ifndef Spectrum_h +#define Spectrum_h + +class Spectrum +{ + public: + + float *s; +}; + +#endif diff --git a/source/Sun.h b/source/Sun.h new file mode 100644 index 0000000000000000000000000000000000000000..49e5c17f6c6a712a3d42bd57afe2e515c3a5bdbb --- /dev/null +++ b/source/Sun.h @@ -0,0 +1,34 @@ +#ifndef Sun_h +#define Sun_h + +#include "Model.h" +#include "Exposure.h" +#include "Counts.h" +#include "Psf.h" +#include "Sources.h" +#include "Parameters.h" +#include "Variables.h" + + +class Sun : public BaseModel { + public: + /** \brief The constructor takes the same parameters as the BaseModel constructor. + * + * The constructor should load and prepare the data used in the model. + */ + + + + Sun(const Exposure &exposure, const CountsMap &counts, const Psf &psf, const Sources &sources, const Parameters &pars, const Skymap &filter, unsigned int configure = 3); //AWS20080717 + + + void getMap(const Variables &vars, Skymap &map); + gradMap getComponents(const Variables & vars, const std::string &prefix); + void getGrads(const Variables &vars, const std::string & varName, Skymap &map); + gradMap getSecDer(const Variables & vars, const std::string & varName1, const std::string & varName2); + private: + + Skymap fSun_disk_Map, fSun_IC_Map; +}; + +#endif diff --git a/source/UnConvolved.h b/source/UnConvolved.h new file mode 100644 index 0000000000000000000000000000000000000000..70b762b5c33b4ca462fabe22f61a13901909216c --- /dev/null +++ b/source/UnConvolved.h @@ -0,0 +1,74 @@ +//energy integrated, NOT convolved model skymaps + +#include"Distribution.h" +#include"HIH2IC.h" //AWS20081208 +#include"Sources.h" //AWS20081208 +#include"Skymap.h" //AWS20081208 + +class UnConvolved +{ + public: + + Distribution EGRET_bremss_skymap; + Distribution *EGRET_IC_iso_skymap; // 3 ISRF components + Distribution *EGRET_IC_aniso_skymap; // 3 ISRF components //AWS20060907 + Distribution EGRET_pi0_decay_skymap; + Distribution EGRET_isotropic_skymap; + Distribution EGRET_total_skymap; + + Distribution EGRET_HIR_skymap; //AWS20041223 + Distribution EGRET_H2R_skymap; //AWS20041223 + + Distribution EGRET_HIT_skymap; //AWS20050103 + Distribution EGRET_H2T_skymap; //AWS20050103 + + + + Distribution GLAST_bremss_skymap; //AWS20080514 + Distribution *GLAST_IC_iso_skymap; // 3 ISRF components + Distribution *GLAST_IC_aniso_skymap; // 3 ISRF components + Distribution GLAST_pi0_decay_skymap; + Distribution GLAST_isotropic_skymap; + Distribution GLAST_total_skymap; + + Distribution GLAST_HIR_skymap; + Distribution GLAST_H2R_skymap; + + Distribution GLAST_HIT_skymap; + Distribution GLAST_H2T_skymap; + +// GLAST models from gardian +// has to be a pointer apparently (otherwise galplot.cc gives compile error: "no matching function for Galplot()" ) +// something to do with how inheritance from BaseModel works ? +// JoinModel.cc also uses (a vector of) pointers to models +// HIH2IC GLAST_model_HIH2IC; //AWS20080718 NO + HIH2IC *GLAST_model_HIH2IC; //AWS20080718 YES + + + + Skymap GLAST_unconvolved_counts; //AWS20080718 + Skymap GLAST_unconvolved_counts_over_exposure;//AWS20080718 + + Skymap GLAST_unconvolved_counts_HI; //AWS20080718 + Skymap GLAST_unconvolved_counts_H2; //AWS20080718 + Skymap GLAST_unconvolved_counts_IC; //AWS20080718 + Skymap GLAST_unconvolved_counts_HII; //AWS20080810 + + Skymap GLAST_unconvolved_counts_sources; //AWS20081209 + Skymap GLAST_unconvolved_intensity_sources; //AWS20081212 + + + Skymap GLAST_unconvolved_counts_sourcepop_sublimit; //AWS20090109 + Skymap GLAST_unconvolved_counts_sourcepop_soplimit; //AWS20090109 + + + Skymap GLAST_unconvolved_intensity_sourcepop_sublimit; //AWS20090109 + Skymap GLAST_unconvolved_intensity_sourcepop_soplimit; //AWS20090109 + + Skymap GLAST_unconvolved_counts_solar_IC; //AWS20100215 + Skymap GLAST_unconvolved_intensity_solar_IC; //AWS20100215 + + Skymap GLAST_unconvolved_counts_solar_disk; //AWS20100215 + Skymap GLAST_unconvolved_intensity_solar_disk; //AWS20100215 + +}; diff --git a/source/Units.h b/source/Units.h new file mode 100644 index 0000000000000000000000000000000000000000..c023e7ecb7605187063fe413e9264c1e4a9a60ce --- /dev/null +++ b/source/Units.h @@ -0,0 +1,278 @@ +#ifndef _utl_SystemOfUnits_h_ +#define _utl_SystemOfUnits_h_ + +// Partition the global namespace in order to avoid ambiguity between +// Geant4 units + +namespace utl { + + /* + The conversion factors defined in this file + convert your data into a set of base units, so that + all dimensional quantities in the code are in a + single system of units. You can also + use the conversions defined here to, for example, + display data with the unit of your choice. For example: + \code + cout << "s = " << s/mm << " mm"; + \endcode + The base units are : + - meter (meter) + - nanosecond (nanosecond) + - electron Volt (eV) + - positron charge (eplus) + - degree Kelvin (kelvin) + - the amount of substance (mole) + - luminous intensity (candela) + - radian (radian) + - steradian (steradian) + + Below is a non-exhaustive list of derived and pratical units + (i.e. mostly the SI units). + + The SI numerical value of the positron charge is defined here, + as it is needed for conversion factor : positron charge = e_SI (coulomb) + + This is a slightly modified version of the units definitions + written by the Geant4 collaboration + + */ + + // + // Length [L] + // + static const double meter = 1.0; + static const double meter2 = meter*meter; + static const double meter3 = meter*meter*meter; + + static const double millimeter = 1.e-3*meter; + static const double millimeter2 = millimeter*millimeter; + static const double millimeter3 = millimeter*millimeter*millimeter; + + static const double centimeter = 10.*millimeter; + static const double centimeter2 = centimeter*centimeter; + static const double centimeter3 = centimeter*centimeter*centimeter; + + static const double kilometer = 1000.*meter; + static const double kilometer2 = kilometer*kilometer; + static const double kilometer3 = kilometer*kilometer*kilometer; + + static const double micrometer = 1.e-6*meter; + static const double micron = 1.e-6*meter; + static const double nanometer = 1.e-9*meter; + static const double angstrom = 1.e-10*meter; + static const double fermi = 1.e-15*meter; + + static const double barn = 1.e-28*meter2; + static const double millibarn = 1.e-3 *barn; + static const double microbarn = 1.e-6 *barn; + static const double nanobarn = 1.e-9 *barn; + static const double picobarn = 1.e-12*barn; + + // symbols + static const double mm = millimeter; + static const double mm2 = millimeter2; + static const double mm3 = millimeter3; + + static const double cm = centimeter; + static const double cm2 = centimeter2; + static const double cm3 = centimeter3; + + static const double m = meter; + static const double m2 = meter2; + static const double m3 = meter3; + + static const double km = kilometer; + static const double km2 = kilometer2; + static const double km3 = kilometer3; + + // + // Angle + // + static const double radian = 1.; + static const double milliradian = 1.e-3*radian; + static const double degree = (3.14159265358979323846/180.0)*radian; + + static const double steradian = 1.; + + // symbols + static const double rad = radian; + static const double mrad = milliradian; + static const double sr = steradian; + static const double deg = degree; + + // + // Time [T] + // + static const double nanosecond = 1.; + static const double second = 1.e+9 *nanosecond; + static const double millisecond = 1.e-3 *second; + static const double microsecond = 1.e-6 *second; + static const double picosecond = 1.e-12*second; + static const double minute = 60*second; + static const double hour = 60*minute; + static const double day = 24*hour; + + static const double hertz = 1./second; + static const double kilohertz = 1.e+3*hertz; + static const double megahertz = 1.e+6*hertz; + + // symbols + static const double ns = nanosecond; + static const double s = second; + static const double ms = millisecond; + + // + // Electric charge [Q] + // + static const double eplus = 1. ; // positron charge + static const double e_SI = 1.602176462e-19; // positron charge in coulomb + static const double coulomb = eplus/e_SI; // coulomb = 6.24150 e+18*eplus + + // + // Energy [E] + // + static const double electronvolt = 1.; + static const double megaelectronvolt = 1.e+6*electronvolt; + static const double kiloelectronvolt = 1.e+3*electronvolt; + static const double gigaelectronvolt = 1.e+9*electronvolt; + static const double teraelectronvolt = 1.e+12*electronvolt; + static const double petaelectronvolt = 1.e+15*electronvolt; + static const double exaelectronvolt = 1.e+18*electronvolt; + static const double zettaelectronvolt= 1.e+21*electronvolt; + + static const double joule = electronvolt/e_SI; // joule = 6.24150 e+12 * MeV + static const double erg = 1.0e-7*joule; + + // symbols + static const double MeV = megaelectronvolt; + static const double eV = electronvolt; + static const double keV = kiloelectronvolt; + static const double GeV = gigaelectronvolt; + static const double TeV = teraelectronvolt; + static const double PeV = petaelectronvolt; + static const double EeV = exaelectronvolt; + static const double ZeV = zettaelectronvolt; + + // + // Mass [E][T^2][L^-2] + // + static const double kilogram = joule*second*second/(meter*meter); + static const double gram = 1.e-3*kilogram; + static const double milligram = 1.e-3*gram; + + // symbols + static const double kg = kilogram; + static const double g = gram; + static const double mg = milligram; + + // + // Power [E][T^-1] + // + static const double watt = joule/second; // watt = 6.24150 e+3 * MeV/ns + + // + // Force [E][L^-1] + // + static const double newton = joule/meter; // newton = 6.24150 e+9 * MeV/mm + + // + // Pressure [E][L^-3] + // + static const double hep_pascal = newton/m2; // pascal = 6.24150 e+3 * MeV/mm3 + static const double bar = 100000*hep_pascal; // bar = 6.24150 e+8 * MeV/mm3 + static const double atmosphere = 101325*hep_pascal; // atm = 6.32420e+8*MeV/mm3 + + // + // Electric current [Q][T^-1] + // + static const double ampere = coulomb/second; // ampere = 6.24150e+9*eplus/ns + static const double milliampere = 1.e-3*ampere; + static const double microampere = 1.e-6*ampere; + static const double nanoampere = 1.e-9*ampere; + + // + // Electric potential [E][Q^-1] + // + static const double megavolt = megaelectronvolt/eplus; + static const double kilovolt = 1.e-3*megavolt; + static const double volt = 1.e-6*megavolt; + + // + // Electric resistance [E][T][Q^-2] + // + static const double ohm = volt/ampere; // ohm = 1.60217e-16*(MeV/eplus)/(eplus/ns) + + // + // Electric capacitance [Q^2][E^-1] + // + static const double farad = coulomb/volt; // farad = 6.24150e+24 * eplus/Megavolt + static const double millifarad = 1.e-3*farad; + static const double microfarad = 1.e-6*farad; + static const double nanofarad = 1.e-9*farad; + static const double picofarad = 1.e-12*farad; + + // + // Magnetic Flux [T][E][Q^-1] + // + static const double weber = volt*second; // weber = 1000*megavolt*ns + + // + // Magnetic Field [T][E][Q^-1][L^-2] + // + static const double tesla = volt*second/meter2; // tesla =0.001*megavolt*ns/mm2 + + static const double gauss = 1.e-4*tesla; + static const double kilogauss = 1.e-1*tesla; + + // + // Inductance [T^2][E][Q^-2] + // + static const double henry = weber/ampere; // henry = 1.60217e-7*MeV*(ns/eplus)**2 + + // + // Temperature + // + static const double kelvin = 1.; + + // + // Amount of substance + // + static const double mole = 1.; + + // + // Activity [T^-1] + // + static const double becquerel = 1./second ; + static const double curie = 3.7e+10 * becquerel; + + // + // Absorbed dose [L^2][T^-2] + // + static const double gray = joule/kilogram ; + + // + // Luminous intensity [I] + // + static const double candela = 1.; + + // + // Luminous flux [I] + // + static const double lumen = candela*steradian; + + // + // Illuminance [I][L^-2] + // + static const double lux = lumen/meter2; + + // + // Miscellaneous + // + static const double perCent = 0.01 ; + static const double perThousand = 0.001; + static const double perMillion = 0.000001; + +} // end of utl namespace + +#endif diff --git a/source/Utils.h b/source/Utils.h new file mode 100644 index 0000000000000000000000000000000000000000..cd7c9ef5b1586708bb8bc2addb73579adf8498a9 --- /dev/null +++ b/source/Utils.h @@ -0,0 +1,79 @@ +#ifndef Utils_h +#define Utils_h + +#include +#include + +/** \brief Several utilities common to many classes are stored in this namespace */ +namespace Utils{ + /** \brief Assuming array is sorted ascendantly, find the index where the array value is just below the searched for value. + * + * \param array the ascendantly sorted array + * \param value the value to search for + * + * \note Returns the lowest index in the array if value is lower than all values in the array + */ + int indexBelow(const std::valarray &array, double value); + /** \brief Assuming array is sorted ascendantly, find the index where the array value is just above the searched for value. + * + * \param array the ascendantly sorted array + * \param value the value to search for + * + * \note Returns the highest index in the array if value is higher than all values in the array + */ + int indexAbove(const std::valarray &array, double value); + /** \brief Find the index range which brackets the range from lowValue to highValue + * + * \param array the ascendantly sorted array to search + * \param lowValue the lower limit of the range to bracket + * \param highValue the upper limit of the range to bracket + * \param lowIndex a reference to a integer, which will contain the lower index in the array + * \param highIndex a reference to a integer, which will contain the higher index in the array + * + * If the range is outside the array span, the index will be in the end closer to the range. The + * index will only be the same if the array size is less than 2. + */ + void findIndexRange(const std::valarray &array, double lowValue, double highValue, int &lowIndex, int &highIndex); + /** \overload */ + void findIndexRange(const std::vector &array, double lowValue, double highValue, int &lowIndex, int &highIndex); + /** Linear interpolation of values + * + * \param x the value at which to evaluate the interpolation + * \param xValues the vector of x values in the data to interpolate + * \param yValues the vector of y values in the data to interpolate + * + * This method extrapolates silently. + */ + double linearInterpolation(double x, const std::vector &xValues, const std::vector &yValues); + /** \brief Calculate the 3rd order Lagrange Coefficients, given x and y values. + * + * \param x the x values to use + * \param y the y values to use + * \param coeff on return, stores the calculated coefficients. The coefficients are stored in order, so the total + * function will be \f$ f = \sum^i coeff_i x^i \f$ + */ + void lagrangeCoeff4(const std::valarray & x,const std::valarray & y, std::valarray & coeff); + + /** \brief Nicely formed status indicator, showing percent completed in both numbers and as a bar. + * + * This class is OMP thread safe. + */ + class StatusIndicator { + public: + /** \brief Construct a status indicator with a string and number of steps + * + * \param marker is the string to print in front of the status indicator + * \param steps the total number of steps expected + */ + StatusIndicator(const std::string & marker, int steps); + /** \brief Increase the number of finished steps by one and reprint. */ + void refresh(); + private: + std::string fmarker; //!< Store the marker + int fsteps, fcurrentStep; //!< Number of steps and current step. + long ftInit; //!< Initial time when counter started, used for time estimates + double ftEstimate; //!< Estimated time needed to complete the project + }; +} + +#endif diff --git a/source/Variables.h b/source/Variables.h new file mode 100644 index 0000000000000000000000000000000000000000..fcc9c6857633e7b74397e7a27cf2bef427a3bfa2 --- /dev/null +++ b/source/Variables.h @@ -0,0 +1,171 @@ +#ifndef Variables_h +#define Variables_h + +#include "Parameters.h" +#include +#include +#include +#include +#include + +/** \brief Handles communication of free variables between models and fitters. + * + * Stores the names, value and limits of variables. + */ +class Variables { + public: + //! Default constructor to initialize the variables + Variables() : fmaxLength(0) {} + /** \brief Add variables from a Parameters object. + * + * \param names is the variable names to look up in the Parameters object + * \param pars is the Parameters object containing the parameters and their value. + * The format for the variables in the parameter object has to be + * \code + * name = value errorEstimate + * \endcode + * The limits are optional but both have to be present if one is. + */ + void add(const std::vector & names, const Parameters &pars); + /** \brief Add variables from another Variables object */ + void add(const Variables & vars); + /** \brief Add a variable without limits + * + * \param name is the variable name + * \param value is its value + * \param error is its error boundary + */ + void add(const std::string &name, double value, double error); + /** \brief Add a variable with limits + * + * \param name is the variable name + * \param value is its value + * \param error is its error boundary + * \param lowerLimit is its lower limit + * \param upperLimit is its upper limit + */ + void add(const std::string &name, double value, double error, double lowerLimit, double upperLimit); + + /** \brief Return a reference to the variable value */ + double & operator [] (const std::string &name); + /** \brief Return a constant reference to the variable value */ + const double & operator [] (const std::string &name) const; + /** \brief Return a reference to the variable value */ + double & value(const std::string &name); + /** \brief Return a constant reference to the variable value */ + const double & value(const std::string &name) const; + + /** \brief Return a reference to the variable error */ + double & error(const std::string &name); + /** \brief Return a constant reference to the variable error */ + const double & error(const std::string &name) const; + + /** \brief Return a vector with all the variable names */ + std::vector getNames() const; + + /** \brief Set the upper limit of a variable + * \param name is the variable name + * \param limit is the new upper limit. + */ + void setUpperLimit(const std::string &name, double limit); + /** \brief Set the lower limit of a variable + * \param name is the variable name + * \param limit is the new lower limit. + */ + void setLowerLimit(const std::string &name, double limit); + + /** \brief Return a constant reference to the variable upper limit, if it exists */ + const double & upperLimit(const std::string &name) const; + /** \brief Return a constant reference to the variable lower limit, if it exists */ + const double & lowerLimit(const std::string &name) const; + + /** \brief Return the values of the variables as a vector. + * + * Uses the alphabetical order of the parameters + */ + std::vector getValueVector() const; + /** \brief Sets the values of the variables from a vector. + * + * Uses the alphabetical order of the parameters and throws an error if the sizes do not match. + */ + void setValueVector(const std::vector & values); + + //! Clear all set variables + void clear(); + + /** \brief Default error class */ + class VariableError : public std::exception { + private: + std::string eMessage; + public: + /** \brief Sets a default error string "Error in Variables class" */ + VariableError() : eMessage("Error in Variables class"){} + /** \brief Sets a specialized error string */ + VariableError(const std::string & message) : eMessage(message){} + ~VariableError() throw(){} + /** \brief Returns the error string */ + virtual const char* what () const throw(){ + return(eMessage.c_str()); + } + }; + /** \brief Thrown when the variable name does not exist */ + class VariableNotFound : public std::exception { + private: + std::string eMessage; + public: + /** \brief Sets a default error string "Variable not found" */ + VariableNotFound() : eMessage("Variable not found"){} + /** \brief Sets a specialized error string */ + VariableNotFound(const std::string & message) : eMessage(message){} + ~VariableNotFound() throw(){} + /** \brief Returns the error string */ + virtual const char* what () const throw(){ + return(eMessage.c_str()); + } + }; + /** \brief Thrown when trying to access a limit that is not there */ + class LimitNotSet : public std::exception { + private: + std::string eMessage; + public: + /** \brief Sets a default error string "Limit not set" */ + LimitNotSet() : eMessage("Limit not set"){} + /** \brief Sets a specialized error string */ + LimitNotSet(const std::string & message) : eMessage(message){} + ~LimitNotSet() throw(){} + /** \brief Returns the error string */ + virtual const char* what () const throw(){ + return(eMessage.c_str()); + } + }; + /** \brief Thrown when trying to create a variable with a name that already exists */ + class VariableAlreadyCreated : public std::exception { + private: + std::string eMessage; + public: + /** \brief Sets a default error string "Variable already created" */ + VariableAlreadyCreated() : eMessage("Variable already created"){} + /** \brief Sets a specialized error string */ + VariableAlreadyCreated(const std::string & message) : eMessage(message){} + ~VariableAlreadyCreated() throw(){} + /** \brief Returns the error string */ + virtual const char* what () const throw(){ + return(eMessage.c_str()); + } + }; + /** \brief Output operator for the Variables class */ + friend std::ostream & operator << (std::ostream &os, const Variables &vars); + + private: + /** \brief A structure to store each variable */ + struct variable { + double value, error, uLimit, lLimit; + bool uSet, lSet; + }; + std::map fvars; + typedef std::map::iterator varIt; + typedef std::map::const_iterator cvarIt; + int fmaxLength; //!< The maximum lenght of a variable name +}; + +#endif diff --git a/source/constants.h b/source/constants.h new file mode 100644 index 0000000000000000000000000000000000000000..9d2063511c4b70811c82f186b9a99a9a45b32c3e --- /dev/null +++ b/source/constants.h @@ -0,0 +1,72 @@ + +//**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****| +// * constants.h * galprop package * 08/16/2001 +//**"****!****"****!****"****!****"****!****"****!****"****!****"****!****"****| + +#ifndef constants_h +#define constants_h + +// aws routines +const double Rsun = 8.5; // kpc - Galactocentric radius of the Sun AWS20100118 for consistency with galprop and used by nH.cc + +const double year2sec=365.25*24.*3600.; +const double kpc2cm =3.08568e21; +const double c =3.0e10; + +const double m_electron = 0.511; // MeV +const double m_proton = 938. ; // MeV + +const double h_planck=6.6260755e-27 ; // erg sec +const double eV_to_erg=1.60217733e-12; +const double erg_to_eV=1./eV_to_erg; + +// imos routines + +const double // PHYSICAL CONSTANTS: + Pi = 3.141592653589793, // Pi number + C = 2.99792458e10, // cm/s, =c speed of light + H = 6.626075540e-27, // erg*s, =h Planck constant + H2Pi= 6.582122020e-22, // MeV*s, =h/(2Pi) Planck constant + ALPHAf=1./137.035989561, // =e^2/(hc/2Pi) fine-structure const. + Rele= 2.8179409238e-13, // cm, =e^2/mc^2 class. electron radius + MEV2ERG = 1.6021773349e-6, // MeV/erg, conversion const. + H2PiC = 0.19732705359e-10,// MeV*cm, =hc/(2Pi) conversion const. + // IONIZATION POTENTIALS: + EiH = 13.6e-6, // MeV, Hydrogen (in electron energy losses) + EiHe = 24.59e-6, // MeV, Helium (in electron energy losses) + EH = 19.e-6, // MeV, H eff. ioniz. potential (in nucl. loss) + EHe= 44.e-6, // MeV, He eff. ioniz. potential (in nucl. loss) + // RADIATIVE LENGHTS: + TH = 62.8, // g/cm^2, Hydrogen + THe= 93.1, // g/cm^2, Helium + // MASSES: Particle Data Group (1990, Phys.Lett. 239) + amu = 0.931494, // GeV, mass 12C /12= atomic mass unit // IMOS20010816 + Mele = 0.5109990615, // MeV/c^2, electron rest mass + MD = 1.232, // GeV, Delta(1232)-isobar mean mass + GD = 0.0575, // GeV, Delta-isobar wight (0.115/2.) + Mp = 0.938, // GeV, Proton rest mass + Mn = 0.9396, // GeV, Neutron rest mass + Md = 1.8756, // GeV, Deutron rest mass + Mpi0 = 0.135, // GeV, Neutral pion rest mass + Mpi1 = 0.1396, // GeV, Charged pion rest mass + Mmu = 0.10566, // GeV, Muon rest mass + MK = 0.49365, // GeV, Kaon rest mass + // REACTION THRESHOLDS: (from kinematics) + Pth0 = 0.78, // GeV/c, min momentum for pp->pi^o (+2p) + Pth1 = 1.65, //(1.22)GeV/c, -"- for pp->pi^-(+2p+pi^+) + Pth2 = 0.80, // GeV/c, -"- for pp->pi^+ (+p+n) + Pth3 = 0.791, // GeV/c, -"- for pp->pi^+ (+d) + Pth4 = 3.302, // GeV/c, -"- for pp->K^- (+2p+K^+) + Pth5 = 1.8332, // GeV/c, -"- for pp->K^+ (+p+n) + // BRANCHING RATIOS: + BR1 = 0.635, // br. ratio for channel K -> mu +nu + BR2 = 0.212, // br. ratio for channel K -> pi^(+/-) +pi^0 + // Helium FORM FACTORS: [Gould 1969, Phys.Rev.185,72 (p.77)] + DD[11] = {0., 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1., 2., 5., 10.}, + PH1[11]= {134.60, 133.85, 133.11, 130.86, 127.17, + 120.35, 104.60, 89.94, 74.19, 54.26, 40.94}, + PH2[11]= {131.40, 130.51, 130.33, 129.26, 126.76, + 120.80, 105.21, 89.46, 73.03, 51.84, 37.24}; + + +#endif diff --git a/source/fort_interface.h b/source/fort_interface.h new file mode 100644 index 0000000000000000000000000000000000000000..a54d4179f1ad12d5c3239739e15172198ce681e9 --- /dev/null +++ b/source/fort_interface.h @@ -0,0 +1,33 @@ + +//**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****| +// * fort_interface.h * galprop package * 2001/05/11 +//**"****!****"****!****"****!****"****!****"****!****"****!****"****!****"****| + +#ifndef fort_interface_h +#define fort_interface_h + +// fort_interface1.cc prototypes + +void set_sigma_cc(); +double e_loss_compton_cc(double,double); +double bremss_spec_cc(double,double,int,int); +double pp_meson_cc(double,double,int,int,int); +int test_e_loss_compton_cc(); +int test_bremss_spec_cc(); +int test_pp_meson_cc(); + +// fort_interface2.cc prototypes + +double wsigma_cc(int,int,int,int,double); // IMOS20020502 +double yieldx_cc(int,int,int,int,float); // IMOS20020502 +void sigtap_cc(int); // IMOS20010511 +double sighad_cc(int,double,double,double,double,double); // IMOS20020502 +double cfactor_cc(int,double,double,double,double,double,double,double,double); +double antiproton_cc(int,double,double,int,int,int,int); // IMOS20010511 +double synchrotron_cc(double,double,double); +int test_antiproton_cc(); +int test_synchrotron_cc(); +int test_cfactor(); + +#endif + diff --git a/source/galplot.h b/source/galplot.h new file mode 100644 index 0000000000000000000000000000000000000000..d7e311f606304e7361fdbdfbdb5acc0a8a560f41 --- /dev/null +++ b/source/galplot.h @@ -0,0 +1,106 @@ + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include //AWS20041021 +#include + +#include "Galplotdef.h" +#include "Data.h" +#include "Convolved.h" +#include "UnConvolved.h" +#include "GCR_data.h" +#include "FITS.h" +#include "Skymap.h" //AWS20140721 for energy dispersion routine + +int plot_HIR(); +int read_IC_skymap(char *IC_type); + +int read_bremss_skymap(); +int read_bremss_HIR_skymap(); //AWS20041223 +int read_bremss_H2R_skymap(); //AWS20041223 + +int read_pi0_decay_skymap(); +int read_pi0_decay_HIR_skymap();//AWS20041223 +int read_pi0_decay_H2R_skymap();//AWS20041223 + +int read_synchrotron_skymap() ;//AWS20050726 +int plot_synchrotron_skymap() ;//AWS20050726 +int read_synchrotron_data () ;//AWS20050728 + +int plot_skymaps(); +int plot_profile(int longprof, int latprof, int ip,int ic,int bremss,int pi0,int total); +int plot_convolved_EGRET_profile(int longprof, int latprof,int ip,int ic,int bremss,int pi0,int total); + +int read_EGRET_data(); +int read_EGRET_psf(); +int plot_EGRET_skymaps(int mode); +int plot_EGRET_profile(int longprof, int latprof,int ip,int mode); +int plot_EGRET_spectrum(int mode); + +int read_COMPTEL_data(); +int plot_COMPTEL_spectrum(int mode); + +int plot_SPI_spectrum (int mode);// AWS20040405 +int plot_SPI_spectrum_spimodfit (int mode);// AWS20050308 +int plot_RXTE_spectrum(int mode);// AWS20040407 +int plot_OSSE_spectrum(int mode);// AWS20041001 +int plot_IBIS_spectrum(int mode);// AWS20041015 +int plot_GINGA_spectrum(int mode);// AWS20041102 +int plot_model_ridge (int mode);// AWS20041102 +int plot_MILAGRO_spectrum (int mode);//AWS20050518 +int plot_Chandra_spectrum(int mode); //AWS20050818 + +int plot_spectrum(int ic,int bremss,int pi0,int total); +int plot_isotropic_sreekumar_spectrum(); +int plot_isotropic_EGRET_spectrum(); +int fit_EGRET(int options); +int fit_EGRET_Xco(int options);//AWS20050104 +int plot_hunter_10_50GeV_spectrum(); + +int plot_gcr_spectra(); +int plot_gcr_spectra_data(int Z,int A,int index_mult);//AWS20091009 +int plot_gcr_spectra_ratios(); +int plot_gcr_spectra_ratios_data(int i_ratio); +int convolve_EGRET(); +int convolve_EGRET_HI_H2();//AWS20041223 +int convolve(Distribution &map, Distribution psf); +double sabin(double Y0,double XBINSZ,double YBINSZ,int IY); + +double energy_integral +(double emin,double emax,double e1,double e2,double intensity_esq1,double intensity_esq2); +int test_energy_integral (); + +void modulate(double *Ekin, float *cr_density, int np, int z, int a,double phi, int key); + + + +int plot_SPI_spiskymax_profile(int longprof,int latprof,int mode);//AWS20050720 +int plot_SPI_model_profile(int longprof, int latprof,double e_min,double e_max,int ic,int bremss,int total);//AWS20050720 +int convolve_SPI(double emin,double emax);//AWS20050721 + + +int test_Skymap(int verbose); //AWS20080516 +int EBV_convert(int verbose); //AWS20090820 + + +int Skymap_Energy_Dispersion(Skymap &skymap, double E_interp_factor, int interpolation, string parameter_file_name, string parameter_file_type, string exposure_file_name, string exposure_file_type, int use_Aeff, double E_true_threshold, double E_meas_threshold, int debug); //AWS20150301 + +int spectrum_Energy_Dispersion(valarray E, valarray &spectrum, double E_interp_factor, int interpolation, string parameter_file_name, string parameter_file_type, string exposure_file_name, string exposure_file_type, int use_Aeff, double E_true_threshold, double E_meas_threshold, int debug); //AWS20150301 diff --git a/source/galprop.h b/source/galprop.h new file mode 100644 index 0000000000000000000000000000000000000000..393b7125369c770f4b62eaa2d764c92f278fe537 --- /dev/null +++ b/source/galprop.h @@ -0,0 +1,176 @@ + +//**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****| +// * galprop.h * galprop package * 08/16/2001 +//**"****!****"****!****"****!****"****!****"****!****"****!****"****!****"****| + +// function prototypes + +#ifndef galprop_h +#define galprop_h +using namespace std; //AWS20050912 +#include //AWS20050912 +#include //AWS20051124 +#include + +#include "fort_interface.h" +#include "constants.h" + +void nucleon_cs(int,double,int,int,int,double*,double*,double*,double*,double*,double*);// IMOS20010511 +void read_nucdata(); +void Kcapture_cs(double,int,int,double*,double*); // IMOS20010816 +double isotope_cs(double,int,int,int,int,int,int*); +double nucdata(int,int,int,int,int,int,int*,int*,double*); // IMOS20010816 +double nucleon_loss(int,int,double,double,double,double, // nucleon + double*, double*); // energy losses +double electron_loss(double,double,double,double,double,double, // electron + double*,double*,double*,double*,double*,double*); // energy losses + +int create_gcr(); +int create_transport_arrays(Particle&); +int create_galaxy(); +int create_SNR(); +int source_SNR_event(Particle&,double); +void decayed_cross_sections(int,int,int,int,double*,int,double*); // IMOS20010816 +double sigma_boron_dec_heinbach_simon(int,int,int,int,double); +int gen_secondary_source(Particle&); +void read_nucdata(); + +int kinematic(int,int,char*,double&,double&,double&,double&,double&,double&,int); +int print_BC(); +int propagate_particles(); +int propel(Particle &particle); +int tridag (float*,float*,float*,float*,float*,int); +int tridag_sym(float*,float*,float*,float*,float*,int); + +void protri(Particle &particle, + +Distribution &alpha1_x,Distribution &alpha1_y,Distribution &alpha1_z,Distribution &alpha1_p, +Distribution &alpha2_x,Distribution &alpha2_y,Distribution &alpha2_z,Distribution &alpha2_p, +Distribution &alpha3_x,Distribution &alpha3_y,Distribution &alpha3_z,Distribution &alpha3_p, + +Distribution &Nx1_, Distribution &Ny1_, Distribution &Nz1_, Distribution &Np1_, +Distribution &Nx2_, Distribution &Ny2_, Distribution &Nz2_, Distribution &Np2_, +Distribution &Nx3_, Distribution &Ny3_, Distribution &Nz3_, Distribution &Np3_, + +Distribution &total_source_function, + +double dt, int nrept_outer ,double f_use); + +int source_SNR_event_vec(Particle &particle,double t, + float *total_source_function_x); + +int store_gcr(); +int store_gcr_full(); +int read_gcr(); + +int read_isrf(); +int read_HIR(); +int read_COR(); + + +int gen_isrf_energy_density(); +int He_to_H_CS(double,int,int,int,int,double*,double*); +double D_pp(double,double,double,double); + +int e_KN_loss(Particle&); +double IC_cross_section(double,double,double); +double IC_anisotropy_factor(double,double,double,int); +double ionization_bethe(int Z, double beta); + +double nHI (double,double); +double nH2 (double,double); +double nHII(double,double); + +double nHI_av (double,double,double,double,double); +double nH2_av (double,double,double,double,double); +double nHII_av(double,double,double,double,double); + + +double source_distribution (double,double,double); +double B_field_model(double,double,int); +double B_field_model(double,double,double,int); + +int propel_diagnostics + (Particle &particle, + Distribution &alpha1_r, + Distribution &alpha1_z, + Distribution &alpha1_p, + Distribution &alpha2_r, + Distribution &alpha2_z, + Distribution &alpha2_p, + Distribution &alpha3_r, + Distribution &alpha3_z, + Distribution &alpha3_p, + Distribution &total_source_function, + double dt + ); + +int propel_diagnostics + (Particle &particle, + Distribution &alpha1_x, + Distribution &alpha1_y, + Distribution &alpha1_z, + Distribution &alpha1_p, + Distribution &alpha2_x, + Distribution &alpha2_y, + Distribution &alpha2_z, + Distribution &alpha2_p, + Distribution &alpha3_x, + Distribution &alpha3_y, + Distribution &alpha3_z, + Distribution &alpha3_p, + Distribution &total_source_function, + double dt + ); + +int electrons_normalize(); +int nuclei_normalize(); +int gen_IC_emiss(); +int gen_IC_skymap(); +int store_IC_skymap(char *IC_type); +int store_IC_skymap_comp(char *IC_type); +int gen_bremss_emiss(); +int store_bremss_emiss(); +int store_bremss_ionized_skymap(); +int gen_bremss_ionized_skymap(); +int gen_bremss_skymap(); +int store_bremss_skymap(); +int gen_pi0_decay_emiss(); +int store_pi0_decay_emiss(); +int gen_pi0_decay_skymap(); +int store_pi0_decay_skymap(); +int HIR_iRing(double RR); + +int gen_synch_emiss(); +int gen_synch_skymap(); +int store_synch_skymap(); +int store_ionization_rate(); + +int gen_secondary_positron_source (Particle &particle); +int gen_tertiary_antiproton_source (Particle &particle); // IMOS20000605.15 +int gen_secondary_antiproton_source(Particle &particle); +int gen_secondary_proton_source (Particle &particle); // IMOS20000605.16 + +int gen_ionization_rate(); +int cr_luminosity(); + +float isrf_energy_density(float rr, float zz); + +void skeleton(); +void skeleton2(); + +double gauss(double mean, double sigma); + +int test_sigma_boron_dec_heinbach_simon(); +int test_kinematic(); +int test_He_to_H_CS(); +int test_nH(); +int test_Distribution(); +int test_Particle(); +int test_isotope_cs(); +int test_float_accuracy(); +int test_source_SNR_event(); +int test_suite(); +int test_gauss(); //AWS20010410 + +#endif diff --git a/source/galprop_classes.h b/source/galprop_classes.h new file mode 100644 index 0000000000000000000000000000000000000000..9e0f3ec929248299e1e398a63f39babe13b0ea50 --- /dev/null +++ b/source/galprop_classes.h @@ -0,0 +1,18 @@ + +//**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****| +// * galprop_classes.h * galprop package * 4/14/2000 +//**"****!****"****!****"****!****"****!****"****!****"****!****"****!****"****| + +#ifndef galprop_classes_h +#define galprop_classes_h + +#include "Configure.h" +#include "Galdef.h" +#include "Spectrum.h" +#include "Distribution.h" +#include "Particle.h" +#include "Galaxy.h" + +#include "galplot.h" + +#endif