Skip to content
Snippets Groups Projects
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;
}