Index: trunk/SophyaProg/PMixer/extractRS.cc
===================================================================
--- trunk/SophyaProg/PMixer/extractRS.cc	(revision 817)
+++ trunk/SophyaProg/PMixer/extractRS.cc	(revision 877)
@@ -12,5 +12,5 @@
 #include "tvector.h"
 #include "vector3d.h"
-
+#include "nbrandom.h"
 #include "timing.h"
 #include "sambainit.h"
@@ -34,7 +34,4 @@
 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,4 +42,5 @@
 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,6 +53,6 @@
 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,5 +88,5 @@
   
   
-  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,8 +95,7 @@
   // 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,8 +117,7 @@
     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,11 +133,20 @@
 	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,5 +174,5 @@
 	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,5 +190,6 @@
 	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,5 +211,5 @@
     TVector<float> outCoeff(nbOfPix);
 
-    for(int pix=1;pix<nbOfPix;pix++)
+    for(int pix=0;pix<nbOfPix;pix++)
       {
 	Vector smallFDeNu(maxOfSkyComp);
@@ -216,18 +222,17 @@
 	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,35 +249,13 @@
 	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,19 +274,35 @@
   
   // 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,5 +325,5 @@
   
   //  Checking datacards
-  if (dc.NbParam("EXT_RS") < 1)   {
+  if (dc.NbParam("EXT_RS") < 2)   {
      rc = 71; 
      msg = "Invalid parameters  - Check @EXT_RS card ";
@@ -337,4 +336,6 @@
   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,5 +372,5 @@
   msg = "OK";
   nskycomp = ncomp;
-  
+  hp_nside = nside;
   // Checking for PATH definition card
   key = "RS_EXT_PATH";
@@ -409,23 +410,4 @@
   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 */
