Minor bug: modifiy to use more than 128 bands. Use only one value for...

Minor bug: modifiy to use more than 128 bands. Use only one value for ERR_SCALE. Be careful at the extreme z limits for limits

sed/QSO/SALVATO2015/best_MAR2010.list~

-MARA23032010/CB1_0_LOIII4.sed A 
-MARA23032010/S0_template_norm.sed A
-MARA23032010/Sb_template_norm.sed A
-MARA23032010/Spi4_template_norm.sed A
-MARA23032010/M82_template_norm.sed A 
-MARA23032010/I22491_template_norm.sed.save A 
-MARA23032010/Sey18_template_norm.sed A
-MARA23032010/Sey2_template_norm.sed A
-MARA23032010/S0_10_QSO2_90.sed A
-MARA23032010/S0_20_QSO2_80.sed A
-MARA23032010/S0_30_QSO2_70.sed A
-MARA23032010/S0_40_QSO2_60.sed A
-MARA23032010/S0_50_QSO2_50.sed A
-MARA23032010/S0_60_QSO2_40.sed A
-MARA23032010/S0_70_QSO2_30.sed A
-MARA23032010/S0_80_QSO2_20.sed A
-MARA23032010/S0_90_QSO2_10.sed A
-MARA23032010/Mrk231_template_norm.sed A
-MARA23032010/I22491_90_TQSO1_10.sed  A
-MARA23032010/I22491_80_TQSO1_20.sed A
-MARA23032010/I22491_70_TQSO1_30.sed A
-MARA23032010/I22491_60_TQSO1_40.sed A
-MARA23032010/I22491_50_TQSO1_50.sed A
-MARA23032010/I22491_40_TQSO1_60.sed A
-MARA23032010/pl_I22491_30_TQSO1_70.sed A
-MARA23032010/pl_I22491_20_TQSO1_80.sed A
-MARA23032010/pl_I22491_10_TQSO1_90.sed A
-MARA23032010/pl_QSOH_template_norm.sed A
-MARA23032010/pl_QSO_DR2_029_t0.spec A
-MARA23032010/pl_TQSO1_template_norm.sed A 

source/SED.cpp

@@ -1246,24 +1246,18 @@ void GalSED::generateEmPhys(double zmet, double qi) {
 
   // For Z < 0.03 in solar unit, Z1=0.0004
   double Z1_line[65] = { 22.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.245, 0.245, 0.295, 0.203, 0.27, 0.015, 0.005, 0.002, 0.256, 0.109, 0.466, 0.036, 0.01, 1.0, 1.097, 3.159, 0.003, 0.0, 0.096, 0.008, 0.009, 0.005, 2.87, 0.015, 0.026, 0.037, 0.029, 0.028, 0.027, 0.012, 0.007, 0.067, 0.0, 0.0, 0.0, 0.0, 0.0, 0.232, 0.352, 0.165, 0.0906, 0.0368, 0.0366, 0.0254, 0.0184, 0.0138, 0.0106, 0.0185, 0.00666, 0.00541, 0.0884, 0.0471, 0.0281, 0.0186, 0.0127, 0.00914, 0.00683, 0.00524 };
-  // Ratio from Patrice Cloudy
-  // double Z1_line[65] = {25.498, 0.0, 0.125, 0.257, 0.003, 0.037, 0.377, 0.416, 0.568, 0.505, 0.209, 0.317, 4026.191, 0.007, 0.002, 0.265, 0.082, 0.478, 0.041, 0.01, 1.0, 2.521, 7.521, 0.001, 0.001, 0.104, 0.019, 0.019, 0.02, 2.779, 0.058, 0.029, 0.076, 0.058, 0.026, 0.093, 0.013, 0.011, 0.022, 0.19, 0.478, 0.002, 0.002, 0.002, 0.222, 0.296, 0.152, 0.086, 0.053, 0.036, 0.025, 0.018, 0.0138, 0.0106, 0.0185, 0.00666, 0.00541, 0.065, 0.04, 0.025, 0.017, 0.012, 0.009, 0.00683, 0.00524 };
   // check
   longEL = sizeof(Z1_line)/sizeof(Z1_line[0]);
   if(longEL!= 65)cout << " Error with size of lines Z1 " << longEL << endl;
 
    // 0.03 < Z <= 0.3 in solar unit, Z2=0.004
    double Z2_line[65] = { 22.2, 0.0, 0.058, 0.0, 0.0, 0.0, 0.31, 0.896, 0.896, 0.416, 0.192, 0.283, 0.015, 0.017, 0.007, 0.256, 0.066, 0.466, 0.036, 0.014, 1.0, 1.617, 4.752, 0.01, 0.0, 0.108, 0.041, 0.017, 0.059, 2.87, 0.175, 0.03, 0.188, 0.138, 0.023, 0.071, 0.027, 0.014, 0.176, 0.2, 0.51, 0.0, 0.0, 0.0, 0.225, 0.352, 0.165, 0.0906, 0.0368, 0.0366, 0.0254, 0.0184, 0.0138, 0.0106, 0.0185, 0.00666, 0.00541, 0.0884, 0.0471, 0.0281, 0.0186, 0.0127, 0.00914, 0.00683, 0.00524 };
-  // Ratio from Patrice Cloudy
-  // double Z2_line[65] = { 23.988, 0.0, 0.018, 0.116, 0.004, 0.073, 0.612, 0.753, 1.023, 0.315, 0.211, 0.26, 0.015, 0.013, 0.004, 0.265, 0.022, 0.477, 0.043, 0.014, 1.0, 1.631, 4.866, 0.005, 0.003, 0.111, 0.036, 0.015, 0.081, 2.801, 0.238, 0.031, 0.196, 0.149, 0.023, 0.117, 0.016, 0.014, 0.028, 0.279, 0.701, 0.003, 0.004, 0.003, 0.218, 0.312, 0.159, 0.089, 0.055, 0.037, 0.026, 0.019, 0.0138, 0.0106, 0.0185, 0.00666, 0.00541, 0.07, 0.042, 0.027, 0.018, 0.012, 0.009, 0.00683, 0.00524};
   // Check
   longEL = sizeof(Z2_line)/sizeof(Z2_line[0]);
   if(longEL!= 65)cout << " Error with size of lines Z2 " << longEL << endl;
 
    // 0.3 < Z in solar unit, Z3-Z5=0.008-0.05
    double Z3_line[65] = { 22.2, 0.11, 0.01, 0.18, 0.01, 0.29, 0.07, 1.505, 1.505, 0.3, 0.107, 0.159, 0.015, 0.029, 0.011, 0.256, 0.01, 0.466, 0.05, 0.0, 1.0, 1.399, 4.081, 0.03, 0.01, 0.14, 0.13, 0.03, 0.136, 2.87, 0.404, 0.03, 0.3, 0.21, 0.04, 0.035, 0.026, 0.014, 0.086, 0.365, 0.945, 0.048, 0.058, 0.054, 0.229, 0.352, 0.165, 0.0906, 0.0368, 0.0366, 0.0254, 0.0184, 0.0138, 0.0106, 0.0185, 0.00666, 0.00541, 0.0884, 0.0471, 0.0281, 0.0186, 0.0127, 0.00914, 0.00683, 0.00524  };
-  // Ratio from Patrice Cloudy
-  // double Z3_line[65] = { 23.637, 0.11, 0.0, 0.004, 0.0, 0.012, 0.113, 0.277, 0.369, 0.043, 0.214, 0.175, 0.015, 0.006, 0.002, 0.263, 0.001, 0.474, 0.048, 0.0, 1.0, 0.289, 0.863, 0.004, 0.001, 0.125, 0.009, 0.003, 0.139, 2.856, 0.409, 0.036, 0.203, 0.158, 0.022, 0.065, 0.002, 0.001, 0.015, 0.196, 0.493, 0.001, 0.002, 0.001, 0.22, 0.337, 0.168, 0.094, 0.058, 0.038, 0.027, 0.019, 0.0138, 0.0106, 0.0185, 0.00666, 0.00541, 0.078, 0.047, 0.029, 0.019, 0.013, 0.01, 0.00683, 0.00524};  
   // Check
   longEL = sizeof(Z3_line)/sizeof(Z3_line[0]);
   if(longEL!= 65)cout << " Error with size of lines Z3 " << longEL << endl;
@@ -2067,7 +2061,6 @@ void GalSED::fillPhys(const vector<double>& phys)
  */
 void GalSED::fit(const onesource& source, const int imagm, const string zfix,  const double consideredZ, const double lmasi, const double lmass,const double magabsB[2], const double magabsF[2], const double funz0, const int bp[2]){
 
-
        // Initialize the chi2 by a crazy high value in case there is no fit
        chi2=1e9;
        // Case ZFIX=YES with the spec-z already set, reject other solutions in the library
@@ -2080,7 +2073,7 @@ void GalSED::fit(const onesource& source, const int imagm, const string zfix,  c
          // Upper-limits
          for (int k=0;k<imagm;k++){
            if (source.bused[k]==1 && source.sab[k]<=0 && (dm*(flux[k]))>source.ab[k]){
-                  reject=1;
+	         reject=1;
                }
          }
          // Use mass scaling for model rejection

source/globals.h

@@ -4,7 +4,7 @@
 #include <string>
 
 #include <bitset>
-#define MAX_CONTEXT 128
+#define MAX_CONTEXT 1024
 #define CHECK_CONTEXT_BIT(context, n) (std::bitset<MAX_CONTEXT>(context).test(n))
 
 using namespace std;

source/onesource.cpp

@@ -69,7 +69,6 @@ void onesource::readsource(const string line,const int cat_fmt,const string catt
      getline( ss, str_inp); // store the end of the file into a string
   }
 
-
   cout << "Read source " << spec << " context " << cont << " spectro-z " << zs << " \r " << flush;
 
   return;
@@ -322,11 +321,15 @@ void onesource::substellar(const string substar, vector<flt> allFilters){
 void onesource::errors(const vector<double> min_err, const vector<double> fac_err,const int imagm){
 
   // check if the number of filters in min_err corresponds to the expected one
-  if( (size_t)imagm== min_err.size()){
+  if( (size_t)imagm== min_err.size() || min_err.size()==1){
     for(int k=0; k<imagm; k++){ // loop over the filters
      if(sab[k]>0){
       sab[k]=1.086*sab[k]/ab[k];
-      sab[k]=sqrt(pow(sab[k],2.)+pow(min_err[k],2.)); // Add in quadrature
+      if(min_err.size()>1){
+	sab[k]=sqrt(pow(sab[k],2.)+pow(min_err[k],2.)); // Add in quadrature
+      }else{
+	sab[k]=sqrt(pow(sab[k],2.)+pow(min_err[0],2.)); // Add in quadrature
+      }	
       // Put abs(ab[k]) if negative flux
       sab[k]=abs(ab[k])*sab[k]/1.086;
      }
@@ -337,9 +340,13 @@ void onesource::errors(const vector<double> min_err, const vector<double> fac_er
 
   // multiply all the errors in flux by a common factor
   // check if the number of filters in fac_err corresponds to the expected one
-  if((size_t)imagm==fac_err.size()){
+  if((size_t)imagm==fac_err.size() || fac_err.size()==1){
     for(int k=0; k<imagm; k++){
-      if(sab[k]>0) sab[k]=sab[k]*fac_err[k];}
+      if(fac_err.size()>1){
+         if(sab[k]>0) sab[k]=sab[k]*fac_err[k];
+       }else{
+         if(sab[k]>0) sab[k]=sab[k]*fac_err[0];}
+    }
   }else{
     cout << " Can not add error in quadrature: dimension of fac_err " << imagm << " " << fac_err.size() << endl;
   }
@@ -1394,7 +1401,7 @@ void onesource::limits(vector<SED*> fulllib, vector<double> limits_zbin,int limi
   double diff=0;
 
   // Check if feasible
-  if(indmin[0]>=0){
+  if((indmin[0]>=0) && (size_t)(indmin[0])<fulllib.size()-1){
 
     // Loop over the redshift bin for the selection filter
     FltSel=-1;
@@ -1412,10 +1419,10 @@ void onesource::limits(vector<SED*> fulllib, vector<double> limits_zbin,int limi
      
     // check if the difference between the observed magnitude and magnitude limit 
     while(diff < (MagSel - mab[FltSel]) ){
-
         // check if the difference between the observed magnitude at the minimum redshift and maximum redshift
         diff = fulllib[indmin[0]+deca]->mag_ori[FltSel] - fulllib[indmin[0]]->mag_ori[FltSel];
         deca+=1;
+	
         // Check that we are below the size of the library
         if((size_t)(indmin[0]+deca)>=fulllib.size()-1){deca-=1 ; break;}
         // Check that we didn't reach the maximum redshift of this model

source/photoz_lib.cpp

@@ -86,7 +86,13 @@ PhotoZ::PhotoZ(keymap& key_analysed)
   int nerr = int(min_err.size());
 
   // ERR_FACTOR Multiply the flux uncertainties by a given factor - 1.0 default
-  fac_err = ((key_analysed["ERR_FACTOR"]).split_double("1.0",nerr)); //should be # of filters instead of nerr
+  fac_err = ((key_analysed["ERR_FACTOR"]).split_double("1.0",-1));
+  int nfac = int(fac_err.size());
+  if((nerr>1) && (nfac>1) && (nfac!=nerr)){
+    cout << "The number of filters in ERR_SCALE and ERR_FACTOR do not correspond." << endl;
+    cout << " Error in the keywords. " << endl;
+    cout << " Stop " << endl;
+  }
 
   // BD_SCALE Bands used for the scaling - 0 default (all bands)
   bdscal = ((key_analysed["BD_SCALE"]).split_long("0",1))[0];
@@ -972,7 +978,6 @@ vector<vector<int>> bestFilter(int nbFlt, vector<double> gridz, vector<SED*> ful
   vector<vector<int>> bestFlt;
   vector<int> bestFlt_z;
   int j=0;
-
   // Loop over the grid in redshift
   for(size_t i=0; i<gridz.size();i++){
 
@@ -1102,7 +1107,6 @@ void minimizekcolor(vector<double> gridz, vector<SED*> fulllib, vector<vector<in
   int imagm=bestFlt[0].size();
   int nbz=bestFlt.size();
   
-  
   // initialize the maximum and minimum values
   for(int i=0; i<nbz;i++){ // Loop over all redshifts 
     for (int k=0; k<imagm;k++){ // loop over all filters
@@ -1127,8 +1131,7 @@ void minimizekcolor(vector<double> gridz, vector<SED*> fulllib, vector<vector<in
    }
   }
   int nbPossi=possiFlt.size();
-  
-  
+    
   // loop over the reference filters to be considered
   for (int r=0; r<imagm; ++r){ 
     
@@ -1278,7 +1281,6 @@ void PhotoZ::run_photoz(const GalMag& opaOut, const vector<double>& a0,
   // MIN_THRES threshold to trigger the detection - 0.1 by default
   double min_thres  =((keys["MIN_THRES"]).split_double("0.1",1))[0];
 
-
   /* Define what are the filters to be used for the absolute magnitude depending on the method adopted */
   // MABS_METHOD method to compute the absolute magnitudes
   int method = ((keys["MABS_METHOD"]).split_int("0",-1))[0];
@@ -1344,6 +1346,7 @@ void PhotoZ::run_photoz(const GalMag& opaOut, const vector<double>& a0,
       onesource oneObj(nobj,gridz);
       // Read the object
       oneObj.readsource(line,cat_fmt,cattyp,imagm);
+      cout << "Entree " << oneObj.spec << endl;
       // Use zspec from external file 
       // open the external file with zspec  
       if ( externalzfile.substr(0, 4)!="NONE" ){
@@ -1399,8 +1402,7 @@ void PhotoZ::run_photoz(const GalMag& opaOut, const vector<double>& a0,
       if((methz[0]=='M' || methz[0]=='m') && (zfix[0]=='N' || zfix[0]=='n')){
 	oneObj.chimin[0]=1.e9;
 	oneObj.fit(fullLib,imagm,"YES",oneObj.zml,lmasi,lmass,magabsB,magabsF,funz0,bp);
-      }
-      
+      }      
       // Interpolation at the new redshift  (only gal for the moment)
       oneObj.interp_lib(fullLib,imagm,lcdm);
       // Compute absolute magnitudes
@@ -1426,6 +1428,7 @@ void PhotoZ::run_photoz(const GalMag& opaOut, const vector<double>& a0,
 	// Compute the IR luminosities
 	oneObj.generatePDF_IR(fullLibIR);
       }
+      cout << "Sortie3 " << oneObj.spec << endl;
       
       // write the object in output
       oneObj.write_out(fullLib,fullLibIR,stout,outkeywords);
@@ -1440,6 +1443,7 @@ void PhotoZ::run_photoz(const GalMag& opaOut, const vector<double>& a0,
       // Computes proba in various intervals of z
       oneObj.probaz();
       */
+      cout << "Sortie " << oneObj.spec << endl;
       
     }
   }