Changeset 877 in Sophya for trunk/SophyaProg/PMixer/extractRS.cc

Timestamp:

Apr 11, 2000, 3:24:27 PM (25 years ago)

Author:

ansari

Message:

Sophie cf. skymixer : idem::

File:

• trunk/SophyaProg/PMixer/extractRS.cc (modified) (18 diffs)

trunk/SophyaProg/PMixer/extractRS.cc

@@ -12 +12 @@
 #include "tvector.h"
 #include "vector3d.h"
-
+#include "nbrandom.h"
 #include "timing.h"
 #include "sambainit.h"
@@ -34 +34 @@
 SpectralResponse * getSpectralResponse(DataCards & dc);
 template <class T>
-void addComponent(SpectralResponse&  sr, PixelMap<T>& finalMap, 
-                  PixelMap<T>& mapToAdd, RadSpectra& rs, double K=1.);
-template <class T>
 void  MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig);
 float ComputeFrequency(DataCards &,string&);
@@ -45 +42 @@
 static bool rdmap = false;    // true -> Read map first 
 static char mapPath[256];     // Path for input maps
+static int  hp_nside = 32;    // HealPix NSide
 static int  nskycomp = 0;     // Number of sky components
 static int  debuglev = 0;     // Debug Level
@@ -55 +53 @@
 int main(int narg, char * arg[])
 {
-  if ((narg < 4) || ((narg > 1) && (strcmp(arg[1], "-h") == 0) )) {
-    cout << "  Usage: extractRS parameterFile outputfitsnameformaps  outputfitsnameforsource [outppfname]" << endl;
+  if ((narg < 3) || ((narg > 1) && (strcmp(arg[1], "-h") == 0) )) {
+    cout << "  Usage: extractRS parameterFile outputfitsnameformaps  [outppfname]" << endl;
     exit(0);
   }
@@ -90 +88 @@
   
   
-  cout << " extractRS/Info : NComp = " <<  nskycomp << endl;
+  cout << " extractRS/Info : NComp = " <<  nskycomp << " SphereGorski_NSide= " << hp_nside << endl;
   cout << "  ... MapPath = " << (string)mapPath << "  DebugLev= " << debuglev 
        << "  PrintLev= " << printlev << endl;
@@ -97 +95 @@
   // We create an output persist file for writing debug objects
   if (debuglev > 0) so = new POutPersist("extrRSdbg.ppf");
-  
-  
-  TMatrix<float> outgs;
-  bool okout = false;
+  SphereGorski<float> SourceMap(hp_nside);
+  SphereGorski<float> OutputMap(hp_nside);
+  bool okout = true;
   
   try {
@@ -120 +117 @@
     int nb_source;
     TMatrix<float> rareSource; //par freq et par numero de source
-        
+    
     for(sk = 0; sk<maxOfSkyComp; sk++) 
       {
-        
         TMatrix<float> ings;    
         // Opening the file containing the map
@@ -137 +133 @@
         if (printlev > 2) 
           {
-            MeanSig(ings.DataBlock(), moy, sig );
-            cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
-          }
-
+            double min = 1000000000;
+            double max = -10000000000;
+            for(int x=0;x<ings.NRows(); x++)
+              {
+                for(int y=0;y<ings.NCols(); y++)
+                  {
+                    if(min>ings(x,y)) min = ings(x,y);
+                    if(max<ings(x,y)) max = ings(x,y);
+                  }
+              }
+            cout << "min and max values for input map" << min << " " << max << endl;
+          }
         // Computing the frequency according to the name of the file
         freq(sk) =  ComputeFrequency(dc,key);
         if (printlev > 2)  cout << "filling spectrum for freq = " << freq(sk) << endl;
+        
         // Opening the file containing the point source
         sprintf(buff, "%d", sk+1);
@@ -169 +174 @@
         for(int x=0;x<ings.NRows(); x++)
           {
-            for(int y=0;x<ings.NCols(); x++)
+            for(int y=0;y<ings.NCols(); y++)
               {
                 spectrum(x,y) = ings(x,y)*1.e-26; //for Jansky->W
@@ -185 +190 @@
         K = dc.DParam(key, 1, 1.);
         if (debuglev > 4) {  // Writing the input map to the outppf
-          so->PutObject(ings, key);    // $CHECK$  Reza 05/04/2000 
+          FIO_TMatrix<float> fiog;
+          fiog.Write(*so, key);
         }
         
@@ -205 +211 @@
     TVector<float> outCoeff(nbOfPix);
 
-    for(int pix=1;pix<nbOfPix;pix++)
+    for(int pix=0;pix<nbOfPix;pix++)
       {
         Vector smallFDeNu(maxOfSkyComp);
@@ -216 +222 @@
         if (printlev > 10)  cout << "RadiationSpectra" << radSV << endl;
 
-        TVector<float> smallFDeNuT(maxOfSkyComp);
-        for(sk = 0; sk<maxOfSkyComp; sk++) 
-          {
-            smallFDeNuT(sk) = smallFDeNu(sk);
-          }
         outCoeff(pix) = radSV.filteredIntegratedFlux(*sr);
       }
-
-    FitsIoServer fiosrs;
-    fiosrs.save(outCoeff,arg[2]);
-    cout << "Output Map (TMatrix<float>) written to FITS file "  
-         << arg[2]<< endl;
-    cout << endl;
-
+    
+    for(int i=0;i<OutputMap.NbPixels();i++)
+      {
+        float random = frand01();
+        if(random == 1) random = frand01();
+        int randomPix = random*OutputMap.NbPixels();
+        int outpix = random*nbOfPix; 
+        OutputMap(i) = outCoeff(outpix)*1.e18;
+      }
+
+    
     // ---------------------------------------------------------------------
     //      Integration in the detector band of the sources
@@ -244 +249 @@
         sourceCoeff(nsource) = radSV_Source.filteredIntegratedFlux(*sr);
       }
-    FitsIoServer fiosrs_source;
-    fiosrs_source.save(sourceCoeff,arg[3]);
-    cout << "Output Map (TMatrix<float>) written to FITS file "  
-         << arg[3]<< endl;
-    cout << endl;
-
-    // ---------------------------------------------------------------------
-    //     Check of write/read on FITS file
-    // ---------------------------------------------------------------------
-    if (printlev > 2)
-    {
-      cout << " Check of write/read on FITS file for skymap!" << endl;
-      TVector<float> outout;
-      FitsIoServer fiotest;
-      fiotest.load(outout, arg[2]);
-      cout << "Coeff(2)" << outout(2) << " ? " << outCoeff(2) << endl;
-      cout << endl;
-    }
-    if (printlev > 2)
-    {
-      cout << " Check of write/read on FITS file for sources!" << endl;
-      FitsIoServer fiotest;
-      TVector<float> outout2;
-      fiotest.load(outout2, arg[3]);
-      cout << "sourceCoeff(2)" << outout2(2) << " ? " << sourceCoeff(2) << endl;
-    }
-
-
-  }   // End of try block
-  
-  
+
+    for(int i=0;i<nb_source;i++)
+      {
+        float random = frand01();
+        int outPix = random*OutputMap.NbPixels(); 
+        SourceMap(outPix) = sourceCoeff(i); 
+        OutputMap(outPix) = OutputMap(i)+ SourceMap(outPix)*1.E18;
+      }
+  }  
   catch (PException exc) {
     msg = exc.Msg();
@@ -291 +274 @@
   
   // Saving the output map in FITS format 
-  if (okout) {
-    FitsIoServer fios;
-    fios.save(outgs, arg[2]);
-    cout << "Output Map (TMatrix<float>) written to FITS file " 
-         << (string)(arg[2]) << endl;
-    if(printlev>2) PrtTim("End of WriteFITS ");
-    // Saving the output map in PPF format 
-    if (narg > 3) {
-      POutPersist s(arg[3]); 
-      s.PutObject(outgs);   //  $CHECK$  Reza 05/04/2000
-      cout << "Output Map (TMatrix<float>) written to POutPersist file "  
-           << (string)(arg[3]) << endl;
-      if(printlev>2) PrtTim("End of WritePPF ");
+  if (okout) 
+    {
+      //       {
+      //        FitsIoServer fios;
+      //        fios.save(SourceMap, arg[3]);
+      //        double moy,sig;
+      //        MeanSig(SourceMap.DataBlock(),moy,sig);
+      //        cout << " MeanSig for Source Map - Mean= " << moy << " Sigma= " << sig << endl;
+      //        cout << "Source Map (SphereGorski<float>) written to FITS file " 
+      //             << (string)(arg[3]) << endl;
+      //       }
+      {
+        FitsIoServer fios;
+        fios.save(OutputMap, arg[2]);
+        double moy,sig;
+        MeanSig(OutputMap.DataBlock(),moy,sig);
+        cout << " MeanSig for Output Map - Mean= " << moy << " Sigma= " << sig << endl;
+        cout << "Output Map (SphereGorski<float>) written to FITS file " 
+             << (string)(arg[2]) << endl;
+      }
+      if(printlev>2) PrtTim("End of WriteFITS ");
+      // Saving the output map in PPF format 
+      if (narg > 2) {
+        POutPersist s(arg[3]); 
+        FIO_SphereGorski<float> fiog(OutputMap);
+        fiog.Write(s);
+        cout << "Output Map (SphereGorski<float>) written to POutPersist file "  
+             << (string)(arg[3]) << endl;
+        if(printlev>2) PrtTim("End of WritePPF ");
+      }
     }
-  }
   if (so) delete so;  //  Closing the debug ppf file 
   return(rc);
@@ -326 +325 @@
   
   //  Checking datacards
-  if (dc.NbParam("EXT_RS") < 1)   {
+  if (dc.NbParam("EXT_RS") < 2)   {
      rc = 71; 
      msg = "Invalid parameters  - Check @EXT_RS card ";
@@ -337 +336 @@
   int ncomp = 0;
   ncomp = dc.IParam("EXT_RS", 0);
+  int nside = 32;
+  nside = dc.IParam("EXT_RS", 1);
   
   for(kc=0; kc<ncomp; kc++) 
@@ -371 +372 @@
   msg = "OK";
   nskycomp = ncomp;
-  
+  hp_nside = nside;
   // Checking for PATH definition card
   key = "RS_EXT_PATH";
@@ -409 +410 @@
   if (gsig >= 0.) gsig = sqrt(gsig);
 }
-/* Nouvelle-Fonction */
-template <class T> 
-void addComponent(SpectralResponse&  sr, PixelMap<T>& finalMap, 
-                  PixelMap<T>& mapToAdd, RadSpectra& rs, double K)
-{
-  // finalMap = finalMap + coeff* mapToAdd
-  // coeff    =  convolution of sr and rs 
-  // compute the coefficient corresponding to mapToAdd
-  if (finalMap.NbPixels() != mapToAdd.NbPixels())   
-    throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
-  double coeff = rs.filteredIntegratedFlux(sr) * K;
-  if (printlev > 1) 
-    cout << " addComponent - Coeff= " << coeff << " (K= " << K << ")" << endl; 
-  for(int i=0; i<finalMap.NbPixels(); i++)
-    {
-      finalMap(i) += coeff * mapToAdd(i);
-    }
-}
-
 
 /* Nouvelle-Fonction */