-
Andrew Strong authored
gen_source_population.cc: if(1) generate Fermi catalogue to file, for possible control later. start and end to cout
Andrew Strong authoredgen_source_population.cc: if(1) generate Fermi catalogue to file, for possible control later. start and end to cout
gen_source_population.cc 40.78 KiB
#include"Galplot.h"
#include"galprop_classes.h"
#include"galprop.h"
#include"constants.h"
#include"fitsio.h"
int Galplot::gen_source_population()
{
int status=0;
int ip,ii;
unsigned seed;
double logL;//AWS20060119
double logL_sensitivity_file; //AWS20200712
cout<<" >>>> gen_source_population "<<endl;
txt_stream<<" >>>> gen_source_population "<<endl;
///////////////////////////////////////////////////////////////////////
cout<<" sourcepop1"<<endl;
// init needs the sensitivity file to read AWS20170113
sourcepop1.Fermi_sensitivity_file = new char[200]; //AWS20170113
strcpy(sourcepop1.Fermi_sensitivity_file,galplotdef.Fermi_sensitivity_file); //AWS20170113
sourcepop1.init();
strcpy(sourcepop1.title,"sourcepop1_population_synthesis"); //AWS20140201 used in filename so use underscores
strcat(sourcepop1.title,"_"); strcat(sourcepop1.title,galplotdef.galdef_ID); //AWS20140202
strcat(sourcepop1.title,"_"); strcat(sourcepop1.title,galplotdef.psfile_tag); //AWS20140201
strcpy(sourcepop1.galdef_ID, galplotdef.galdef_ID); //AWS20200107
strcpy(sourcepop1.psfile_tag,galplotdef.psfile_tag); //AWS20200107
sourcepop1.txt_stream=&txt_stream; //NB can only copy as pointer
// these parameters not from galplotdef file
sourcepop1.luminosity_function_type=2;
// used only for profiles in SourcePopulation
sourcepop1.longprof_b_min=galplotdef. lat_min1;
sourcepop1.longprof_b_max=galplotdef. lat_max2;
// SourcePopulation uses -180<l<180
sourcepop1. latprof_l_min=galplotdef. long_min1-180.;
sourcepop1. latprof_l_max=galplotdef. long_max2-180.;
// used for skymap so must have value from galaxy
sourcepop1.longprof_dlong=galaxy. d_long;
sourcepop1. latprof_dlat =galaxy. d_lat;
// allow for FITS convention
// in galplotdef
sourcepop1.long_min1=galplotdef.long_min1 - galdef.d_long/2.;
sourcepop1.long_min2=galplotdef.long_min2 - galdef.d_long/2;
sourcepop1.long_max1=galplotdef.long_max1 + galdef.d_long/2;
sourcepop1.long_max2=galplotdef.long_max2 + galdef.d_long/2;
sourcepop1. lat_min1=galplotdef. lat_min1 - galdef.d_lat /2.;
sourcepop1. lat_min2=galplotdef. lat_min2 - galdef.d_lat /2.;
sourcepop1. lat_max1=galplotdef. lat_max1 + galdef.d_lat /2.;
sourcepop1. lat_max2=galplotdef. lat_max2 + galdef.d_lat /2.;
sourcepop1.L=1.e34;
sourcepop1.E_ref_low = galplotdef.sourcepop1_E_ref_low ; //AWS20060109
sourcepop1.E_ref_high = galplotdef.sourcepop1_E_ref_high; //AWS20060109
sourcepop1.n_E = galaxy.n_E_gammagrid;
sourcepop1. E = new double[sourcepop1.n_E];
for (ip=0; ip<sourcepop1.n_E;ip++) sourcepop1.E[ip]=galaxy.E_gamma[ip];
read_EGRET_data();//AWS20111013
cout<<"data.n_E_EGRET="<<data.n_E_EGRET<<endl;
sourcepop1.n_E_bin = data.n_E_EGRET;
if(data.n_E_EGRET>=10)//AWS20111013 avoid segfault if EGRET data not read AWS20140114: was ==13 but should work for other values
{
sourcepop1.E_bin_low =new double[ sourcepop1.n_E_bin ];
sourcepop1.E_bin_high =new double[ sourcepop1.n_E_bin ];
for (ip=0; ip< data.n_E_EGRET; ip++)
{
cout<<"data.n_E_EGRET="<<data.n_E_EGRET<<" ip="<<ip<<endl;
sourcepop1.E_bin_low [ip]=data.E_EGRET[ip ];
sourcepop1.E_bin_high[ip]=data.E_EGRET[ip+1];
}
}//if
// these parameters from galplotdef file
sourcepop1.verbose =galplotdef.sourcepop1_verbose;
sourcepop1.density0 =galplotdef.sourcepop1_density0;
sourcepop1.oversample=galplotdef.sourcepop1_oversample; //AWS20060118
sourcepop1.L_min =galplotdef.sourcepop1_L_min;
sourcepop1.L_max =galplotdef.sourcepop1_L_max;
sourcepop1.alpha_L =galplotdef.sourcepop1_alpha_L;
sourcepop1.flux_detection_limit
=galplotdef.sourcepop1_fluxlimit;
sourcepop1.alpha_R= galplotdef.sourcepop1_alpha_R;
sourcepop1. beta_R= galplotdef.sourcepop1_beta_R;
sourcepop1.zscale = galplotdef.sourcepop1_zscale;
sourcepop1.alpha_z= galplotdef.sourcepop1_alpha_z;
sourcepop1.spectral_model = galplotdef.sourcepop1_spectral_model;//AWS20100824
sourcepop1.spectrum_g= galplotdef.sourcepop1_spectrumg;
sourcepop1.spectrum_norm_E=galplotdef.sourcepop1_specnormE;
sourcepop1.spectrum_g_0= galplotdef.sourcepop1_spect_g_0; //AWS20051219
sourcepop1.spectrum_g_1= galplotdef.sourcepop1_spect_g_1; //AWS20051219
sourcepop1.spectrum_g_2= galplotdef.sourcepop1_spect_g_2; //AWS20051219
sourcepop1.spectrum_br0= galplotdef.sourcepop1_spect_br0; //AWS20051219
sourcepop1.spectrum_br1= galplotdef.sourcepop1_spect_br1; //AWS20051219
sourcepop1.spectrum_g_0_sigma=galplotdef.sourcepop1_spectrum_g_0_sigma; //AWS20100824
sourcepop1.spectrum_br0_sigma=galplotdef.sourcepop1_spectrum_br0_sigma; //AWS20100824
sourcepop1.spectrum_g_1_sigma=galplotdef.sourcepop1_spectrum_g_1_sigma; //AWS20100824
sourcepop1.spectrum_br1_sigma=galplotdef.sourcepop1_spectrum_br1_sigma; //AWS20100824
sourcepop1.spectrum_g_2_sigma=galplotdef.sourcepop1_spectrum_g_2_sigma; //AWS20100824
// same order as galprop maps AWS20131124
// sourcepop1.healpix_skymap_order = galaxy.pi0_decay_hp_skymap.Order(); //AWS20131124
// in case this routine is called before galprop maps are read, use Galdef value //AWS20140127
sourcepop1.healpix_skymap_order = galdef.healpix_order; //AWS20140127
cout<<"gen_source_population: skymap order as input galprop maps: order = "<< sourcepop1.healpix_skymap_order<<endl; //AWS20131124
sourcepop1.print();
seed=1234;
//---------------------------------------------------------
// test MSPs
// int MSPs=0;
// if(galplotdef.sourcepop_ext_model==2) MSPs=1; //AWS20120530
// removed at r1247
//----------------------------------------------------------
cout<<"Fermi catalogue 1 analysis: sourcepop4"<<endl;
txt_stream<<"Fermi catalogue 1 analysis: sourcepop4"<<endl;
// init needs the sensitivity file to read. AWS20170113
// "=" does not copy the char it seems so do it explicitly instead of putting "=" before init AWS20170113
sourcepop4.Fermi_sensitivity_file = new char[200]; //AWS20170113
strcpy(sourcepop4.Fermi_sensitivity_file,galplotdef.Fermi_sensitivity_file); //AWS20170113
sourcepop4.init(); //AWS20081215
sourcepop4=sourcepop1;
strcpy(sourcepop4.title,"sourcepop4_Fermi_Catalogue_1"); //AWS20140201 used in filename so use underscores
strcat(sourcepop4.title,"_"); strcat(sourcepop4.title,galplotdef.galdef_ID); //AWS20140202
strcat(sourcepop4.title,"_"); strcat(sourcepop4.title,galplotdef.psfile_tag); //AWS20140201
strcpy(sourcepop4.galdef_ID, galplotdef.galdef_ID); //AWS20200107
strcpy(sourcepop4.psfile_tag,galplotdef.psfile_tag); //AWS20200107
sourcepop4.n_E_bin = data.GLAST_counts_healpix .getEMin().size(); //AWS20090105
sourcepop4.E_bin_low = new double[ sourcepop4.n_E_bin ]; //AWS20090105
sourcepop4.E_bin_high = new double[ sourcepop4.n_E_bin ]; //AWS20090105
for (ip=0; ip<sourcepop4.n_E_bin ; ip++) //AWS20090105
{
sourcepop4.E_bin_low [ip]= data.GLAST_counts_healpix.getEMin()[ip]; //AWS20090105
sourcepop4.E_bin_high[ip]= data.GLAST_counts_healpix.getEMax()[ip]; //AWS20090105
}
sourcepop4.verbose=4;
// Fermi 2FGL has only >100 MeV and > 1000 MeV, so choose on basis of population synthesis value
int options= 100; // use flux> 100 MeV //AWS20110808
if(galplotdef.sourcepop1_E_ref_low > 999.) //AWS20110830
options=1000; // use flux>1000 MeV
if(galplotdef.sourcepop1_E_ref_low >9999.) //AWS20120503
options=10000; // use flux>10000 MeV
cout<<" sourcepop4: using options = "<<options<< " MeV since galplotdef.sourcepop1_E_ref_low = "<< galplotdef.sourcepop1_E_ref_low<<endl;
sourcepop4.gen_sample_Fermi_catalogue(configure.fits_directory,galplotdef.Fermi_cat_file, options); //AWS20110808
cout<<"after sourcepop4.gen_sample_Fermi_catalogue"<<endl;
sourcepop4.print();
sourcepop4.IDL_plot_control=1; //AWS20200623 first plot, no oplot
sourcepop4.IDL_psym =6; //AWS20200623 open squares no line
sourcepop4.analyse_sample();
//----------------------------------------------------------
// additional catalogues AWS20140114
cout<<"Fermi catalogue 2 analysis: sourcepop6"<<endl;
txt_stream<<"Fermi catalogue 2 analysis: sourcepop6"<<endl;
// init needs the sensitivity file to read. AWS20170113
// "=" does not copy the char it seems so do it explicitly instead of putting "=" before init AWS20170113
sourcepop6.Fermi_sensitivity_file = new char[200]; //AWS20170113
strcpy(sourcepop6.Fermi_sensitivity_file,galplotdef.Fermi_sensitivity_file); //AWS20170113
sourcepop6.init();
sourcepop6=sourcepop1;
strcpy(sourcepop6.title,"sourcepop6_Fermi_Catalogue_2"); //AWS20140201 used in filename so use underscores
strcat(sourcepop6.title,"_"); strcat(sourcepop6.title,galplotdef.galdef_ID); //AWS20140202
strcat(sourcepop6.title,"_"); strcat(sourcepop6.title,galplotdef.psfile_tag); //AWS20140201
strcpy(sourcepop6.galdef_ID, galplotdef.galdef_ID); //AWS20200107
strcpy(sourcepop6.psfile_tag,galplotdef.psfile_tag); //AWS20200107
sourcepop6.n_E_bin = data.GLAST_counts_healpix .getEMin().size();
sourcepop6.E_bin_low = new double[ sourcepop6.n_E_bin ];
sourcepop6.E_bin_high = new double[ sourcepop6.n_E_bin ];
for (ip=0; ip<sourcepop6.n_E_bin ; ip++)
{
sourcepop6.E_bin_low [ip]= data.GLAST_counts_healpix.getEMin()[ip];
sourcepop6.E_bin_high[ip]= data.GLAST_counts_healpix.getEMax()[ip];
}
sourcepop6.verbose=4;
// Fermi 2FGL has only >100 MeV and > 1000 MeV, so choose on basis of population synthesis value
options= 100; // use flux> 100 MeV
if(galplotdef.sourcepop1_E_ref_low > 999.) //AWS20110830
options=1000; // use flux>1000 MeV
if(galplotdef.sourcepop1_E_ref_low >9999.)
options=10000; // use flux>10000 MeV
cout<<" sourcepop6: using options = "<<options<< " MeV since galplotdef.sourcepop1_E_ref_low = "<< galplotdef.sourcepop1_E_ref_low<<endl;
sourcepop6.gen_sample_Fermi_catalogue(configure.fits_directory,galplotdef.Fermi_cat_file2, options);
cout<<"after sourcepop6.gen_sample_Fermi_catalogue"<<endl;
sourcepop6.print();
sourcepop6.IDL_plot_control=2; //AWS20200623 oplot on first plot
sourcepop6.IDL_psym =5; //AWS20200623 open triangles no line
sourcepop6.analyse_sample();
// additional catalogue AWS20140115
cout<<"Fermi catalogue 3 analysis: sourcepop7"<<endl;
txt_stream<<"Fermi catalogue 3 analysis: sourcepop7"<<endl;
// init needs the sensitivity file to read. AWS20170113
// "=" does not copy the char it seems so do it explicitly instead of putting "=" before init AWS20170113
sourcepop7.Fermi_sensitivity_file = new char[200]; //AWS20170113
strcpy(sourcepop7.Fermi_sensitivity_file,galplotdef.Fermi_sensitivity_file); //AWS20170113
sourcepop7.init();
sourcepop7=sourcepop1;
strcpy(sourcepop7.title,"sourcepop7_Fermi_Catalogue_3"); //AWS20140201 used in filename so use underscores
strcat(sourcepop7.title,"_"); strcat(sourcepop7.title,galplotdef.galdef_ID); //AWS20140202
strcat(sourcepop7.title,"_"); strcat(sourcepop7.title,galplotdef.psfile_tag); //AWS20140201
strcpy(sourcepop7.galdef_ID, galplotdef.galdef_ID); //AWS20200107
strcpy(sourcepop7.psfile_tag,galplotdef.psfile_tag); //AWS20200107
sourcepop7.n_E_bin = data.GLAST_counts_healpix .getEMin().size();
sourcepop7.E_bin_low = new double[ sourcepop7.n_E_bin ];
sourcepop7.E_bin_high = new double[ sourcepop7.n_E_bin ];
for (ip=0; ip<sourcepop7.n_E_bin ; ip++)
{
sourcepop7.E_bin_low [ip]= data.GLAST_counts_healpix.getEMin()[ip];
sourcepop7.E_bin_high[ip]= data.GLAST_counts_healpix.getEMax()[ip];
}
sourcepop7.verbose=4;
// Fermi 2FGL has only >100 MeV and > 1000 MeV, so choose on basis of population synthesis value
options= 100; // use flux> 100 MeV
if(galplotdef.sourcepop1_E_ref_low > 999.) //AWS20110830
options=1000; // use flux>1000 MeV
if(galplotdef.sourcepop1_E_ref_low >9999.)
options=10000; // use flux>10000 MeV
cout<<" sourcepop7: using options = "<<options<< " MeV since galplotdef.sourcepop1_E_ref_low = "<< galplotdef.sourcepop1_E_ref_low<<endl;
sourcepop7.gen_sample_Fermi_catalogue(configure.fits_directory,galplotdef.Fermi_cat_file3, options);
cout<<"after sourcepop7.gen_sample_Fermi_catalogue"<<endl;
sourcepop7.print();
sourcepop7.IDL_plot_control =2; //AWS20200623 oplot on first plot
sourcepop7.IDL_psym =4; //AWS20200623 open diamonds no line
sourcepop7.analyse_sample();
//----------------------------------------------------------
cout<<"Fermi population synthesis: sourcepop5"<<endl; //AWS20090107
txt_stream<<"Fermi population synthesis: sourcepop5"<<endl; //AWS20090107
// init needs the sensitivity file to read. AWS20170113
// "=" does not copy the char it seems so do it explicitly instead of putting "=" before init AWS20170113
sourcepop5.Fermi_sensitivity_file = new char[200]; //AWS20170113
strcpy(sourcepop5.Fermi_sensitivity_file,galplotdef.Fermi_sensitivity_file); //AWS20170113
sourcepop5.init(); //AWS20080107
sourcepop5=sourcepop1;
strcpy(sourcepop5.title,"sourcepop5_Fermi_population_synthesis"); //AWS20140201 used in filename so use underscores
strcat(sourcepop5.title,"_"); strcat(sourcepop5.title,galplotdef.galdef_ID); //AWS20140202
strcat(sourcepop5.title,"_"); strcat(sourcepop5.title,galplotdef.psfile_tag); //AWS20140201
strcpy(sourcepop5.galdef_ID, galplotdef.galdef_ID); //AWS20200107
strcpy(sourcepop5.psfile_tag,galplotdef.psfile_tag); //AWS20200107
sourcepop5.n_E_bin = data.GLAST_counts_healpix .getEMin().size();
sourcepop5.E_bin_low = new double[ sourcepop5.n_E_bin ];
sourcepop5.E_bin_high = new double[ sourcepop5.n_E_bin ];
for (ip=0; ip<sourcepop5.n_E_bin ; ip++)
{
sourcepop5.E_bin_low [ip]= data.GLAST_counts_healpix.getEMin()[ip];
sourcepop5.E_bin_high[ip]= data.GLAST_counts_healpix.getEMax()[ip];
}
sourcepop5.verbose=4;
cout<<"before sourcepop5.gen_sample_Fermi_catalogue"<<endl;
sourcepop5.gen_sample(seed);
cout<<"after sourcepop5.gen_sample_Fermi_catalogue"<<endl;
sourcepop5.print();
sourcepop5.IDL_plot_control=3; //AWS20200714 oplot on first plot model NS, soplimit, sublimit
sourcepop5.IDL_psym =0; //AWS20200714 just line
sourcepop5.analyse_sample();
char Fermi_catalogue_from_popsyn_sample_file[500]; //AWS20111011
strcpy(Fermi_catalogue_from_popsyn_sample_file,"plots/");
strcat(Fermi_catalogue_from_popsyn_sample_file,"Fermi_catalogue_from_popsyn_sample_");
strcat(Fermi_catalogue_from_popsyn_sample_file,galdef.galdef_ID);
strcat(Fermi_catalogue_from_popsyn_sample_file,"_" );
strcat(Fermi_catalogue_from_popsyn_sample_file,galplotdef.psfile_tag);
strcat(Fermi_catalogue_from_popsyn_sample_file,".fits");
cout<<"starting gen_Fermi_catalogue_from_sample"<<endl;
if(1) //AWS20200721
sourcepop5.gen_Fermi_catalogue_from_sample(configure.fits_directory,Fermi_catalogue_from_popsyn_sample_file, options); //AWS20110829
cout<<"finished gen_Fermi_catalogue_from_sample"<<endl;
// copy healpix skymaps to standard place for use in Fermi profiles.
unconvolved.GLAST_unconvolved_intensity_sourcepop_sublimit = sourcepop5.healpix_skymap_intensity_sublimit;
unconvolved.GLAST_unconvolved_intensity_sourcepop_soplimit = sourcepop5.healpix_skymap_intensity_soplimit;
// convolved=unconvolved at this stage until convolution is implemented
convolved. GLAST_convolved_intensity_sourcepop_sublimit = unconvolved.GLAST_unconvolved_intensity_sourcepop_sublimit;
convolved. GLAST_convolved_intensity_sourcepop_soplimit = unconvolved.GLAST_unconvolved_intensity_sourcepop_soplimit;
//---------------------------------------------------------
/*AWS20110818
txt_stream<<endl<<"Comparison source synthesis : EGRET catalogue"<<endl;
txt_stream<<"source counts for region "<<sourcepop1.long_min1 <<" < l < "<<sourcepop1.long_max1 <<", "<<sourcepop1.long_min2 <<" < l < "<<sourcepop1.long_max2<<" / "
<<sourcepop1. lat_min1 <<" < b < "<<sourcepop1.lat_max1 <<", "<< sourcepop1. lat_min2 <<" < b < "<<sourcepop1. lat_max2<<endl;
sourcepop1.print(); sourcepop2.print();
logL=0.0;
for (ii=0;ii< sourcepop1.n_dlnN_dlnS;ii++)
{
if(pow(10, sourcepop1.lnS_min+ii*sourcepop1.dlnS)>=sourcepop1.flux_detection_limit
&&sourcepop1.dlnN_dlnS[ii]>0.0)
logL+= - sourcepop1.dlnN_dlnS [ii] + sourcepop2.dlnN_dlnS [ii] * log(sourcepop1.dlnN_dlnS[ii]);
}
for (ii=0;ii< sourcepop1.n_dlnN_dlnS;ii++)
{
if( sourcepop1.dlnN_dlnS[ii]>0 || sourcepop2.dlnN_dlnS[ii]>0 )
txt_stream <<"ii="<<ii
<<" lnS="<< sourcepop1.lnS_min+ii*sourcepop1.dlnS <<"/"<< sourcepop1.lnS_min+(ii+1)*sourcepop1.dlnS
<<" S=" <<pow(10, sourcepop1.lnS_min+ii*sourcepop1.dlnS) <<"/"<<pow(10, sourcepop1.lnS_min+(ii+1)*sourcepop1.dlnS)
<<" N(S)= "<< sourcepop1.dlnN_dlnS [ii]<<":" << sourcepop2.dlnN_dlnS [ii]
<<" N(>S)= "<< sourcepop1.dlnN_dlnS_int[ii]<<":" << sourcepop2.dlnN_dlnS_int[ii]
<<" F(S)= "<< sourcepop1.FS [ii]<<":" << sourcepop2.FS [ii]
<<" F(>S)= "<< sourcepop1.FS_int [ii]<<":" << sourcepop2.FS_int [ii]
<<endl;
}
txt_stream <<"log likelihood "<< logL<<endl;
*/
//---------------------------------------------------------
txt_stream<<endl<<"Comparison source synthesis : Fermi catalogue 1"<<endl; //AWS20081216
txt_stream<<"source counts for region "<<sourcepop5.long_min1 <<" < l < "<<sourcepop5.long_max1 <<", "<<sourcepop5.long_min2 <<" < l < "<<sourcepop5.long_max2<<" / "
<<sourcepop5. lat_min1 <<" < b < "<<sourcepop5.lat_max1 <<", "<< sourcepop5. lat_min2 <<" < b < "<<sourcepop5. lat_max2<<endl; // AWS20111011
sourcepop5.print(); sourcepop4.print();// AWS20111011
logL=0.0;
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++) // AWS20111011
{
if(pow(10, sourcepop5.lnS_min+ii*sourcepop5.dlnS)>=sourcepop5.flux_detection_limit
&&sourcepop5.dlnN_dlnS[ii]>0.0) // AWS20111011
logL+= - sourcepop5.dlnN_dlnS [ii] + sourcepop4.dlnN_dlnS [ii] * log(sourcepop5.dlnN_dlnS[ii]);// AWS20111011
}
logL_sensitivity_file=0.0;
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++) // AWS20200712
{
if(sourcepop5.dlnN_dlnS_soplimit[ii]>0.0)
logL_sensitivity_file+= - sourcepop5.dlnN_dlnS_soplimit [ii] + sourcepop4.dlnN_dlnS [ii] * log(sourcepop5.dlnN_dlnS_soplimit[ii]);
}
txt_stream << "sourcepop5.n_dlnN_dlnS="<< sourcepop5.n_dlnN_dlnS<<endl;//AWS20111011
txt_stream << "sourcepop4.n_dlnN_dlnS="<< sourcepop4.n_dlnN_dlnS<<endl;//AWS20111011
txt_stream << "sourcepop5 model all sources, Fermi catalogue 1 "<<endl;
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++)
{
if( sourcepop5.dlnN_dlnS[ii]>0 || sourcepop4.dlnN_dlnS[ii]>0 )
txt_stream <<"ii="<<ii
<<" lnS="<< sourcepop5.lnS_min+ii*sourcepop5.dlnS <<"/"<< sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS
<<" S=" <<pow(10, sourcepop5.lnS_min+ii*sourcepop5.dlnS) <<"/"<<pow(10, sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS)
<<" N(S)= "<< sourcepop5.dlnN_dlnS [ii]<<":" << sourcepop4.dlnN_dlnS [ii]
<<" N(>S)= "<< sourcepop5.dlnN_dlnS_int[ii]<<":" << sourcepop4.dlnN_dlnS_int[ii]
<<" F(S)= "<< sourcepop5.FS [ii]<<":" << sourcepop4.FS [ii]
<<" F(>S)= "<< sourcepop5.FS_int [ii]<<":" << sourcepop4.FS_int [ii]
<<endl;
}
txt_stream <<endl;
txt_stream << "sourcepop5: model sources above threshold, sourcepop4: Fermi catalogue 1"<<endl;
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++)
{
if( sourcepop5.dlnN_dlnS[ii]>0 || sourcepop4.dlnN_dlnS[ii]>0 )
txt_stream <<"ii="<<ii
<<" lnS="<< sourcepop5.lnS_min+ii*sourcepop5.dlnS <<"/"<< sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS
<<" S=" <<pow(10, sourcepop5.lnS_min+ii*sourcepop5.dlnS) <<"/"<<pow(10, sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS)
<<" N(S)= "<< sourcepop5.dlnN_dlnS_soplimit [ii]<<":" << sourcepop4.dlnN_dlnS [ii]
<<" N(>S)= "<< sourcepop5.dlnN_dlnS_int_soplimit[ii]<<":" << sourcepop4.dlnN_dlnS_int[ii]
<<" F(S)= "<< sourcepop5.FS_soplimit [ii]<<":" << sourcepop4.FS [ii]
<<" F(>S)= "<< sourcepop5.FS_int_soplimit [ii]<<":" << sourcepop4.FS_int [ii]
<<endl;
}
txt_stream <<endl;
txt_stream << "sourcepop5: model sources below threshold, sourcepop4: Fermi catalogue 1"<<endl;
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++)
{
if( sourcepop5.dlnN_dlnS[ii]>0 || sourcepop4.dlnN_dlnS[ii]>0 )
txt_stream <<"ii="<<ii
<<" lnS="<< sourcepop5.lnS_min+ii*sourcepop5.dlnS <<"/"<< sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS
<<" S=" <<pow(10, sourcepop5.lnS_min+ii*sourcepop5.dlnS) <<"/"<<pow(10, sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS)
<<" N(S)= "<< sourcepop5.dlnN_dlnS_sublimit [ii]<<":" << sourcepop4.dlnN_dlnS [ii]
<<" N(>S)= "<< sourcepop5.dlnN_dlnS_int_sublimit[ii]<<":" << sourcepop4.dlnN_dlnS_int[ii]
<<" F(S)= "<< sourcepop5.FS_sublimit [ii]<<":" << sourcepop4.FS [ii]
<<" F(>S)= "<< sourcepop5.FS_int_sublimit [ii]<<":" << sourcepop4.FS_int [ii]
<<endl;
}
// one-line summary statistics
txt_stream << " galdef_ID = " << galplotdef.galdef_ID;
txt_stream << " psfile_tag = " << galplotdef.psfile_tag;
txt_stream << " sourcepop4 vs sourcepop5: cat 1 versus model: ";
txt_stream << " density = " << sourcepop5.density0;
txt_stream << " L_min = " << sourcepop5.L_min;
txt_stream << " L_max = " << sourcepop5.L_max;
txt_stream << " alpha_L = " << sourcepop5.alpha_L;
txt_stream << " long_min1=" << galplotdef.long_min1;
txt_stream << " long_max1=" << galplotdef.long_max1;
txt_stream << " long_min2=" << galplotdef.long_min2;
txt_stream << " long_max2=" << galplotdef.long_max2;
txt_stream << " lat_min1=" << galplotdef.lat_min1;
txt_stream << " lat_max1=" << galplotdef.lat_max1;
txt_stream << " lat_min2=" << galplotdef.lat_min2;
txt_stream << " lat_max2=" << galplotdef.lat_max2;
txt_stream <<" log likelihood using user-defined sensitivity = "<< logL;
txt_stream <<" log likelihood using sensitivity file = " << logL_sensitivity_file<<endl;
// short form
txt_stream << " " << sourcepop5.density0;
txt_stream << " " << sourcepop5.L_min;
txt_stream << " " << sourcepop5.L_max;
txt_stream << " " << sourcepop5.alpha_L;
txt_stream << " " << galplotdef.long_min1;
txt_stream << " " << galplotdef.long_max1;
txt_stream << " " << galplotdef.long_min2;
txt_stream << " " << galplotdef.long_max2;
txt_stream << " " << galplotdef.lat_min1;
txt_stream << " " << galplotdef.lat_max1;
txt_stream << " " << galplotdef.lat_min2;
txt_stream << " " << galplotdef.lat_max2;
txt_stream <<" " << logL;
txt_stream <<" " << logL_sensitivity_file;
txt_stream <<" logL statistics" <<endl;
//---------------------------------------------------------
txt_stream<<endl<<"Comparison source synthesis : Fermi catalogue 2"<<endl;
txt_stream<<"source counts for region "<<sourcepop5.long_min1 <<" < l < "<<sourcepop5.long_max1 <<", "<<sourcepop5.long_min2 <<" < l < "<<sourcepop5.long_max2<<" / "
<<sourcepop5. lat_min1 <<" < b < "<<sourcepop5.lat_max1 <<", "<< sourcepop5. lat_min2 <<" < b < "<<sourcepop5. lat_max2<<endl; // AWS20111011
sourcepop5.print(); sourcepop6.print();// AWS20111011
logL=0.0;
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++) // AWS20111011
{
if(pow(10, sourcepop5.lnS_min+ii*sourcepop5.dlnS)>=sourcepop5.flux_detection_limit
&&sourcepop5.dlnN_dlnS[ii]>0.0) // AWS20111011
logL+= - sourcepop5.dlnN_dlnS [ii] + sourcepop6.dlnN_dlnS [ii] * log(sourcepop5.dlnN_dlnS[ii]);// AWS20111011
}
logL_sensitivity_file=0.0;
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++) // AWS20200712
{
if(sourcepop5.dlnN_dlnS_soplimit[ii]>0.0)
logL_sensitivity_file+= - sourcepop5.dlnN_dlnS_soplimit [ii] + sourcepop6.dlnN_dlnS [ii] * log(sourcepop5.dlnN_dlnS_soplimit[ii]);
}
txt_stream << "sourcepop5.n_dlnN_dlnS="<< sourcepop5.n_dlnN_dlnS<<endl;//AWS20111011
txt_stream << "sourcepop6.n_dlnN_dlnS="<< sourcepop6.n_dlnN_dlnS<<endl;//AWS20111011
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++)
{
if( sourcepop5.dlnN_dlnS[ii]>0 || sourcepop6.dlnN_dlnS[ii]>0 )
txt_stream <<"ii="<<ii
<<" lnS="<< sourcepop5.lnS_min+ii*sourcepop5.dlnS <<"/"<< sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS
<<" S=" <<pow(10, sourcepop5.lnS_min+ii*sourcepop5.dlnS) <<"/"<<pow(10, sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS)
<<" N(S)= "<< sourcepop5.dlnN_dlnS [ii]<<":" << sourcepop6.dlnN_dlnS [ii]
<<" N(>S)= "<< sourcepop5.dlnN_dlnS_int[ii]<<":" << sourcepop6.dlnN_dlnS_int[ii]
<<" F(S)= "<< sourcepop5.FS [ii]<<":" << sourcepop6.FS [ii]
<<" F(>S)= "<< sourcepop5.FS_int [ii]<<":" << sourcepop6.FS_int [ii]
<<endl;
}
txt_stream <<endl;
txt_stream << "sourcepop5: model sources above threshold, sourcepop6: Fermi catalogue 2"<<endl;
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++)
{
if( sourcepop5.dlnN_dlnS[ii]>0 || sourcepop6.dlnN_dlnS[ii]>0 )
txt_stream <<"ii="<<ii
<<" lnS="<< sourcepop5.lnS_min+ii*sourcepop5.dlnS <<"/"<< sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS
<<" S=" <<pow(10, sourcepop5.lnS_min+ii*sourcepop5.dlnS) <<"/"<<pow(10, sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS)
<<" N(S)= "<< sourcepop5.dlnN_dlnS_soplimit [ii]<<":" << sourcepop6.dlnN_dlnS [ii]
<<" N(>S)= "<< sourcepop5.dlnN_dlnS_int_soplimit[ii]<<":" << sourcepop6.dlnN_dlnS_int[ii]
<<" F(S)= "<< sourcepop5.FS_soplimit [ii]<<":" << sourcepop6.FS [ii]
<<" F(>S)= "<< sourcepop5.FS_int_soplimit [ii]<<":" << sourcepop6.FS_int [ii]
<<endl;
}
txt_stream <<endl;
txt_stream << "sourcepop5: model sources below threshold, sourcepop6: Fermi catalogue 2"<<endl;
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++)
{
if( sourcepop5.dlnN_dlnS[ii]>0 || sourcepop6.dlnN_dlnS[ii]>0 )
txt_stream <<"ii="<<ii
<<" lnS="<< sourcepop5.lnS_min+ii*sourcepop5.dlnS <<"/"<< sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS
<<" S=" <<pow(10, sourcepop5.lnS_min+ii*sourcepop5.dlnS) <<"/"<<pow(10, sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS)
<<" N(S)= "<< sourcepop5.dlnN_dlnS_sublimit [ii]<<":" << sourcepop6.dlnN_dlnS [ii]
<<" N(>S)= "<< sourcepop5.dlnN_dlnS_int_sublimit[ii]<<":" << sourcepop6.dlnN_dlnS_int[ii]
<<" F(S)= "<< sourcepop5.FS_sublimit [ii]<<":" << sourcepop6.FS [ii]
<<" F(>S)= "<< sourcepop5.FS_int_sublimit [ii]<<":" << sourcepop6.FS_int [ii]
<<endl;
}
// one-line summary statistics
txt_stream << " galdef_ID = " << galplotdef.galdef_ID;
txt_stream << " psfile_tag = " << galplotdef.psfile_tag;
txt_stream << " sourcepop6 vs sourcepop5: cat 2 versus model: ";
txt_stream << " density = " << sourcepop5.density0;
txt_stream << " L_min = " << sourcepop5.L_min;
txt_stream << " L_max = " << sourcepop5.L_max;
txt_stream << " alpha_L = " << sourcepop5.alpha_L;
txt_stream << " long_min1=" << galplotdef.long_min1;
txt_stream << " long_max1=" << galplotdef.long_max1;
txt_stream << " long_min2=" << galplotdef.long_min2;
txt_stream << " long_max2=" << galplotdef.long_max2;
txt_stream << " lat_min1=" << galplotdef.lat_min1;
txt_stream << " lat_max1=" << galplotdef.lat_max1;
txt_stream << " lat_min2=" << galplotdef.lat_min2;
txt_stream << " lat_max2=" << galplotdef.lat_max2;
txt_stream <<" log likelihood using user-defined sensitivity = "<< logL;
txt_stream <<" log likelihood using sensitivity file = " << logL_sensitivity_file<<endl;
// short form
txt_stream << " " << sourcepop5.density0;
txt_stream << " " << sourcepop5.L_min;
txt_stream << " " << sourcepop5.L_max;
txt_stream << " " << sourcepop5.alpha_L;
txt_stream << " " << galplotdef.long_min1;
txt_stream << " " << galplotdef.long_max1;
txt_stream << " " << galplotdef.long_min2;
txt_stream << " " << galplotdef.long_max2;
txt_stream << " " << galplotdef.lat_min1;
txt_stream << " " << galplotdef.lat_max1;
txt_stream << " " << galplotdef.lat_min2;
txt_stream << " " << galplotdef.lat_max2;
txt_stream <<" " << logL;
txt_stream <<" " << logL_sensitivity_file;
txt_stream <<" logL statistics" <<endl;
//---------------------------------------------------------
txt_stream<<endl<<"Comparison source synthesis : Fermi catalogue 3"<<endl;
txt_stream<<"source counts for region "<<sourcepop5.long_min1 <<" < l < "<<sourcepop5.long_max1 <<", "<<sourcepop5.long_min2 <<" < l < "<<sourcepop5.long_max2<<" / "
<<sourcepop5. lat_min1 <<" < b < "<<sourcepop5.lat_max1 <<", "<< sourcepop5. lat_min2 <<" < b < "<<sourcepop5. lat_max2<<endl; // AWS20111011
sourcepop5.print(); sourcepop7.print();// AWS20111011
logL=0.0;
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++) // AWS20111011
{
if(pow(10, sourcepop5.lnS_min+ii*sourcepop5.dlnS)>=sourcepop5.flux_detection_limit
&&sourcepop5.dlnN_dlnS[ii]>0.0) // AWS20111011
logL+= - sourcepop5.dlnN_dlnS [ii] + sourcepop7.dlnN_dlnS [ii] * log(sourcepop5.dlnN_dlnS[ii]);// AWS20111011
}
logL_sensitivity_file=0.0;
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++) // AWS20200712
{
if(sourcepop5.dlnN_dlnS_soplimit[ii]>0.0)
logL_sensitivity_file+= - sourcepop5.dlnN_dlnS_soplimit [ii] + sourcepop7.dlnN_dlnS [ii] * log(sourcepop5.dlnN_dlnS_soplimit[ii]);
}
txt_stream << "sourcepop5.n_dlnN_dlnS="<< sourcepop5.n_dlnN_dlnS<<endl;//AWS20111011
txt_stream << "sourcepop7.n_dlnN_dlnS="<< sourcepop7.n_dlnN_dlnS<<endl;//AWS20111011
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++)
{
if( sourcepop5.dlnN_dlnS[ii]>0 || sourcepop7.dlnN_dlnS[ii]>0 )
txt_stream <<"ii="<<ii
<<" lnS="<< sourcepop5.lnS_min+ii*sourcepop5.dlnS <<"/"<< sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS
<<" S=" <<pow(10, sourcepop5.lnS_min+ii*sourcepop5.dlnS) <<"/"<<pow(10, sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS)
<<" N(S)= "<< sourcepop5.dlnN_dlnS [ii]<<":" << sourcepop7.dlnN_dlnS [ii]
<<" N(>S)= "<< sourcepop5.dlnN_dlnS_int[ii]<<":" << sourcepop7.dlnN_dlnS_int[ii]
<<" F(S)= "<< sourcepop5.FS [ii]<<":" << sourcepop7.FS [ii]
<<" F(>S)= "<< sourcepop5.FS_int [ii]<<":" << sourcepop7.FS_int [ii]
<<endl;
}
txt_stream <<endl;
txt_stream << "sourcepop5: model sources above threshold, sourcepop7: Fermi catalogue 3"<<endl;
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++)
{
if( sourcepop5.dlnN_dlnS[ii]>0 || sourcepop7.dlnN_dlnS[ii]>0 )
txt_stream <<"ii="<<ii
<<" lnS="<< sourcepop5.lnS_min+ii*sourcepop5.dlnS <<"/"<< sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS
<<" S=" <<pow(10, sourcepop5.lnS_min+ii*sourcepop5.dlnS) <<"/"<<pow(10, sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS)
<<" N(S)= "<< sourcepop5.dlnN_dlnS_soplimit [ii]<<":" << sourcepop7.dlnN_dlnS [ii]
<<" N(>S)= "<< sourcepop5.dlnN_dlnS_int_soplimit[ii]<<":" << sourcepop7.dlnN_dlnS_int[ii]
<<" F(S)= "<< sourcepop5.FS_soplimit [ii]<<":" << sourcepop7.FS [ii]
<<" F(>S)= "<< sourcepop5.FS_int_soplimit [ii]<<":" << sourcepop7.FS_int [ii]
<<endl;
}
txt_stream <<endl;
txt_stream << "sourcepop5: model sources below threshold, sourcepop7: Fermi catalogue 3"<<endl;
for (ii=0;ii< sourcepop5.n_dlnN_dlnS;ii++)
{
if( sourcepop5.dlnN_dlnS[ii]>0 || sourcepop7.dlnN_dlnS[ii]>0 )
txt_stream <<"ii="<<ii
<<" lnS="<< sourcepop5.lnS_min+ii*sourcepop5.dlnS <<"/"<< sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS
<<" S=" <<pow(10, sourcepop5.lnS_min+ii*sourcepop5.dlnS) <<"/"<<pow(10, sourcepop5.lnS_min+(ii+1)*sourcepop5.dlnS)
<<" N(S)= "<< sourcepop5.dlnN_dlnS_sublimit [ii]<<":" << sourcepop7.dlnN_dlnS [ii]
<<" N(>S)= "<< sourcepop5.dlnN_dlnS_int_sublimit[ii]<<":" << sourcepop7.dlnN_dlnS_int[ii]
<<" F(S)= "<< sourcepop5.FS_sublimit [ii]<<":" << sourcepop7.FS [ii]
<<" F(>S)= "<< sourcepop5.FS_int_sublimit [ii]<<":" << sourcepop7.FS_int [ii]
<<endl;
}
// one-line summary statistics
txt_stream << " galdef_ID = " << galplotdef.galdef_ID;
txt_stream << " psfile_tag = " << galplotdef.psfile_tag;
txt_stream << " sourcepop7 vs sourcepop5: cat 3 versus model: ";
txt_stream << " density = " << sourcepop5.density0;
txt_stream << " L_min = " << sourcepop5.L_min;
txt_stream << " L_max = " << sourcepop5.L_max;
txt_stream << " alpha_L = " << sourcepop5.alpha_L;
txt_stream << " long_min1=" << galplotdef.long_min1;
txt_stream << " long_max1=" << galplotdef.long_max1;
txt_stream << " long_min2=" << galplotdef.long_min2;
txt_stream << " long_max2=" << galplotdef.long_max2;
txt_stream << " lat_min1=" << galplotdef.lat_min1;
txt_stream << " lat_max1=" << galplotdef.lat_max1;
txt_stream << " lat_min2=" << galplotdef.lat_min2;
txt_stream << " lat_max2=" << galplotdef.lat_max2;
txt_stream <<" log likelihood using user-defined sensitivity = " << logL;
txt_stream <<" log likelihood using sensitivity file = " << logL_sensitivity_file<<endl;
// short form
txt_stream << " " << sourcepop5.density0;
txt_stream << " " << sourcepop5.L_min;
txt_stream << " " << sourcepop5.L_max;
txt_stream << " " << sourcepop5.alpha_L;
txt_stream << " " << galplotdef.long_min1;
txt_stream << " " << galplotdef.long_max1;
txt_stream << " " << galplotdef.long_min2;
txt_stream << " " << galplotdef.long_max2;
txt_stream << " " << galplotdef.lat_min1;
txt_stream << " " << galplotdef.lat_max1;
txt_stream << " " << galplotdef.lat_min2;
txt_stream << " " << galplotdef.lat_max2;
txt_stream <<" " << logL;
txt_stream <<" " << logL_sensitivity_file;
txt_stream <<" logL statistics" <<endl;
//---------------------------------------------------------
/*AWS20110818
txt_stream<<endl<<"Comparison source synthesis : SPI catalogue"<<endl;
txt_stream<<"source counts for region "<<sourcepop1.long_min1 <<" < l < "<<sourcepop1.long_max1 <<", "<<sourcepop1.long_min2 <<" < l < "<<sourcepop1.long_max2<<" / "
<<sourcepop1. lat_min1 <<" < b < "<<sourcepop1.lat_max1 <<", "<< sourcepop1. lat_min2 <<" < b < "<<sourcepop1. lat_max2<<endl;
// sourcepop1.print(); sourcepop3.print();
logL=0.0;
for (ii=0;ii< sourcepop1.n_dlnN_dlnS;ii++)
{
if(pow(10, sourcepop1.lnS_min+ii*sourcepop1.dlnS)>=sourcepop1.flux_detection_limit
&&sourcepop1.dlnN_dlnS[ii]>0.0)
logL+= - sourcepop1.dlnN_dlnS [ii] + sourcepop3.dlnN_dlnS [ii] * log(sourcepop1.dlnN_dlnS[ii]);
}
for (ii=0;ii< sourcepop1.n_dlnN_dlnS;ii++)
{
if( sourcepop1.dlnN_dlnS[ii]>0 || sourcepop3.dlnN_dlnS[ii]>0 )
txt_stream <<"ii="<<ii
<<" lnS="<< sourcepop1.lnS_min+ii*sourcepop1.dlnS <<"/"<< sourcepop1.lnS_min+(ii+1)*sourcepop1.dlnS
<<" S=" <<pow(10, sourcepop1.lnS_min+ii*sourcepop1.dlnS) <<"/"<<pow(10, sourcepop1.lnS_min+(ii+1)*sourcepop1.dlnS)
<<" N(S)= "<< sourcepop1.dlnN_dlnS [ii]<<":" << sourcepop3.dlnN_dlnS [ii]
<<" N(>S)= "<< sourcepop1.dlnN_dlnS_int[ii]<<":" << sourcepop3.dlnN_dlnS_int[ii]
<<" F(S)= "<< sourcepop1.FS [ii]<<":" << sourcepop3.FS [ii]
<<" F(>S)= "<< sourcepop1.FS_int [ii]<<":" << sourcepop3.FS_int [ii]
<<endl;
}
txt_stream <<"log likelihood "<< logL<<endl;
*/
cout<<" <<<< gen_source_population "<<endl;
return status;
}