source/*h

source/ArraySlice.h

+#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 <iostream>
+#include <stdexcept>
+#include <valarray>
+
+template<typename T>
+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<T> &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 <typename C>
+	 ArraySlice<T> & 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<T> & operator = (const ArraySlice<T> & slice) {
+	 if (fsize != slice.fsize) {
+	    std::string errMsg = "Slices need to have same size in assignment";
+	    std::cerr<<errMsg<<std::endl;
+	    throw(std::invalid_argument(errMsg));
+	 }
+	 for (size_t i = 0; i < fsize; ++i) {
+	    *(fbegin+i) = slice[i];
+	 }
+	 return (*this);
+      }
+      template <typename C>
+	 ArraySlice<T> & operator = (const ArraySlice<C> & slice) {
+	    if (fsize != slice.fsize) {
+	       std::string errMsg = "Slices need to have same size in assignment";
+	       std::cerr<<errMsg<<std::endl;
+	       throw(std::invalid_argument(errMsg));
+	    }
+	    for (size_t i = 0; i < fsize; ++i) {
+	       *(fbegin+i) = slice[i];
+	    }
+	    return (*this);
+	 }
+      template <typename C>
+	 ArraySlice<T> & operator = (const std::valarray<C> & valarray) {
+	    if (fsize != valarray.size()) {
+	       std::string errMsg = "Valarray and slice need to have same size in assignment";
+	       std::cerr<<errMsg<<std::endl;
+	       throw(std::invalid_argument(errMsg));
+	    }
+	    for (size_t i = 0; i < fsize; ++i){
+	       *(fbegin+i) = valarray[i];
+	    }
+	    return (*this);
+	 }
+
+
+      //Mathematical operators should be self-explanatory
+      //Includes operation with valarrays and regular numbers
+      
+      //Addition with a single number
+      template <typename C>
+	 ArraySlice<T> & operator += (C number) {
+	    for (size_t i = 0; i < fsize; ++i){
+	       *(fbegin+i) += number;
+	    }
+	    return (*this);
+	 }
+      //Subtraction with a single number
+      template <typename C>
+	 ArraySlice<T> & operator -= (C number) {
+	    for (size_t i = 0; i < fsize; ++i){
+	       *(fbegin+i) -= number;
+	    }
+	    return (*this);
+	 }
+      //Multiplication with a single number
+      template <typename C>
+	 ArraySlice<T> & operator *= (C number) {
+	    for (size_t i = 0; i < fsize; ++i){
+	       *(fbegin+i) *= number;
+	    }
+	    return (*this);
+	 }
+      //Division with a single number
+      template <typename C>
+	 ArraySlice<T> & operator /= (C number) {
+	    for (size_t i = 0; i < fsize; ++i){
+	       *(fbegin+i) /= number;
+	    }
+	    return (*this);
+	 }
+
+      //Addition with another slice
+      template <typename C>
+	 ArraySlice<T> & operator += (const ArraySlice<C> & slice) {
+	    if (fsize != slice.fsize) {
+	       std::string errMsg = "Slices need to have same size in addition";
+	       std::cerr<<errMsg<<std::endl;
+	       throw(std::invalid_argument(errMsg));
+	    }
+	    for (size_t i = 0; i < fsize; ++i){
+	       *(fbegin+i) += slice[i];
+	    }
+	    return (*this);
+	 }
+      //Subtraction with another slice
+      template <typename C>
+	 ArraySlice<T> & operator -= (const ArraySlice<C> & slice) {
+	    if (fsize != slice.fsize) {
+	       std::string errMsg = "Slices need to have same size in subtraction";
+	       std::cerr<<errMsg<<std::endl;
+	       throw(std::invalid_argument(errMsg));
+	    }
+	    for (size_t i = 0; i < fsize; ++i){
+	       *(fbegin+i) -= slice[i];
+	    }
+	    return (*this);
+	 }
+      //Multiplication with another slice
+      template <typename C>
+	 ArraySlice<T> & operator *= (const ArraySlice<C> & slice) {
+	    if (fsize != slice.fsize) {
+	       std::string errMsg = "Slices need to have same size in multiplication";
+	       std::cerr<<errMsg<<std::endl;
+	       throw(std::invalid_argument(errMsg));
+	    }
+	    for (size_t i = 0; i < fsize; ++i){
+	       *(fbegin+i) *= slice[i];
+	    }
+	    return (*this);
+	 }
+      //Division with another slice
+      template <typename C>
+	 ArraySlice<T> & operator /= (const ArraySlice<C> & slice) {
+	    if (fsize != slice.fsize) {
+	       std::string errMsg = "Slices need to have same size in division";
+	       std::cerr<<errMsg<<std::endl;
+	       throw(std::invalid_argument(errMsg));
+	    }
+	    for (size_t i = 0; i < fsize; ++i){
+	       *(fbegin+i) /= slice[i];
+	    }
+	    return (*this);
+	 }
+
+      //Addition with another valarray
+      template <typename C>
+	 ArraySlice<T> & operator += (const std::valarray<C> & valarray) {
+	    if (fsize != valarray.size()) {
+	       std::string errMsg = "Valarray and slice need to have same size in addition";
+	       std::cerr<<errMsg<<std::endl;
+	       throw(std::invalid_argument(errMsg));
+	    }
+	    for (size_t i = 0; i < fsize; ++i){
+	       *(fbegin+i) += valarray[i];
+	    }
+	    return (*this);
+	 }
+      //Subtraction with another valarray
+      template <typename C>
+	 ArraySlice<T> & operator -= (const std::valarray<C> & valarray) {
+	    if (fsize != valarray.size()) {
+	       std::string errMsg = "Valarray and slice need to have same size in subtraction";
+	       std::cerr<<errMsg<<std::endl;
+	       throw(std::invalid_argument(errMsg));
+	    }
+	    for (size_t i = 0; i < fsize; ++i){
+	       *(fbegin+i) -= valarray[i];
+	    }
+	    return (*this);
+	 }
+      //Multiplication with another valarray
+      template <typename C>
+	 ArraySlice<T> & operator *= (const std::valarray<C> & valarray) {
+	    if (fsize != valarray.size()) {
+	       std::string errMsg = "Valarry and slice need to have same size in multiplication";
+	       std::cerr<<errMsg<<std::endl;
+	       throw(std::invalid_argument(errMsg));
+	    }
+	    for (size_t i = 0; i < fsize; ++i){
+	       *(fbegin+i) *= valarray[i];
+	    }
+	    return (*this);
+	 }
+      //Division with another valarray
+      template <typename C>
+	 ArraySlice<T> & operator /= (const std::valarray<C> & valarray) {
+	    if (fsize != valarray.size()) {
+	       std::string errMsg = "Valarray and slice need to have same size in division";
+	       std::cerr<<errMsg<<std::endl;
+	       throw(std::invalid_argument(errMsg));
+	    }
+	    for (size_t i = 0; i < fsize; ++i){
+	       *(fbegin+i) /= valarray[i];
+	    }
+	    return (*this);
+	 }
+
+   private:
+      //Basic constuctor is private
+      ArraySlice () {}
+      const size_t fsize; //The size should not change
+      T * fbegin;
+      const T * fbeginConst;
+};
+
+template<typename T>
+std::ostream & operator << (std::ostream &os, const ArraySlice<T> & arr){
+   os<<"[";
+   for (size_t i = 0; i < arr.size()-1; ++i){
+      os<<arr[i]<<",";
+   }
+   os<<arr[arr.size()-1]<<"]";
+   return os;
+}
+#endif

source/CAT.h

+
+//#include"FITS.h"
+#include"galprop_classes.h" // since FITS already here and included in galplot
+
+class CAT
+{
+ public:
+
+
+
+  long n_sources;
+  double *ra  ,*dec  ;
+  double *longitude,*latitude;
+  char **name;
+
+  double *flux,*flux_err;                      //AWS20060110 
+
+  int read (char *directory);     
+  int read (char *directory,char*filename);    //AWS20060110     
+  void print();
+
+  void find (char *name_,
+             double *ra_, double *dec_, double *longitude_, double *latitude_);
+};

source/Configure.h

+
+//**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****|
+// * Configure.h *                                 galprop package * 4/14/2000 
+//**"****!****"****!****"****!****"****!****"****!****"****!****"****!****"****|
+
+#ifndef Configure_h
+#define Configure_h
+
+using namespace std;
+#include<iostream>//AWS20050624
+#include<cmath>   //AWS20050624
+#include<string> //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
+
+

source/Convolved.h

+//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<double> GLAST_convolved_counts;              //AWS20080604
+   Skymap<double> GLAST_convolved_counts_over_exposure;//AWS20080604
+
+   Skymap<double> GLAST_convolved_counts_HI;           //AWS20080606
+   Skymap<double> GLAST_convolved_counts_H2;           //AWS20080606
+   Skymap<double> GLAST_convolved_counts_IC;           //AWS20080606
+   Skymap<double> GLAST_convolved_counts_HII;          //AWS20090810
+
+
+   Skymap<double> GLAST_convolved_counts_sources;      //AWS20081209
+   Skymap<double> GLAST_convolved_intensity_sources;   //AWS20081212
+
+   Skymap<double> GLAST_convolved_counts_sourcepop_sublimit; //AWS20090109
+   Skymap<double> GLAST_convolved_counts_sourcepop_soplimit; //AWS20090109
+
+
+   Skymap<double> GLAST_convolved_intensity_sourcepop_sublimit; //AWS20090109
+   Skymap<double> GLAST_convolved_intensity_sourcepop_soplimit; //AWS20090109
+
+
+   Skymap<double> GLAST_convolved_counts_solar_IC;        //AWS20100215
+   Skymap<double> GLAST_convolved_intensity_solar_IC;     //AWS20100215
+
+   Skymap<double> GLAST_convolved_counts_solar_disk;      //AWS20100215
+   Skymap<double> GLAST_convolved_intensity_solar_disk;   //AWS20100215
+
+};

source/Coordinate.h

+#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

source/Counts.h

+#ifndef COUNTSMAP_H
+#define COUNTSMAP_H
+
+#include <string>
+#include <valarray>
+#include "Skymap.h"
+#include <exception>
+#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<double> & eMin, std::valarray<double> & 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<double> & eMin, std::valarray<double> & 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<double> &eMin, std::valarray<double> & 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<long> & getCountsMap() const;
+		const Skymap<int> & getCountsMap() const;                              //AWS20091002 
+		/** \brief Return a const reference to the lower energy boundaries */
+		const std::valarray<double> & getEMin() const;
+		/** \brief Return a const reference to the upper energy boundaries */
+		const std::valarray<double> & 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<SM::Coordinate, double> & photon);
+		//		Skymap<long> fDataSkymap;  //!< Stores the counts data
+		Skymap<int> fDataSkymap;  //!< Stores the counts data                    AWS20091002
+		std::valarray<double> 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

source/Data.h

+
+
+
+ 
+#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    *