Commit 6624368a authored by Matthieu Tristram's avatar Matthieu Tristram
Browse files

WARNING: change of parameters !!!

No need for all foreground parameters:
- when considering TE: only AdustTP
- when considering EE: only AdustPP
- when considering TT: AdustTT, Asz, Acib, Aksz, Aszxcib, and all Aps
parent 63469ed9
...@@ -472,17 +472,18 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC ...@@ -472,17 +472,18 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC
for(unsigned int i = 0; i < _nMap; i++) {sprintf(tmpPar,"c%d",i); _n.push_back(tmpPar);} for(unsigned int i = 0; i < _nMap; i++) {sprintf(tmpPar,"c%d",i); _n.push_back(tmpPar);}
// End: Calibration // End: Calibration
//SP 9/3/15: global calibration factor //global calibration factor
_n.push_back("A_planck"); _n.push_back("A_planck");
//SP comment
// Begin: First beam eigenmode amplitudes // Begin: First beam eigenmode amplitudes
//for(unsigned int c = 0; c < _nXSpectra; c++) {sprintf(tmpPar,"beta%d",c); _n.push_back(tmpPar);} //for(unsigned int c = 0; c < _nXSpectra; c++) {sprintf(tmpPar,"beta%d",c); _n.push_back(tmpPar);}
// End: First beam eigenmode amplitudes // End: First beam eigenmode amplitudes
// Begin: Amplitudes of the point sources power spectrum // Begin: Amplitudes of the point sources power spectrum
for(unsigned int f = 0; f < _nXFreq; f++) {sprintf(tmpPar,"Aps%s",_XFreq[f].c_str()); _n.push_back(tmpPar);} if( _typeStatus[0]) {
for(unsigned int f = 0; f < _nXFreq; f++) {sprintf(tmpPar,"Aps%s",_XFreq[f].c_str()); _n.push_back(tmpPar);}
}
// End: Amplitudes of the point sources power spectrum // End: Amplitudes of the point sources power spectrum
...@@ -497,36 +498,38 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC ...@@ -497,36 +498,38 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC
if( it != freqPlanckHFI.end() ) _fnu[f][i] = fnu[it-freqPlanckHFI.begin()]; if( it != freqPlanckHFI.end() ) _fnu[f][i] = fnu[it-freqPlanckHFI.begin()];
} }
} }
// tSZ spectrum template, in Dl=l(l+1)/2pi Cl, units uK at 143GHz // tSZ spectrum template, in Dl=l(l+1)/2pi Cl, units uK at 143GHz
planck_assert(file_present(SZFile),string("missing file : ")+SZFile); if( _typeStatus[0]) {
planck_assert(file_present(SZFile),string("missing file : ")+SZFile);
fitshandle * SZInput = new fitshandle();
SZInput->open(SZFile); fitshandle * SZInput = new fitshandle();
SZInput->goto_hdu(2); SZInput->open(SZFile);
SZInput->goto_hdu(2);
unsigned int sizeSZ = SZInput->nelems(2);
double * ellSZ = new double[sizeSZ]; unsigned int sizeSZ = SZInput->nelems(2);
double * Cl1haloSZ = new double[sizeSZ]; double * ellSZ = new double[sizeSZ];
double * Cl2haloSZ = new double[sizeSZ]; double * Cl1haloSZ = new double[sizeSZ];
SZInput->read_column_raw(1,ellSZ,sizeSZ); double * Cl2haloSZ = new double[sizeSZ];
SZInput->read_column_raw(2,Cl1haloSZ,sizeSZ); SZInput->read_column_raw(1,ellSZ,sizeSZ);
SZInput->read_column_raw(3,Cl2haloSZ,sizeSZ); SZInput->read_column_raw(2,Cl1haloSZ,sizeSZ);
SZInput->read_column_raw(3,Cl2haloSZ,sizeSZ);
_ClSZ.resize(_nXFreq);
for(unsigned int f = 0; f < _nXFreq; f++) { _ClSZ.resize(_nXFreq);
_ClSZ[f].resize(_maxOflmax[0]+1,0); for(unsigned int f = 0; f < _nXFreq; f++) {
for(unsigned int l = 0; l <= _maxOflmax[0]-ellSZ[0]; l++) _ClSZ[f][ellSZ[l]] = 2*M_PI/ellSZ[l]/(ellSZ[l]+1)*(Cl1haloSZ[l]+Cl2haloSZ[l])*_fnu[f][0]*_fnu[f][1]/fnu[1]/fnu[1]; _ClSZ[f].resize(_maxOflmax[0]+1,0);
for(unsigned int l = 0; l <= _maxOflmax[0]-ellSZ[0]; l++) _ClSZ[f][ellSZ[l]] = 2*M_PI/ellSZ[l]/(ellSZ[l]+1)*(Cl1haloSZ[l]+Cl2haloSZ[l])*_fnu[f][0]*_fnu[f][1]/fnu[1]/fnu[1];
}
SZInput->close();
delete [] ellSZ;
delete [] Cl1haloSZ;
delete [] Cl2haloSZ;
delete SZInput;
_n.push_back("Asz");
} }
SZInput->close();
delete [] ellSZ;
delete [] Cl1haloSZ;
delete [] Cl2haloSZ;
delete SZInput;
_n.push_back("Asz");
// End : tSZ // End : tSZ
...@@ -543,39 +546,41 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC ...@@ -543,39 +546,41 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC
} }
// CIB spectra // CIB spectra
_ClCIB.resize(_nXFreq); if( _typeStatus[0]) {
for(unsigned int f = 0; f < _nXFreq; f++) { _ClCIB.resize(_nXFreq);
//sprintf(tmpFile,"%s_%s.fits",pathToCIBSpectra.c_str(),_XFreq[f].c_str()); for(unsigned int f = 0; f < _nXFreq; f++) {
//planck_assert(file_present(tmpFile),string("missing file : ")+tmpFile); //sprintf(tmpFile,"%s_%s.fits",pathToCIBSpectra.c_str(),_XFreq[f].c_str());
//planck_assert(file_present(tmpFile),string("missing file : ")+tmpFile);
ostringstream os ;
os << pathToCIBSpectra << "_" << _XFreq[f] << ".fits"; ostringstream os ;
tmpFile=Parser::CheckPath(os.str()); os << pathToCIBSpectra << "_" << _XFreq[f] << ".fits";
tmpFile=Parser::CheckPath(os.str());
fitshandle * CIBInput = new fitshandle();
CIBInput->open(tmpFile); fitshandle * CIBInput = new fitshandle();
CIBInput->goto_hdu(2); CIBInput->open(tmpFile);
CIBInput->goto_hdu(2);
unsigned int sizeCIB = CIBInput->nelems(2);
double * ellCIB = new double[sizeCIB]; unsigned int sizeCIB = CIBInput->nelems(2);
double * Cl1haloCIB = new double[sizeCIB]; double * ellCIB = new double[sizeCIB];
double * Cl2haloCIB = new double[sizeCIB]; double * Cl1haloCIB = new double[sizeCIB];
CIBInput->read_column_raw(1,ellCIB,sizeCIB); double * Cl2haloCIB = new double[sizeCIB];
CIBInput->read_column_raw(2,Cl1haloCIB,sizeCIB); CIBInput->read_column_raw(1,ellCIB,sizeCIB);
CIBInput->read_column_raw(3,Cl2haloCIB,sizeCIB); CIBInput->read_column_raw(2,Cl1haloCIB,sizeCIB);
CIBInput->read_column_raw(3,Cl2haloCIB,sizeCIB);
_ClCIB[f].resize(_maxOflmax[0]+1);
for(unsigned int l = 0; l <= _maxOflmax[0]-ellCIB[0]; l++) _ClCIB[f][ellCIB[l]] = (Cl1haloCIB[l]+Cl2haloCIB[l])/_gnu[f][0]/_gnu[f][1]; _ClCIB[f].resize(_maxOflmax[0]+1);
for(unsigned int l = 0; l <= _maxOflmax[0]-ellCIB[0]; l++) _ClCIB[f][ellCIB[l]] = (Cl1haloCIB[l]+Cl2haloCIB[l])/_gnu[f][0]/_gnu[f][1];
CIBInput->close();
CIBInput->close();
delete [] ellCIB;
delete [] Cl1haloCIB;
delete [] Cl2haloCIB;
delete CIBInput;
}
delete [] ellCIB; _n.push_back("Acib");
delete [] Cl1haloCIB;
delete [] Cl2haloCIB;
delete CIBInput;
} }
_n.push_back("Acib");
// End : CIB // End : CIB
...@@ -630,71 +635,75 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC ...@@ -630,71 +635,75 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC
delete DustInput; delete DustInput;
} }
_n.push_back("AdustTT"); if( _typeStatus[0]) _n.push_back("AdustTT");
_n.push_back("AdustPP"); if( _typeStatus[1] || _typeStatus[2]) _n.push_back("AdustPP");
_n.push_back("AdustTP"); if( _typeStatus[3]) _n.push_back("AdustTP");
// End: Dust // End: Dust
// Begin : kSZ // Begin : kSZ
planck_assert(file_present(kSZFile),string("missing file : ")+kSZFile); if( _typeStatus[0]) {
planck_assert(file_present(kSZFile),string("missing file : ")+kSZFile);
fitshandle * kSZInput = new fitshandle();
kSZInput->open(kSZFile); fitshandle * kSZInput = new fitshandle();
kSZInput->goto_hdu(2); kSZInput->open(kSZFile);
kSZInput->goto_hdu(2);
//kSZ template, Dl=l(l+1)/2pi Cl, units uK
unsigned int sizekSZ = kSZInput->nelems(2); //kSZ template, Dl=l(l+1)/2pi Cl, units uK
double * ellkSZ = new double[sizekSZ]; unsigned int sizekSZ = kSZInput->nelems(2);
double * ClkSZ = new double[sizekSZ]; double * ellkSZ = new double[sizekSZ];
kSZInput->read_column_raw(1,ellkSZ,sizekSZ); double * ClkSZ = new double[sizekSZ];
kSZInput->read_column_raw(2,ClkSZ,sizekSZ); kSZInput->read_column_raw(1,ellkSZ,sizekSZ);
kSZInput->read_column_raw(2,ClkSZ,sizekSZ);
_ClkSZ.resize(_maxOflmax[0]+1,0);
for(unsigned int l = 0; l <= _maxOflmax[0]-ellkSZ[0]; l++) _ClkSZ[ellkSZ[l]] = 2*M_PI/ellkSZ[l]/(ellkSZ[l]+1)*ClkSZ[l]; _ClkSZ.resize(_maxOflmax[0]+1,0);
for(unsigned int l = 0; l <= _maxOflmax[0]-ellkSZ[0]; l++) _ClkSZ[ellkSZ[l]] = 2*M_PI/ellkSZ[l]/(ellkSZ[l]+1)*ClkSZ[l];
kSZInput->close();
kSZInput->close();
delete [] ellkSZ;
delete [] ClkSZ; delete [] ellkSZ;
delete kSZInput; delete [] ClkSZ;
delete kSZInput;
_n.push_back("Aksz");
_n.push_back("Aksz");
}
// End : kSZ // End : kSZ
// Begin: tSZxCIB spectra // Begin: tSZxCIB spectra
_ClSZxCIB.resize(_nXFreq); if( _typeStatus[0]) {
for(unsigned int f = 0; f < _nXFreq; f++) { _ClSZxCIB.resize(_nXFreq);
//sprintf(tmpFile,"%s_%s.fits",pathToSZxCIBSpectra.c_str(),_XFreq[f].c_str()); for(unsigned int f = 0; f < _nXFreq; f++) {
//planck_assert(file_present(tmpFile),string("missing file : ")+tmpFile); //sprintf(tmpFile,"%s_%s.fits",pathToSZxCIBSpectra.c_str(),_XFreq[f].c_str());
//planck_assert(file_present(tmpFile),string("missing file : ")+tmpFile);
ostringstream os ;
os << pathToSZxCIBSpectra << "_" << _XFreq[f] << ".fits"; ostringstream os ;
tmpFile=Parser::CheckPath(os.str()); os << pathToSZxCIBSpectra << "_" << _XFreq[f] << ".fits";
tmpFile=Parser::CheckPath(os.str());
fitshandle * SZxCIBInput = new fitshandle();
SZxCIBInput->open(tmpFile); fitshandle * SZxCIBInput = new fitshandle();
SZxCIBInput->goto_hdu(2); SZxCIBInput->open(tmpFile);
SZxCIBInput->goto_hdu(2);
unsigned int sizeSZxCIB = SZxCIBInput->nelems(2);
double * ellSZxCIB = new double[sizeSZxCIB]; unsigned int sizeSZxCIB = SZxCIBInput->nelems(2);
double * ClSZxCIB = new double[sizeSZxCIB]; double * ellSZxCIB = new double[sizeSZxCIB];
SZxCIBInput->read_column_raw(1,ellSZxCIB,sizeSZxCIB); double * ClSZxCIB = new double[sizeSZxCIB];
SZxCIBInput->read_column_raw(2,ClSZxCIB,sizeSZxCIB); SZxCIBInput->read_column_raw(1,ellSZxCIB,sizeSZxCIB);
SZxCIBInput->read_column_raw(2,ClSZxCIB,sizeSZxCIB);
_ClSZxCIB[f].resize(_maxOflmax[0]+1);
for(unsigned int l = 0; l <= _maxOflmax[0]-ellSZxCIB[0]; l++) _ClSZxCIB[f][ellSZxCIB[l]] = ClSZxCIB[l]/_gnu[f][0]/_gnu[f][1]; _ClSZxCIB[f].resize(_maxOflmax[0]+1);
for(unsigned int l = 0; l <= _maxOflmax[0]-ellSZxCIB[0]; l++) _ClSZxCIB[f][ellSZxCIB[l]] = ClSZxCIB[l]/_gnu[f][0]/_gnu[f][1];
SZxCIBInput->close();
SZxCIBInput->close();
delete [] ellSZxCIB;
delete [] ClSZxCIB; delete [] ellSZxCIB;
delete SZxCIBInput; delete [] ClSZxCIB;
delete SZxCIBInput;
}
_n.push_back("Aszxcib");
} }
_n.push_back("Aszxcib");
// End : SZxCIB // End : SZxCIB
#ifdef PolarizationEfficiency #ifdef PolarizationEfficiency
...@@ -944,6 +953,9 @@ void HiLLiPOP::computeResiduals( vector<double> ClCMB, vector<double>& nuisance, ...@@ -944,6 +953,9 @@ void HiLLiPOP::computeResiduals( vector<double> ClCMB, vector<double>& nuisance,
double foregrounds = 0; double foregrounds = 0;
unsigned int modeOffset = 0; unsigned int modeOffset = 0;
char tmpPar[32]; char tmpPar[32];
vector<double> Aps(_nXFreq);
double Asz, Acib, Aksz, Aszxcib;
double AdustTT, AdustTP, AdustPP;
_residual.resize(5); _residual.resize(5);
_instrumental.resize(5); _instrumental.resize(5);
...@@ -965,28 +977,29 @@ void HiLLiPOP::computeResiduals( vector<double> ClCMB, vector<double>& nuisance, ...@@ -965,28 +977,29 @@ void HiLLiPOP::computeResiduals( vector<double> ClCMB, vector<double>& nuisance,
#endif #endif
// Amplitudes of the point sources // Amplitudes of the point sources
vector<double> Aps(_nXFreq); if( _typeStatus[0]) {
for(unsigned int f = 0; f < _nXFreq; f++) { for(unsigned int f = 0; f < _nXFreq; f++) {
sprintf(tmpPar,"Aps%s",_XFreq[f].c_str()); sprintf(tmpPar,"Aps%s",_XFreq[f].c_str());
Aps[f] = nuisance[getIndex(tmpPar)]; Aps[f] = nuisance[getIndex(tmpPar)];
}
} }
// tSZ amplitude // tSZ amplitude
double Asz = nuisance[getIndex("Asz")]; if( _typeStatus[0]) Asz = nuisance[getIndex("Asz")];
// CIB amplitude // CIB amplitude
double Acib = nuisance[getIndex("Acib")]; if( _typeStatus[0]) Acib = nuisance[getIndex("Acib")];
// Dust amplitude // Dust amplitude
double AdustTT = nuisance[getIndex("AdustTT")]; if( _typeStatus[0]) AdustTT = nuisance[getIndex("AdustTT")];
double AdustPP = nuisance[getIndex("AdustPP")]; if( _typeStatus[1] || _typeStatus[2]) AdustPP = nuisance[getIndex("AdustPP")];
double AdustTP = nuisance[getIndex("AdustTP")]; if( _typeStatus[3]) AdustTP = nuisance[getIndex("AdustTP")];
// kSZ amplitude // kSZ amplitude
double Aksz = nuisance[getIndex("Aksz")]; if( _typeStatus[0]) Aksz = nuisance[getIndex("Aksz")];
// tSZxCIB amplitude // tSZxCIB amplitude
double Aszxcib = nuisance[getIndex("Aszxcib")]; if( _typeStatus[0]) Aszxcib = nuisance[getIndex("Aszxcib")];
#ifdef PolarizationEfficiency #ifdef PolarizationEfficiency
// Polarization efficiency // Polarization efficiency
...@@ -1298,14 +1311,16 @@ void HiLLiPOP::WriteOutput(const vector<double>& ClCMB,const vector<double>& nui ...@@ -1298,14 +1311,16 @@ void HiLLiPOP::WriteOutput(const vector<double>& ClCMB,const vector<double>& nui
of << "\t Residual" << CMBmode[m] << _XFreq[f]; of << "\t Residual" << CMBmode[m] << _XFreq[f];
//of << "\t Instrumental" << CMBmode[m] << _XFreq[f]; //of << "\t Instrumental" << CMBmode[m] << _XFreq[f];
} }
for(unsigned int f = 0; f < _nXFreq; f++) of << "\t PS" << _XFreq[f]; if( _typeStatus[0]) {
for(unsigned int f = 0; f < _nXFreq; f++) of << "\t SZ" << _XFreq[f]; for(unsigned int f = 0; f < _nXFreq; f++) of << "\t PS" << _XFreq[f];
for(unsigned int f = 0; f < _nXFreq; f++) of << "\t CIB" << _XFreq[f]; for(unsigned int f = 0; f < _nXFreq; f++) of << "\t SZ" << _XFreq[f];
for(unsigned int f = 0; f < _nXFreq; f++) of << "\t CIB" << _XFreq[f];
of << "\t kSZ";
for(unsigned int f = 0; f < _nXFreq; f++) of << "\t SZxCIB" << _XFreq[f];
}
for(unsigned int m = 0; m < 4; m++) for(unsigned int m = 0; m < 4; m++)
if( _modeStatus[m] || (m==3 && (_modeStatus[3] || _modeStatus[4]))) if( _modeStatus[m] || (m==3 && (_modeStatus[3] || _modeStatus[4])))
for(unsigned int f = 0; f < _nXFreq; f++) of << "\t Dust" << CMBmode[m] << _XFreq[f]; for(unsigned int f = 0; f < _nXFreq; f++) of << "\t Dust" << CMBmode[m] << _XFreq[f];
of << "\t kSZ";
for(unsigned int f = 0; f < _nXFreq; f++) of << "\t SZxCIB" << _XFreq[f];
of << endl; of << endl;
// End: header // End: header
...@@ -1319,12 +1334,16 @@ void HiLLiPOP::WriteOutput(const vector<double>& ClCMB,const vector<double>& nui ...@@ -1319,12 +1334,16 @@ void HiLLiPOP::WriteOutput(const vector<double>& ClCMB,const vector<double>& nui
else of << " " << 0;// << " " << 0; else of << " " << 0;// << " " << 0;
} }
} }
for(unsigned int f = 0; f < _nXFreq; f++) { if( _typeStatus[0]) {
sprintf(tmpPar,"Aps%s",_XFreq[f].c_str()); for(unsigned int f = 0; f < _nXFreq; f++) {
of << " " << nuisance[getIndex(tmpPar)]; sprintf(tmpPar,"Aps%s",_XFreq[f].c_str());
of << " " << nuisance[getIndex(tmpPar)];
}
for(unsigned int f = 0; f < _nXFreq; f++) of << " " << nuisance[getIndex("Asz")]*_ClSZ[f][l];
for(unsigned int f = 0; f < _nXFreq; f++) of << " " << nuisance[getIndex("Acib")]*_ClCIB[f][l];
of << " " << nuisance[getIndex("Aksz")]*_ClkSZ[l];
for(unsigned int f = 0; f < _nXFreq; f++) of << " " << nuisance[getIndex("Aszxcib")]*_ClSZxCIB[f][l];
} }
for(unsigned int f = 0; f < _nXFreq; f++) of << " " << nuisance[getIndex("Asz")]*_ClSZ[f][l];
for(unsigned int f = 0; f < _nXFreq; f++) of << " " << nuisance[getIndex("Acib")]*_ClCIB[f][l];
for(unsigned int m = 0; m < 4; m++) { for(unsigned int m = 0; m < 4; m++) {
if( _modeStatus[m] || (m==3 && (_modeStatus[3] || _modeStatus[4]))) if( _modeStatus[m] || (m==3 && (_modeStatus[3] || _modeStatus[4])))
for(unsigned int f = 0; f < _nXFreq; f++) { for(unsigned int f = 0; f < _nXFreq; f++) {
...@@ -1335,8 +1354,6 @@ void HiLLiPOP::WriteOutput(const vector<double>& ClCMB,const vector<double>& nui ...@@ -1335,8 +1354,6 @@ void HiLLiPOP::WriteOutput(const vector<double>& ClCMB,const vector<double>& nui
if( m == 3 && _modeStatus[3] && _modeStatus[4] ) of << " " << nuisance[getIndex("AdustTP")]*(_ClDust[f][3*(_maxOflmax[0]+1)+l]+_ClDust[f][4*(_maxOflmax[0]+1)+l])/2.; if( m == 3 && _modeStatus[3] && _modeStatus[4] ) of << " " << nuisance[getIndex("AdustTP")]*(_ClDust[f][3*(_maxOflmax[0]+1)+l]+_ClDust[f][4*(_maxOflmax[0]+1)+l])/2.;
} }
} }
of << " " << nuisance[getIndex("Aksz")]*_ClkSZ[l];
for(unsigned int f = 0; f < _nXFreq; f++) of << " " << nuisance[getIndex("Aszxcib")]*_ClSZxCIB[f][l];
of << endl; of << endl;
} }
of.close(); of.close();
......
...@@ -498,97 +498,103 @@ void Hillipop::ProcessNuisanceVariables(const string SZFile,const string pathToC ...@@ -498,97 +498,103 @@ void Hillipop::ProcessNuisanceVariables(const string SZFile,const string pathToC
#endif #endif
// Begin: Amplitudes of the point sources power spectrum (MJy^2.sr-1) // Begin: Amplitudes of the point sources power spectrum (MJy^2.sr-1)
double radio[6] = { 7.76, 5.36, 4.26, 4.83, 3.6, 3.22}; if( _typeStatus[0]) {
double dusty[6] = { 0.18179145, 0.47110054, 1.90295795, 1.24375538, 5.06409993, 20.99528025}; double radio[6] = { 7.76, 5.36, 4.26, 4.83, 3.6, 3.22};
_ClpsRadio.resize( _nXFreq); double dusty[6] = { 0.18179145, 0.47110054, 1.90295795, 1.24375538, 5.06409993, 20.99528025};
_ClpsDusty.resize( _nXFreq); _ClpsRadio.resize( _nXFreq);
for(unsigned int f = 0; f < _nXFreq; f++) { _ClpsDusty.resize( _nXFreq);
_ClpsRadio[f] = radio[f] / _gnu[f][0] / _gnu[f][1]; for(unsigned int f = 0; f < _nXFreq; f++) {
_ClpsDusty[f] = dusty[f] / _gnu[f][0] / _gnu[f][1]; _ClpsRadio[f] = radio[f] / _gnu[f][0] / _gnu[f][1];
_ClpsDusty[f] = dusty[f] / _gnu[f][0] / _gnu[f][1];
}
_n.push_back( "Aradio");
_n.push_back( "Adusty");
} }
_n.push_back( "Aradio");
_n.push_back( "Adusty");
// End: Amplitudes of the point sources power spectrum // End: Amplitudes of the point sources power spectrum
// Begin : tSZ // Begin : tSZ
// frequency dependence (100,143,217,353,545,857) GHz of the SZ effect // frequency dependence (100,143,217,353,545,857) GHz of the SZ effect
double fnu[6] = {-4.031, -2.785, 0.187, 6.205, 14.455, 26.335}; if( _typeStatus[0]) {
_fnu.resize(_nXFreq); double fnu[6] = {-4.031, -2.785, 0.187, 6.205, 14.455, 26.335};
for(unsigned int f = 0; f < _nXFreq; f++) { _fnu.resize(_nXFreq);
_fnu[f].resize(2); for(unsigned int f = 0; f < _nXFreq; f++) {
for(unsigned int i = 0; i < 2; i++) { _fnu[f].resize(2);
it = find(freqPlanckHFI.begin(),freqPlanckHFI.end(),_XFreq2Freq[f][i].c_str()); for(unsigned int i = 0; i < 2; i++) {
if( it != freqPlanckHFI.end() ) _fnu[f][i] = fnu[it-freqPlanckHFI.begin()]; it = find(freqPlanckHFI.begin(),freqPlanckHFI.end(),_XFreq2Freq[f][i].c_str());
if( it != freqPlanckHFI.end() ) _fnu[f][i] = fnu[it-freqPlanckHFI.begin()];
}
} }
}
// tSZ spectrum template, in Dl=l(l+1)/2pi Cl, units uK at 143GHz
// tSZ spectrum template, in Dl=l(l+1)/2pi Cl, units uK at 143GHz planck_assert(file_present(SZFile),string("missing file : ")+SZFile);
planck_assert(file_present(SZFile),string("missing file : ")+SZFile);
fitshandle * SZInput = new fitshandle(); fitshandle * SZInput = new fitshandle();
SZInput->open(SZFile); SZInput->open(SZFile);
SZInput->goto_hdu(2); SZInput->goto_hdu(2);
unsigned int sizeSZ = SZInput->nelems(2);
double * ellSZ = new double[sizeSZ];
double * Cl1haloSZ = new double[sizeSZ];
double * Cl2haloSZ = new double[sizeSZ];
SZInput->read_column_raw(1,ellSZ,sizeSZ);
SZInput->read_column_raw(2,Cl1haloSZ,sizeSZ);
SZInput->read_column_raw(3,Cl2haloSZ,sizeSZ);
_ClSZ.resize(_nXFreq);
for(unsigned int f = 0; f < _nXFreq; f++) {
_ClSZ[f].resize(_maxOflmax[0]+1,0);
for(unsigned int l = 0; l <= _maxOflmax[0]-ellSZ[0]; l++) _ClSZ[f][ellSZ[l]] = 2*M_PI/ellSZ[l]/(ellSZ[l]+1)*(Cl1haloSZ[l]+Cl2haloSZ[l])*_fnu[f][0]*_fnu[f][1]/fnu[1]/fnu[1];
}
unsigned int sizeSZ = SZInput->nelems(2); SZInput->close();
double * ellSZ = new double[sizeSZ];
double * Cl1haloSZ = new double[sizeSZ]; delete [] ellSZ;
double * Cl2haloSZ = new double[sizeSZ]; delete [] Cl1haloSZ;
SZInput->read_column_raw(1,ellSZ,sizeSZ); delete [] Cl2haloSZ;
SZInput->read_column_raw(2,Cl1haloSZ,sizeSZ); delete SZInput;
SZInput->read_column_raw(3,Cl2haloSZ,sizeSZ);
_ClSZ.resize(_nXFreq); _n.push_back("Asz");
for(unsigned int f = 0; f < _nXFreq; f++) {
_ClSZ[f].resize(_maxOflmax[0]+1,0);
for(unsigned int l = 0; l <= _maxOflmax[0]-ellSZ[0]; l++) _ClSZ[f][ellSZ[l]] = 2*M_PI/ellSZ[l]/(ellSZ[l]+1)*(Cl1haloSZ[l]+Cl2haloSZ[l])*_fnu[f][0]*_fnu[f][1]/fnu[1]/fnu[1];
} }
SZInput->close();
delete [] ellSZ;
delete [] Cl1haloSZ;
delete [] Cl2haloSZ;
delete SZInput;
_n.push_back("Asz");
// End : tSZ // End : tSZ
// Begin : CIB // Begin : CIB
_ClCIB.resize(_nXFreq); if( _typeStatus[0]) {
for(unsigned int f = 0; f < _nXFreq; f++) { _ClCIB.resize(_nXFreq);
//sprintf(tmpFile,"%s_%s.fits",pathToCIBSpectra.c_str(),_XFreq[f].c_str()); for(unsigned int f = 0; f < _nXFreq; f++) {
//planck_assert(file_present(tmpFile),string("missing file : ")+tmpFile); //sprintf(tmpFile,"%s_%s.fits",pathToCIBSpectra.c_str(),_XFreq[f].c_str());
//planck_assert(file_present(tmpFile),string("missing file : ")+tmpFile);
ostringstream os ;
os << pathToCIBSpectra << "_" << _XFreq[f] << ".fits"; ostringstream os ;
tmpFile=Parser::CheckPath(os.str()); os << pathToCIBSpectra << "_" << _XFreq[f] << ".fits";
tmpFile=Parser::CheckPath(os.str());
fitshandle * CIBInput = new fitshandle();
CIBInput->open(tmpFile); fitshandle * CIBInput = new fitshandle();
CIBInput->goto_hdu(2); CIBInput->open(tmpFile);
CIBInput