Index: trunk/SophyaProg/PMixer/skymixer.cc
===================================================================
--- trunk/SophyaProg/PMixer/skymixer.cc	(revision 885)
+++ trunk/SophyaProg/PMixer/skymixer.cc	(revision 886)
@@ -17,5 +17,5 @@
 #include "spherehealpix.h"
 #include "fiospherehealpix.h"
-
+#include "fitsspherehealpix.h"
 #include "radspecvector.h"
 #include "blackbody.h"
@@ -101,21 +101,21 @@
 so = NULL;
 
-try {
-  string dcard = arg[1];
-  cout << " Decoding parameters from file " << dcard << endl;
-  dc.ReadFile(dcard);
-
-  rc = CheckCards(dc, msg);
-  if (rc) { 
-    cerr << " Error condition -> Rc= " << rc << endl;
-    cerr << " Msg= " << msg << endl;
-    return(rc);
-  }
-}
-catch (PException exc) {
-  msg = exc.Msg();
-  cerr << " !!!! skymixer/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
-  return(90);
-}   
+ try {
+   string dcard = arg[1];
+   cout << " Decoding parameters from file " << dcard << endl;
+   dc.ReadFile(dcard);
+  
+   rc = CheckCards(dc, msg);
+   if (rc) { 
+     cerr << " Error condition -> Rc= " << rc << endl;
+     cerr << " Msg= " << msg << endl;
+     return(rc);
+   }
+ }
+ catch (PException exc) {
+   msg = exc.Msg();
+   cerr << " !!!! skymixer/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
+   return(90);
+ }   
 
 
@@ -127,141 +127,141 @@
 if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
 
-SphereHEALPix<float> outgs(hp_nside);
-bool okout = false;
-
-try {
-  if (rdmap) {  // Reading map from FITS file
-    FitsIoServer fios;
-    char ifnm[256];
-    strncpy(ifnm, dc.SParam("READMAP", 0).c_str(), 255);
-    ifnm[255] = '\0';
-    cout << " Reading output HealPix map from FITS file " << (string)ifnm << endl;
-    fios.load(outgs, ifnm, 2);
-    cout << " Output HealPIx Map read - NbPixels= " << outgs.NbPixels() << endl;
-    if (printlev > 0) {
-      MeanSig(outgs.DataBlock(), moy, sig );
-      cout << " MeanSig for outpout map - Mean= " << moy << " Sigma= " << sig << endl;
-    }
-  }
-  else {
-    cout << " Output HealPix Map  created - NbPixels= " << outgs.NbPixels() << endl;
-    outgs.SetPixels(0.);
-  }
-
-  // Decoding detection pass-band filter 
-  sr = getSpectralResponse(dc);
-  PrtTim(" After FilterCreation ");
-
-  FitsIoServer fios; // Our FITS IO Server
-  char * flnm, buff[90];
-  string key;
-
-  SphereHEALPix<float> ings(hp_nside);  // The input map
-  double K = 1.;
-  double freqOfMap = 1.;
-  // Loop over sky component 
-  int sk;
-  for(sk = 0; sk<nskycomp; sk++) {
-    cout << " Processing sky component No " << sk+1 << endl;
+ SphereHEALPix<float> outgs(hp_nside);
+ try{
+   if (rdmap) {  // Reading map from FITS file
+     char ifnm[256];
+     strncpy(ifnm, dc.SParam("READMAP", 0).c_str(), 255);
+     ifnm[255] = '\0';
+     cout << " Reading output HealPix map from FITS file " << (string)ifnm << endl;
+     FITS_SphereHEALPix<float> fios(ifnm);
+     outgs = (SphereHEALPix<float>)fios;
+      cout << " Output HealPIx Map read - NbPixels= " << outgs.NbPixels() << endl;
+      if (printlev > 0) {
+	MeanSig(outgs.DataBlock(), moy, sig );
+	cout << " MeanSig for outpout map - Mean= " << moy << " Sigma= " << sig << endl;
+      }
+    }
+    else {
+      cout << " Output HealPix Map  created - NbPixels= " << outgs.NbPixels() << endl;
+      outgs.SetPixels(0.);
+    }
     
-    sprintf(buff, "%d", sk+1);
-    key = (string)"MAPFITSFILE" + buff;
-    flnm = BuildFITSFileName(dc.SParam(key, 0));
+    // Decoding detection pass-band filter 
+    sr = getSpectralResponse(dc);
+    PrtTim(" After FilterCreation ");
     
-    K = dc.DParam(key, 1, 1.);
+    char * flnm, buff[90];
+    string key;
     
-    bool mapDependentOfFreq = false;
-    key = (string)"BETAFITSFILE"+ buff;
-    if(dc.HasKey(key)) 
-      {
-	mapDependentOfFreq = true;
+    //    SphereHEALPix<float> ings(hp_nside);  // The input map
+    double K = 1.;
+    double freqOfMap = 1.;
+    // Loop over sky component 
+    int sk;
+    for(sk = 0; sk<nskycomp; sk++) {
+      cout << " Processing sky component No " << sk+1 << endl;
+      
+      sprintf(buff, "%d", sk+1);
+      key = (string)"MAPFITSFILE" + buff;
+      flnm = BuildFITSFileName(dc.SParam(key, 0));
+      
+      K = dc.DParam(key, 1, 1.);
+      
+      bool mapDependentOfFreq = false;
+      key = (string)"BETAFITSFILE"+ buff;
+      if(dc.HasKey(key)) 
+	{
+	  mapDependentOfFreq = true;
+	}
+      cout << "the map has a mapDependentOfFreq = " << 
+	mapDependentOfFreq << endl;
+      
+      cout << " Reading Input FITS map " << (string)flnm << endl;
+      FITS_SphereHEALPix<float> fiosIn(flnm);
+      SphereHEALPix<float> ings = (SphereHEALPix<float>)fiosIn;
+      
+      if (debuglev > 4) {  // Writing the input map to the outppf
+	FIO_SphereHEALPix<float> fiog(ings);
+	fiog.Write(*so, key);
       }
-    cout << "the map has a mapDependentOfFreq = " << 
-      mapDependentOfFreq << endl;
-
-    cout << " Reading Input FITS map " << (string)flnm << endl;
-    fios.load(ings, flnm, 2);
-    if (debuglev > 4) {  // Writing the input map to the outppf
-      FIO_SphereHEALPix<float> fiog(ings);
-      fiog.Write(*so, key);
-    }
-    if (printlev > 2) {
-      MeanSig(ings.DataBlock(), moy, sig );
-      cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
-    }
-    // getting Emission spectra  
-    if(!mapDependentOfFreq)
-      {
-	if (es) { delete es; es = NULL; }
-	es = getEmissionSpectra(dc, sk);
-	addComponent(*sr, outgs, ings, *es, K);
+      if (printlev > 2) {
+	MeanSig(ings.DataBlock(), moy, sig );
+	cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
       }
-    else
-      {
-	key = (string)"BETAFITSFILE"+ buff;
-	SphereHEALPix<float> betaMap;
-	flnm = BuildFITSFileName(dc.SParam(key, 0));
-	double normFreq = dc.DParam(key, 1, 1.);
-	if (printlev > 4) cout << "....BetaFits... normalization Freq = " << normFreq << endl;
-	int nSideForInt = dc.DParam(key, 2, 1.);
-	if (printlev > 4) cout << "....BetaFits... NSide for Integration map = " << nSideForInt << endl;
-	cout << "....BetaFits...  Reading Beta FITS map " << (string)flnm << endl;
-	fios.load(betaMap, flnm, 2);
-	if(nSideForInt<0) nSideForInt = sqrt((double)betaMap.NbPixels()/12);
-
-	bool bydefault = true;
-	if(!bydefault)
-	  addComponentBeta(outgs,ings,*sr,betaMap,normFreq, K);
-	else
-	  {
-	    // integrate the betamap over the SpectralResponse
-	    SphereHEALPix<float> intBetaMap(nSideForInt);
-	    integratedMap(*sr, betaMap, normFreq, intBetaMap);
-	    betaMap.Resize(8);
-	    MeanSig(intBetaMap.DataBlock(), moy, sig );
-	    if (printlev > 4) cout << "....BetaFits... MeanSig for intBetaMap - Mean= " << moy << " Sigma= " << sig << endl;
-	    // add the integrated beta map 
-	    addComponentBeta(outgs,ings,intBetaMap, K);
-	  }
+      // getting Emission spectra  
+      if(!mapDependentOfFreq)
+	{
+	  if (es) { delete es; es = NULL; }
+	  es = getEmissionSpectra(dc, sk);
+	  addComponent(*sr, outgs, ings, *es, K);
+	}
+      else
+	{
+	  key = (string)"BETAFITSFILE"+ buff;
+	  //SphereHEALPix<float> betaMap;
+	  flnm = BuildFITSFileName(dc.SParam(key, 0));
+	  double normFreq = dc.DParam(key, 1, 1.);
+	  if (printlev > 4) cout << "....BetaFits... normalization Freq = " << normFreq << endl;
+	  int nSideForInt = dc.DParam(key, 2, 1.);
+	  if (printlev > 4) cout << "....BetaFits... NSide for Integration map = " << nSideForInt << endl;
+	  cout << "....BetaFits...  Reading Beta FITS map " << (string)flnm << endl;
+	  
+	  FITS_SphereHEALPix<float> fiosBM(flnm);
+	  SphereHEALPix<float> betaMap = (SphereHEALPix<float>)fiosBM;
+	  
+	  if(nSideForInt<0) nSideForInt = sqrt((double)betaMap.NbPixels()/12);
+	  
+	  bool bydefault = true;
+	  if(!bydefault)
+	    addComponentBeta(outgs,ings,*sr,betaMap,normFreq, K);
+	  else
+	    {
+	      // integrate the betamap over the SpectralResponse
+	      SphereHEALPix<float> intBetaMap(nSideForInt);
+	      integratedMap(*sr, betaMap, normFreq, intBetaMap);
+	      betaMap.Resize(8);
+	      MeanSig(intBetaMap.DataBlock(), moy, sig );
+	      if (printlev > 4) cout << "....BetaFits... MeanSig for intBetaMap - Mean= " << moy << " Sigma= " << sig << endl;
+	      // add the integrated beta map 
+	      addComponentBeta(outgs,ings,intBetaMap, K);
+	    }
+	  MeanSig(outgs.DataBlock(), moy, sig );
+	  cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
+	}
+      
+      if (printlev > 1) {
 	MeanSig(outgs.DataBlock(), moy, sig );
 	cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
       }
-
-    okout = true;
-    if (printlev > 1) {
-      MeanSig(outgs.DataBlock(), moy, sig );
-      cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
-    }
-    sprintf(buff, "End of Proc. Comp. %d ", sk+1);
-    cout << " --------------------------------------- \n" << endl;    
-    PrtTim(buff);
-  }   // End of sky component loop 
-}   // End of try block
-
-catch (PException exc) {
-  msg = exc.Msg();
-  cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
-  rc = 50;
-  }   
-
-// Saving the output map in FITS format 
-if (okout) {
-  FitsIoServer fios;
-  fios.save(outgs, arg[2]);
-  cout << "Output Map (SphereHEALPix<float>) written to FITS file " 
-       << (string)(arg[2]) << endl;
-  PrtTim("End of WriteFITS ");
-  // Saving the output map in PPF format 
-  if (narg > 3) {
+      sprintf(buff, "End of Proc. Comp. %d ", sk+1);
+      PrtTim(buff);
+    } 
+ }
+ catch(PException exc) 
+   {
+     cout << "catched PException" << endl;
+     msg = exc.Msg();
+     cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
+     rc = 50;
+     return(50);
+   } 
+ 
+ // Saving the output map in FITS format 
+ cout << "Output Map (SphereHEALPix<float>) written to FITS file " 
+      << (string)(arg[2]) << endl;
+ FITS_SphereHEALPix<float> fios2(outgs);
+ fios2.Write(arg[2]);
+ PrtTim("End of WriteFITS ");
+ // Saving the output map in PPF format 
+ if (narg > 3) {
    POutPersist s(arg[3]); 
    FIO_SphereHEALPix<float> fiog(&outgs) ;
    fiog.Write(s);
    cout << "Output Map (SphereHEALPix<float>) written to POutPersist file "  
-        << (string)(arg[3]) << endl;
+  	<< (string)(arg[3]) << endl;
    PrtTim("End of WritePPF ");
-   }
-  }
-  if (so) delete so;  //  Closing the debug ppf file 
-  return(rc);
+ }
+ if (so) delete so;  //  Closing the debug ppf file 
+ return(rc);
 }
 
