Upgrade spectra process to select only some of the cross-spectra

src/camel/CMB/HiLLiPOP.cc

@@ -126,6 +126,8 @@ void HiLLiPOP::ProcessMultipoles(const string multipolesFile)
 {
   planck_assert(file_present(multipolesFile),string("missing file : ")+multipolesFile);
 
+  string gModeName[] = {"TT","EE","BB","TE","ET"};
+
   // Cross-spectra
   _lminXSpectra.resize(5);
   _lmaxXSpectra.resize(5);
@@ -165,7 +167,6 @@ void HiLLiPOP::ProcessMultipoles(const string multipolesFile)
   }
   input.close();
 
-
   // Cross-frequency
   _lminXFreq.resize(5);
   _lmaxXFreq.resize(5);
@@ -174,25 +175,46 @@ void HiLLiPOP::ProcessMultipoles(const string multipolesFile)
   for(unsigned int m = 0; m < 5; m++) {
     _lminXFreq[m].resize(_nXFreq,0);
     _lmaxXFreq[m].resize(_nXFreq,0);
-    _XFreqStatus[m].resize(_nXFreq);
+    _XFreqStatus[m].resize(_nXFreq,0);
     for(unsigned int f = 0; f < _nXFreq; f++) {
       unsigned int lminXFreq = 0;
       unsigned int lmaxXFreq = 0;
       for(unsigned int c = 0; c < _XFreq2XSpectra[f].size(); c++) {
-        if( _XSpectraStatus[m][_XFreq2XSpectra[f][c]] ) _XFreqStatus[m][f] = 1;
-        if( c == 0 ) {
-          lminXFreq = _lminXSpectra[m][_XFreq2XSpectra[f][c]];
-          lmaxXFreq = _lmaxXSpectra[m][_XFreq2XSpectra[f][c]];
-        }
-        else {
-          lminXFreq = min(lminXFreq,_lminXSpectra[m][_XFreq2XSpectra[f][c]]);
-          lmaxXFreq = max(lmaxXFreq,_lmaxXSpectra[m][_XFreq2XSpectra[f][c]]);
-        }
+        if( _XSpectraStatus[m][_XFreq2XSpectra[f][c]] ) {
+	  if( _XFreqStatus[m][f] == 0) {
+	    lminXFreq = _lminXSpectra[m][_XFreq2XSpectra[f][c]];
+	    lmaxXFreq = _lmaxXSpectra[m][_XFreq2XSpectra[f][c]];
+	  } else {
+	    lminXFreq = min(lminXFreq,_lminXSpectra[m][_XFreq2XSpectra[f][c]]);
+	    lmaxXFreq = max(lmaxXFreq,_lmaxXSpectra[m][_XFreq2XSpectra[f][c]]);
+	  }
+	  _XFreqStatus[m][f] = 1;
+	}
       }
       _lminXFreq[m][f] = lminXFreq;
       _lmaxXFreq[m][f] = lmaxXFreq;
     }
   }
+
+
+  for(unsigned int m = 0; m < 4; m++) {
+    if( _typeStatus[m]) {
+      cout << gModeName[m] << "XspectraStatus: ";
+      for(unsigned int c = 0; c < _nXSpectra; c++) cout << " " << _XSpectraStatus[m][c];
+      cout << endl;
+      cout << gModeName[m] << " XFreqStatus: ";
+      for(unsigned int c = 0; c < _nXFreq; c++) cout << " " << _XFreqStatus[m][c];
+      cout << endl;
+      cout << gModeName[m] << " lmin ";
+      for( unsigned int f=0; f<_nXFreq; f++)
+	cout << "\t" << _lminXFreq[m][f];
+      cout << endl;
+      cout << gModeName[m] << " lmax ";
+      for( unsigned int f=0; f<_nXFreq; f++)
+	cout << "\t" << _lmaxXFreq[m][f];
+      cout << endl;
+    }
+  }
 }
 
 
@@ -385,13 +407,20 @@ void HiLLiPOP::ProcessCovMatrix(const string pathToCovMatrix)
   }
   // End: Multipole bounds
 
+//   cout << "side " << side << endl;
+//   cout << "lrange ";
+//   for( unsigned int f=0; f<_nXFreq; f++)
+//     cout << "\t" << lrangeXFreq[0][f];
+//   cout << endl;
+
   fitshandle * input = new fitshandle();
   input->open(covMatrixFile);
   input->goto_hdu(2);
 
   unsigned int size = input->nelems(1);
-  double * cov = new double[size];
+  if( size != long(side)*long(side)) { cout << "Wrong Matrix size ! (read:" << size << " need:" << side*side << ")" << endl; exit(-1);};
 
+  double * cov = new double[size];
   input->read_column_raw(1,cov,size);
 
   for(unsigned int m1 = 0; m1 < 4; m1++) {
@@ -408,7 +437,8 @@ void HiLLiPOP::ProcessCovMatrix(const string pathToCovMatrix)
                   _invCovMat[m1][m2][f1][f2].resize(_lmaxXFreq[m1][f1]-_lminXFreq[m1][f1]+1);
                   for(unsigned int l1 = 0; l1 <= _lmaxXFreq[m1][f1]-_lminXFreq[m1][f1]; l1++) {
                     _invCovMat[m1][m2][f1][f2][l1].resize(_lmaxXFreq[m2][f2]-_lminXFreq[m2][f2]+1,0);
-                    for(unsigned int l2 = 0; l2 <= _lmaxXFreq[m2][f2]-_lminXFreq[m2][f2]; l2++) _invCovMat[m1][m2][f1][f2][l1][l2] = 1.0/4.0/M_PI/M_PI*(_lminXFreq[m1][f1]+l1)*(_lminXFreq[m1][f1]+l1+1)*(_lminXFreq[m2][f2]+l2)*(_lminXFreq[m2][f2]+l2+1)*cov[(l1+lrangeXFreq[m1][f1])*side+(l2+lrangeXFreq[m2][f2])]/1e24;
+                    for(unsigned int l2 = 0; l2 <= _lmaxXFreq[m2][f2]-_lminXFreq[m2][f2]; l2++) 
+		      _invCovMat[m1][m2][f1][f2][l1][l2] = 1.0/4.0/M_PI/M_PI*(_lminXFreq[m1][f1]+l1)*(_lminXFreq[m1][f1]+l1+1)*(_lminXFreq[m2][f2]+l2)*(_lminXFreq[m2][f2]+l2+1)*cov[(l1+lrangeXFreq[m1][f1])*side+(l2+lrangeXFreq[m2][f2])]/1e24;
                   }
                 }
               }
@@ -422,14 +452,7 @@ void HiLLiPOP::ProcessCovMatrix(const string pathToCovMatrix)
   delete [] cov;
   input->close();
   delete input;
-  /*
-   cout << "side " << side << endl;
-   cout << "lmin " << _lminXFreq[0][0] << endl;
-   cout << "lmin " << _lminXFreq[1][0] << endl;
-   cout << "lrange " << lrangeXFreq[0][_nXFreq] << endl;
-   cout << "lrange " << lrangeXFreq[1][_nXFreq] << endl;
-   cout << "matrix " << _invCovMat[0][1][0][0][200][200] << endl; // good i.e. correctly read
- */
+
 }