From 2693d53c5f40c9965a71b9ab24809ff64eb81682 Mon Sep 17 00:00:00 2001
From: Andrew Strong <aws@olga4.opt.rzg.mpg.de>
Date: Sat, 30 Nov 2019 13:13:30 +0100
Subject: [PATCH] integration down in energy. v2.0
---
healpix_fermi_intensity_skymap.cc | 291 ++++++++++++++++++++++++++++--
1 file changed, 271 insertions(+), 20 deletions(-)
diff --git a/healpix_fermi_intensity_skymap.cc b/healpix_fermi_intensity_skymap.cc
index 2b7690b..aeae287 100644
--- a/healpix_fermi_intensity_skymap.cc
+++ b/healpix_fermi_intensity_skymap.cc
@@ -68,11 +68,12 @@ NB OpenMP is not actually used here, but is required when HEALPix libraries are
The maps can be viewed, rotated etc with the Aladin programme provided by CDS Strasbourg. Energies are selected by columns using 'Edit>Preferences'.
Note Aladin will read .gz files.
- Version 1.0
+ Version 2.0
A. W. Strong, MPE, 12 June 2019
+20191128 integration from low to high energies added
20190820 EBOUNDS extension for all intensity maps
20190730 moved to gitlab. Corrected exposure units and gtbin: algorithm=HEALPIX
20190612 first version from healpix_skymap_convert.cc
@@ -153,13 +154,21 @@ int main(int argc, char*argv[]){
// use prefix not suffix to preserve .fits if present in input file
string USING = "using_";
- string outfile_energies = prefix+"_healpix_intensity_" +USING +infile1;
- string outfile_energies_smoothed = prefix+"_healpix_intensity_smoothed_" +USING +infile1;
- string outfile_integrated_energies= prefix+"_healpix_intensity_integrated_" +USING +infile1;
- string infile_integrated_energies= prefix+"_healpix_intensity_integrated_" +USING +infile1;
- string outfile_integrated_energies_smoothed=prefix+"_healpix_intensity_integrated_smoothed_" +USING +infile1;
- string outfile_Alm =prefix+"_healpix_intensity_Alm_" +USING +infile1;
- string outfile_2_rebinned =prefix+"_healpix_exposure_rebinned_" +USING +infile2;
+ string outfile_energies = prefix+"_healpix_intensity_" +USING +infile1;
+ string outfile_energies_smoothed = prefix+"_healpix_intensity_smoothed_" +USING +infile1;
+
+ string outfile_integrated_energies_up= prefix+"_healpix_intensity_integrated_up_" +USING +infile1;
+ string infile_integrated_energies_up= prefix+"_healpix_intensity_integrated_up_" +USING +infile1;
+ string outfile_integrated_energies_up_smoothed= prefix+"_healpix_intensity_integrated_up_smoothed_" +USING +infile1;
+
+ string outfile_integrated_energies_down= prefix+"_healpix_intensity_integrated_down_" +USING +infile1;//AWS20191128
+ string infile_integrated_energies_down= prefix+"_healpix_intensity_integrated_down_" +USING +infile1;//AWS20191128
+ string outfile_integrated_energies_down_smoothed=prefix+"_healpix_intensity_integrated_down_smoothed_" +USING +infile1;//AWS20191128
+
+
+
+ string outfile_Alm =prefix+"_healpix_intensity_Alm_" +USING +infile1;
+ string outfile_2_rebinned =prefix+"_healpix_exposure_rebinned_" +USING +infile2;
double fwhm_deg=1.0;
@@ -263,6 +272,13 @@ int main(int argc, char*argv[]){
arr<double> integrated_energies_input_1_Emax;
integrated_energies_input_1_Emax.alloc(nenergy);
+ arr<double> integrated_energies_down_input_1_Emin; //AWS20191128
+ integrated_energies_down_input_1_Emin.alloc(nenergy); //AWS20191128
+
+ arr<double> integrated_energies_down_input_1_Emax; //AWS20191128
+ integrated_energies_down_input_1_Emax.alloc(nenergy); //AWS20191128
+
+
hdu=3;
inputhandle.goto_hdu(hdu);
colnum=1;
@@ -335,7 +351,11 @@ int main(int argc, char*argv[]){
for( j=0;j<nenergy;j++) integrated_energies_input_1_Emin[j]=energies_input_1_Emin[j]; // min energy as for bands
- for( j=0;j<nenergy;j++) integrated_energies_input_1_Emax[j]=energies_input_1_Emax[nenergy-1]; // upper energy always global max
+ for( j=0;j<nenergy;j++) integrated_energies_input_1_Emax[j]=energies_input_1_Emax[nenergy-1]; // upper energy is global max
+
+ for( j=0;j<nenergy;j++) integrated_energies_down_input_1_Emin[j]=energies_input_1_Emin[0]; // min energy is global min AWS20191128
+ for( j=0;j<nenergy;j++) integrated_energies_down_input_1_Emax[j]=energies_input_1_Emax[j]; // upper energy as for bands AWS20191128
+
cout<<"end of read counts"<<endl;
@@ -637,14 +657,14 @@ int main(int argc, char*argv[]){
// FURTHER OPERATIONS //
//////////////////////////////////////////////////////////////////////
- ///////////////////////////////////
- // energy integrated intensities //
- ///////////////////////////////////
+ //////////////////////////////////////
+ // energy integrated intensities up //
+ //////////////////////////////////////
- cout<<"==== energy integrated intensities "<<endl;
+ cout<<"==== energy integrated intensities up "<<endl;
- outfile=outfile_integrated_energies;
+ outfile=outfile_integrated_energies_up;
cout<<"writing to "<<outfile<<endl;
@@ -713,6 +733,86 @@ int main(int argc, char*argv[]){
out.close();
+ ////////////////////////////////////////
+ // energy integrated intensities down // AWS20191128
+ ////////////////////////////////////////
+
+ cout<<"==== energy integrated intensities down "<<endl;
+
+
+ outfile=outfile_integrated_energies_down;
+ cout<<"writing to "<<outfile<<endl;
+
+
+ out.create("!"+outfile);// ! to overwrite
+
+
+
+ prepare_Healpix_fitsmap(out,map_standard, datatype, colname);//healpix_map_fitsio.h
+
+ out.set_key(string("EXTNAME"),string("SKYMAP"));// extension HDU name
+
+ for( j=0;j<nenergy;j++)
+ {
+
+ for( i=0;i<npix;i++)
+ {
+ data[i]=0.;
+ for (int jj=0;jj<=j;jj++) data[i]+=val[i][jj];
+ }
+
+
+
+ colnum=j+1;
+ out.write_column(colnum,data); // fitshandle.h
+ }
+
+
+ // add EBOUNDS extension table with energies Emin Emax
+ // column names come from input counts file defined above
+
+ cout<<"adding EBOUNDS extension"<<endl;
+
+
+ ncols=3;
+
+ energy_table.resize(ncols);
+ for(colnum=1;colnum<=ncols;colnum++)
+ {
+ if(colnum==1) energy_table[colnum-1]=fitscolumn(input_colname[colnum-1]," ",1,PLANCK_INT32 ); //CHANNEL
+ if(colnum==2|colnum==3) energy_table[colnum-1]=fitscolumn(input_colname[colnum-1],"MeV",1,PLANCK_FLOAT64); //EMIN, EMAX
+ }
+
+ out.insert_bintab(energy_table,"EBOUNDS" );
+
+ // hdu=3;
+ // out.goto_hdu (hdu); not required, automatic after insert_bintab
+ // CHANNEL column for compatibility with counts format
+
+
+ channel.alloc(nenergy);
+ cout<<"nenergy="<<nenergy<<endl;
+
+ for(j=0;j<nenergy;j++) channel[j]=j+1;
+
+ colnum=1;
+ out.write_column(colnum,channel);
+
+ colnum=2;
+ out.write_column(colnum,integrated_energies_down_input_1_Emin);
+
+ colnum=3;
+ out.write_column(colnum,integrated_energies_down_input_1_Emax);
+
+
+
+ out.close();
+
+
+
+ ////////////////////////
+ // smoothing all maps //
+ ////////////////////////
//////////////////////////////////////////////////////////////////////////
// energy bins with smoothing, using already converted format as input //
@@ -873,13 +973,13 @@ int main(int argc, char*argv[]){
- ////////////////////////////////////////////////////////////////////////////////
- // energy integrated with smoothing, using already converted format as input //
- ////////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////////////
+ // energy integrated up with smoothing, using already converted format as input //
+ //////////////////////////////////////////////////////////////////////////////////
{
- cout<<"==== smoothing fluxes integrated over energy"<<endl;
+ cout<<"==== smoothing fluxes integrated over energy up"<<endl;
cout<<" Gaussian FWHM = "<<fwhm_deg<< " deg"<<endl;
@@ -887,7 +987,7 @@ int main(int argc, char*argv[]){
Healpix_Map<double> map;
Healpix_Map<double> map_RING_out(order,RING);
- outfile= outfile_integrated_energies_smoothed;
+ outfile= outfile_integrated_energies_up_smoothed;
cout<<"writing to "<<outfile<<endl;
out.create("!"+outfile); // ! to overwrite
@@ -895,7 +995,7 @@ int main(int argc, char*argv[]){
out.set_key(string("EXTNAME"),string("SKYMAP"));// extension HDU name
- infile=outfile_integrated_energies;
+ infile=outfile_integrated_energies_up;
int ncolnum=0;
cout<< "reading "<<infile<<endl;
@@ -1021,6 +1121,157 @@ int main(int argc, char*argv[]){
}
+
+ /////////////////////////////////////////////////////////////////////////////////////
+ // energy integrated down with smoothing, using already converted format as input // AWS20191128
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ {
+
+ cout<<"==== smoothing fluxes integrated over energy down"<<endl;
+
+
+ cout<<" Gaussian FWHM = "<<fwhm_deg<< " deg"<<endl;
+
+ Healpix_Map<double> map;
+ Healpix_Map<double> map_RING_out(order,RING);
+
+ outfile= outfile_integrated_energies_down_smoothed;
+ cout<<"writing to "<<outfile<<endl;
+
+ out.create("!"+outfile); // ! to overwrite
+ prepare_Healpix_fitsmap(out,map_RING_out, datatype, colname);
+
+ out.set_key(string("EXTNAME"),string("SKYMAP"));// extension HDU name
+
+ infile=outfile_integrated_energies_down;
+ int ncolnum=0;
+
+ cout<< "reading "<<infile<<endl;
+
+
+
+ fitshandle in;
+ in.open(infile);
+ in.goto_hdu(2);
+
+
+ ncolnum=in.ncols();
+ cout<<" number of columns ="<<ncolnum<<endl;
+
+
+ for (colnum=1;colnum<=ncolnum;colnum++)
+
+ {
+ if(debug==1) cout<< "reading "<<infile<<", column number: "<< colnum;
+
+ read_Healpix_map_from_fits(infile,map,colnum,2);// HDU 2 (1=primary)
+
+ if(debug==1)
+ {
+ cout<<" Npix = "<<map.Npix();
+ cout<<" Nside = "<<map.Nside();
+ cout<<" Order = "<<map.Order();
+ cout<<" Scheme= "<<map.Scheme()<<endl; // 0 = RING, 1 = NESTED
+ }
+
+ /*
+ // convert to RING (WMAP is NESTED, galprop is RING already)
+ Healpix_Map<double> map_RING(map.Order(),RING); // RING needed for Alm
+ map_RING.Import_nograde(map);
+ */
+
+ // write_Healpix_map_to_fits(outfile,map_RING , datatype);
+
+ int lmax=3*map.Nside();// too large values give wrong results, this is recommended value
+ int mmax=3*map.Nside();
+
+ Alm<xcomplex<double> > alm(lmax,mmax);
+
+ arr<double> weight;
+ weight.alloc(2*map.Nside());// must have at least 2*map.Nside() entries
+ weight.fill(1.0);
+
+ bool add_alm=false;
+
+ // map2alm(map,alm,weight,add_alm); // NB consider iter version
+ // (0 equivalent to map2alm)
+ int num_iter = 3; // AWS20111102
+ map2alm_iter(map,alm,num_iter,weight); // AWS20111102
+
+ double dtr=acos(-1.0)/180.;
+
+
+ double fwhm = fwhm_deg*dtr;// radian !!
+ smoothWithGauss(alm,fwhm);
+
+ // filter
+ int filter = 0;
+ int lfilter_min=0;
+ int lfilter_max=30;
+ int mfilter_min=0;
+ int mfilter_max=1000;
+
+ // for m<0 defined as (-1)^m Al|m| so not explicitly required ? Healpix manual A.4 eq. (22)
+ if(filter==1)
+ for(int l=0;l<=lmax;l++)
+ for(int m=0;m<=l; m++)if(l<lfilter_min || l>lfilter_max || m<mfilter_min || m>mfilter_max) alm(l,m)=0.;
+
+ alm2map(alm,map_RING_out);
+
+ cout<<"energy integrated map average before smoothing = "<< map.average() << " after smoothing = " <<map_RING_out.average()
+ <<" after/before = "<<map_RING_out.average()/map.average() <<endl;
+
+ out.write_column(colnum,map_RING_out.Map()); // fitshandle.h
+
+
+
+
+ }
+
+
+// add EBOUNDS extension table with energies Emin Emax
+ // column names come from input counts file defined above
+
+ cout<<"adding EBOUNDS extension"<<endl;
+
+
+ ncols=3;
+
+ energy_table.resize(ncols);
+ for(colnum=1;colnum<=ncols;colnum++)
+ {
+ if(colnum==1) energy_table[colnum-1]=fitscolumn(input_colname[colnum-1]," ",1,PLANCK_INT32 ); //CHANNEL
+ if(colnum==2|colnum==3) energy_table[colnum-1]=fitscolumn(input_colname[colnum-1],"MeV",1,PLANCK_FLOAT64); //EMIN, EMAX
+ }
+
+ out.insert_bintab(energy_table,"EBOUNDS" );
+
+ // hdu=3;
+ // out.goto_hdu (hdu); not required, automatic after insert_bintab
+ // CHANNEL column for compatibility with counts format
+
+
+ channel.alloc(nenergy);
+ cout<<"nenergy="<<nenergy<<endl;
+
+ for(j=0;j<nenergy;j++) channel[j]=j+1;
+
+ colnum=1;
+ out.write_column(colnum,channel);
+
+ colnum=2;
+ out.write_column(colnum,integrated_energies_down_input_1_Emin);
+
+ colnum=3;
+ out.write_column(colnum,integrated_energies_down_input_1_Emax);
+
+ out.close();
+
+}
+
+
+
} // if(0) file 1 AWS20190906
cout<<endl<<" ==== healpix_fermi_intensity_skymap complete"<<endl;
--
GitLab