Changeset 899 in Sophya for trunk/SophyaProg/PMixer/skymixer.cc

Timestamp:

Apr 13, 2000, 9:44:30 AM (25 years ago)

Author:

ansari

Message:

Sophie: processing the components with a spectral index dependent
of the pixel of the map

File:

• trunk/SophyaProg/PMixer/skymixer.cc (modified) (13 diffs)

trunk/SophyaProg/PMixer/skymixer.cc

@@ -103 +103 @@
  try {
    string dcard = arg[1];
-   cout << " Decoding parameters from file " << dcard << endl;
+   if(printlev > 1) cout << " Decoding parameters from file " << dcard << endl;
    dc.ReadFile(dcard);
   
@@ -134 +134 @@
      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;
+     {
+       FITS_SphereHEALPix<float> fios(ifnm);
+       outgs = (SphereHEALPix<float>)fios;
+     }
+     if(printlev>0)
+       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;
+        cout << " MeanSig for outpout map - Mean= " << 
+          moy << " Sigma= " << sig << endl;
       }
     }
     else {
-      cout << " Output HealPix Map  created - NbPixels= " << outgs.NbPixels() << endl;
+      if(printlev>0)
+        cout << " Output HealPix Map  created - NbPixels= " << 
+          outgs.NbPixels() << endl;
       outgs.SetPixels(0.);
     }
@@ -160 +167 @@
     int sk;
     for(sk = 0; sk<nskycomp; sk++) {
+      cout << "------------------------------------" << endl;
       cout << " Processing sky component No " << sk+1 << endl;
       
@@ -168 +176 @@
       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;
-      
+      SphereHEALPix<float> ings(hp_nside);
+      {
+        FITS_SphereHEALPix<float> fiosIn(flnm);
+        ings = (SphereHEALPix<float>)fiosIn;
+      }
       if (debuglev > 4) {  // Writing the input map to the outppf
         FIO_SphereHEALPix<float> fiog(ings);
@@ -189 +191 @@
         cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
       }
+      bool mapDependentOfFreq = false;
+      key = (string)"BETAFITSFILE"+ buff;
+      if(dc.HasKey(key)) 
+        {
+          mapDependentOfFreq = true;
+        }
+
       // getting Emission spectra  
       if(!mapDependentOfFreq)
@@ -206 +215 @@
           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;
-          
+          SphereHEALPix<float> betaMap(hp_nside);
+          {
+            FITS_SphereHEALPix<float> fiosBM(flnm);
+            betaMap = (SphereHEALPix<float>)fiosBM;
+          }
+          if (debuglev > 4) {  // Writing the input map to the outppf
+            FIO_SphereHEALPix<float> fiogs(betaMap);
+            fiogs.Write(*so, key);
+          }
+
           if(nSideForInt<0) nSideForInt = sqrt((double)betaMap.NbPixels()/12);
-          
           bool bydefault = true;
           if(!bydefault)
@@ -220 +234 @@
               SphereHEALPix<float> intBetaMap(nSideForInt);
               integratedMap(*sr, betaMap, normFreq, intBetaMap);
+              if (debuglev > 4) {  // Writing the input map to the outppf
+                FIO_SphereHEALPix<float> fiogs2(intBetaMap);
+                fiogs2.Write(*so, "INTBETAMAP");
+              }
+
               betaMap.Resize(8);
               MeanSig(intBetaMap.DataBlock(), moy, sig );
@@ -226 +245 @@
               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;
-      }
-      sprintf(buff, "End of Proc. Comp. %d ", sk+1);
+      MeanSig(outgs.DataBlock(), moy, sig );
+      cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
+      cout << "-------------------------------------------------" << endl;
+      
+      sprintf(buff, "End of Processing Component %d ", sk+1);
       PrtTim(buff);
     } 
@@ -250 +267 @@
  cout << "Output Map (SphereHEALPix<float>) written to FITS file " 
       << (string)(arg[2]) << endl;
- FITS_SphereHEALPix<float> fios2(outgs);
- fios2.Write(arg[2]);
+ {
+   FITS_SphereHEALPix<float> fios2(outgs);
+   fios2.Write(arg[2]);
+ }
  PrtTim("End of WriteFITS ");
  // Saving the output map in PPF format 
@@ -417 +436 @@
     }
   if (filt == NULL) throw ParmError("datacard error ! No detection filter");
-  cout << " Filter decoded - Created " << endl;
-  cout << *filt << endl;
-
-// for debug
+  if(printlev>0)
+    {
+      cout << endl;
+      cout << " Filter decoded - Created " << endl;
+      cout << *filt << endl;
+    }
+  // for debug
   if (debuglev > 1)  SpectralResponse2Nt(*filt, *so, ppfname);
   return(filt);
@@ -542 +564 @@
                    SphereHEALPix<T>& intBetaMap)
 {
-  //  int nbpix = intBetaMap.NbPixels();
-  //  SphereHEALPix<T> newMap(sqrt((double)nbpix/12));
-  Sph2Sph(betaMap,intBetaMap);
-  for(int i=0; i<intBetaMap.NbPixels(); i++)
-    {
-      RadSpectra* rs =  new PowerLawSpectra
-        (1.,-intBetaMap(i), 0., normFreq);
-      double coeff = rs->filteredIntegratedFlux(sr);
-      intBetaMap(i) = coeff;
+  PowerLawSpectra rs(1.,-2., 0., normFreq);
+  
+  if(betaMap.NbPixels()!=intBetaMap.NbPixels()) 
+    {
+      Sph2Sph(betaMap,intBetaMap);
+      for(int i=0; i<intBetaMap.NbPixels(); i++)
+        {
+          rs.setExp(-intBetaMap(i));
+          double coeff = rs.filteredIntegratedFlux(sr);
+          intBetaMap(i) = coeff;
+        }
+    }
+  else
+    {      
+      for(int i=0; i<intBetaMap.NbPixels(); i++)
+        {
+          rs.setExp(-betaMap(i));
+          double coeff = rs.filteredIntegratedFlux(sr);
+          intBetaMap(i) = coeff;
+        }
     }
 }
@@ -575 +608 @@
   if(nbpix != intBetaMap.NbPixels())
     {
-      cout << "nbpix != intBetaMap.NbPixels()" << endl;
-      SphereHEALPix<T> bigBetaMap(sqrt((double)nbpix/12));
-      cout << "new map with size corresponding to mapToAdd" <<  sqrt((double)nbpix/12) << endl;
-      Sph2Sph(intBetaMap,bigBetaMap);
       for(int i=0; i<finalMap.NbPixels();i++)
         {
-          finalMap(i) += coeff*mapToAdd(i)*bigBetaMap(i);
+          double teta,phi,val;
+          finalMap.PixThetaPhi(i, teta, phi);
+          int pixel = intBetaMap.PixIndexSph(teta,phi);
+          val = intBetaMap.PixVal(pixel);
+          finalMap(i) += coeff*mapToAdd(i)*val;
         }
     }
@@ -638 +671 @@
   delete[] cnt;
   
-  printf("Sph2Sph_Info NbPix In= %d  Out= %d CntMoy,Sig= %g %g\n",
-         in.NbPixels(), out.NbPixels(), moy, sig);
-  PrtTim("End of  Sph2Sph() ");
+  if(printlev>2)
+    {
+      printf("Sph2Sph_Info NbPix In= %d  Out= %d CntMoy,Sig= %g %g\n",
+             in.NbPixels(), out.NbPixels(), moy, sig);
+      PrtTim("End of  Sph2Sph() ");
+    }
 }