Skip to content

GitLab

CAMEL
CAMEL

📝 User registration policy changed. Please read the IN2P3 Gitlab documentation and more specifically the User Sign-up section. ✅

Commit 57b8f7d1 authored by Matthieu Tristram's avatar Matthieu Tristram
Browse files
Options

Merge branch 'fix-warnings' into 'master' 

Fix warnings

See merge request !5
parents 2cef189c 28ed2caa
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 547 additions and 538 deletions
src/AbsRand/OofNoise.cc
View file @ 57b8f7d1
......@@ -21,68 +21,67 @@
using namespace std;
//--------------------
// C
// C
//----------------
//---------------
// Constructors --
//----------------
OofNoise::OofNoise(AbsGauss* r,double slope,double fmin,double fknee,double fsample):AbsGauss(0.,1.),_r(r),_sigma(r->sigma())
OofNoise::OofNoise(AbsGauss* r,double slope,double fmin,double fknee,double fsample):AbsGauss(0.,1.),_r(r),_sigma(r->sigma())
{
//TW:
double w0=2*M_PI*fmin;
double w1=2*M_PI*fknee;
double Wmax=log10(w1);
double Wmin=log10(w0);
//TW:
double w0=2*M_PI*fmin;
double w1=2*M_PI*fknee;
double Wmax=log10(w1);
double Wmin=log10(w0);
//!atention! j'utilise l'inverse
double a=-slope;
//!atention! j'utilise l'inverse
double a=-slope;
//roughly 1.5pole /per decade
int Nproc=(int)((Wmax-Wmin)*2.+log10(fsample));
//roughly 1.5pole /per decade
int Nproc=(int)((Wmax-Wmin)*2.+log10(fsample));
vector<double> p(Nproc);
vector<double> z(Nproc);
double dp=(Wmax-Wmin)/Nproc;
p[0]=Wmin+(1-a/2)/2*dp;
z[0]=p[0]+a/2*dp;
vector<double> p(Nproc);
vector<double> z(Nproc);
for(int i=1;i<Nproc;i++){
p[i]=p[i-1]+dp;
z[i]=p[i]+a/2*dp;
}
double dp=(Wmax-Wmin)/Nproc;
p[0]=Wmin+(1-a/2)/2*dp;
z[0]=p[0]+a/2*dp;
vector<double> pole(Nproc);
vector<double> zero(Nproc);
for (int i=0;i<Nproc;i++){
pole[i]=POW(10.,p[i]);
zero[i]=POW(10.,z[i]);
}
//filters
for (int i=0;i<Nproc;i++) {
for(int i=1;i<Nproc;i++){
p[i]=p[i-1]+dp;
z[i]=p[i]+a/2*dp;
}
Oof2Noise* f2=new Oof2Noise(0,pole[i]/(2*M_PI),zero[i]/(2*M_PI),fsample);
ff.push_back(f2);
printf("OOFNOISE: filter %d : pole=%f zero=%f\n",i,pole[i],zero[i]);
vector<double> pole(Nproc);
vector<double> zero(Nproc);
for (int i=0;i<Nproc;i++){
pole[i]=POW(10.,p[i]);
zero[i]=POW(10.,z[i]);
}
}
//filters
for (int i=0;i<Nproc;i++) {
//_sigma*=sqrt(fsample);
Oof2Noise* f2=new Oof2Noise(0,pole[i]/(2*M_PI),zero[i]/(2*M_PI),fsample);
ff.push_back(f2);
printf("OOFNOISE: filter %d : pole=%f zero=%f\n",i,pole[i],zero[i]);
}
//_sigma*=sqrt(fsample);
}
OofNoise::~OofNoise(){
for (unsigned i=0;i<ff.size();i++) delete ff[i];
if (_r) delete _r; _r=NULL;
};
OofNoise::~OofNoise(){
for (unsigned i=0;i<ff.size();i++) delete ff[i];
if (_r) delete _r; _r=NULL;
}
double
double
OofNoise::normal(){
double x2=_r->normal();
for (unsigned i=0;i<ff.size();i++) x2=ff[i]->filter(x2);
return x2*_sigma;
};
}
src/camel/CMB/HiLLiPOP.cc
View file @ 57b8f7d1
......@@ -38,7 +38,7 @@ void HiLLiPOP::Init(const string fileName)
cout << "Using: " << fileName << endl;
paramfile parameters(fileName);
string multipolesFile = Parser::CheckPath(parameters.find<string>("MultipolesRange"));
string SZFile = Parser::CheckPath(parameters.find<string>("SZ"));
string kSZFile = Parser::CheckPath(parameters.find<string>("kSZ"));
......@@ -67,7 +67,7 @@ void HiLLiPOP::Init(const string fileName)
_typeStatus[0] = _modeStatus[0];
_typeStatus[1] = _modeStatus[1];
_typeStatus[2] = _modeStatus[2];
if( _modeStatus[3] == 1 || _modeStatus[4] == 1) _typeStatus[3] = 1; else _typeStatus[3] = 0;
if( _modeStatus[3] == 1 || _modeStatus[4] == 1) _typeStatus[3] = 1; else _typeStatus[3] = 0;
ProduceList();
......@@ -141,7 +141,7 @@ void HiLLiPOP::ProcessMultipoles(const string multipolesFile)
for(unsigned int m = 0; m < 5; m++) {
if( m < 4) type = m; else type = 3;
input.goto_hdu(type+2);
input.goto_hdu(type+2);
_lminXSpectra[m].resize(_nXSpectra);
_lmaxXSpectra[m].resize(_nXSpectra);
......@@ -256,7 +256,7 @@ void HiLLiPOP::ProcessXSpectra(const string pathToXSpectra)
// TT - EE - BB - TE
//sprintf(XSpectrumFile,"%s_%d_%d.fits",pathToXSpectra.c_str(),_XSpectra2Maps[c][0],_XSpectra2Maps[c][1]);
//planck_assert(file_present(XSpectrumFile),string("missing file : ")+XSpectrumFile);
ostringstream os ;
os << pathToXSpectra << "_" << _XSpectra2Maps[c][0] << "_" << _XSpectra2Maps[c][1] << ".fits";
string XSpectrumFile=Parser::CheckPath(os.str());
......@@ -321,7 +321,7 @@ void HiLLiPOP::ProcessXSpectraErrors(const string pathToXSpectraErrors)
//sprintf(XSpectrumFile,"%s_%d_%d.fits",pathToXSpectraErrors.c_str(),_XSpectra2Maps[c][0],_XSpectra2Maps[c][1]);
//planck_assert(file_present(XSpectrumFile),string("missing file : ")+XSpectrumFile);
ostringstream os ;
os << pathToXSpectraErrors << "_" << _XSpectra2Maps[c][0] << "_" << _XSpectra2Maps[c][1] << ".fits";
string XSpectrumFile=Parser::CheckPath(os.str());
......@@ -333,7 +333,7 @@ void HiLLiPOP::ProcessXSpectraErrors(const string pathToXSpectraErrors)
// ET
//sprintf(XSpectrumFile,"%s_%d_%d.fits",pathToXSpectraErrors.c_str(),_XSpectra2Maps[c][1],_XSpectra2Maps[c][0]);
//planck_assert(file_present(XSpectrumFile),string("missing file : ")+XSpectrumFile);
os.clear();os.str("");
os << pathToXSpectraErrors << "_" << _XSpectra2Maps[c][1] << "_" << _XSpectra2Maps[c][0] << ".fits";
XSpectrumFile=Parser::CheckPath(os.str());
......@@ -362,7 +362,7 @@ void HiLLiPOP::ProcessXSpectraErrors(const string pathToXSpectraErrors)
if( ell[size-1]-_lmaxXSpectra[type][c] < 0 ) status = false;
planck_assert(status,string("multipoles mismatch"));
for(unsigned int l = 0; l <= _maxOflmax[0]-ell[0]; l++) _ClWeightData[c][m*(_maxOflmax[0]+1)+ell[l]] = pow(2*M_PI/ell[l]/(ell[l]+1)*ClErr[l]*1e12,-2);
for(unsigned int l = 0; l <= _maxOflmax[0]-ell[0]; l++) _ClWeightData[c][m*(_maxOflmax[0]+1)+ell[l]] = pow(2*M_PI/ell[l]/(ell[l]+1)*ClErr[l]*1e12,-2);
delete [] ell;
delete [] ClErr;
......@@ -437,7 +437,7 @@ 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++)
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;
}
}
......@@ -498,15 +498,15 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC
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
if( _typeStatus[0]) {
planck_assert(file_present(SZFile),string("missing file : ")+SZFile);
fitshandle * SZInput = new fitshandle();
SZInput->open(SZFile);
SZInput->goto_hdu(2);
unsigned int sizeSZ = SZInput->nelems(2);
double * ellSZ = new double[sizeSZ];
double * Cl1haloSZ = new double[sizeSZ];
......@@ -514,20 +514,20 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC
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];
}
SZInput->close();
delete [] ellSZ;
delete [] Cl1haloSZ;
delete [] Cl2haloSZ;
delete SZInput;
_n.push_back("Asz");
}
// End : tSZ
......@@ -551,15 +551,15 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC
for(unsigned int f = 0; f < _nXFreq; f++) {
//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";
tmpFile=Parser::CheckPath(os.str());
fitshandle * CIBInput = new fitshandle();
CIBInput->open(tmpFile);
CIBInput->goto_hdu(2);
unsigned int sizeCIB = CIBInput->nelems(2);
double * ellCIB = new double[sizeCIB];
double * Cl1haloCIB = new double[sizeCIB];
......@@ -567,12 +567,12 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC
CIBInput->read_column_raw(1,ellCIB,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];
CIBInput->close();
delete [] ellCIB;
delete [] Cl1haloCIB;
delete [] Cl2haloCIB;
......@@ -606,7 +606,7 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC
for(unsigned int f = 0; f < _nXFreq; f++) {
//sprintf(tmpFile,"%s_%s.fits",pathToDustSpectra.c_str(),_XFreq[f].c_str());
//planck_assert(file_present(tmpFile),string("missing file : ")+tmpFile);
ostringstream os ;
os << pathToDustSpectra << "_" << _XFreq[f] << ".fits";
tmpFile=Parser::CheckPath(os.str());
......@@ -644,27 +644,27 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC
// Begin : kSZ
if( _typeStatus[0]) {
planck_assert(file_present(kSZFile),string("missing file : ")+kSZFile);
fitshandle * kSZInput = new fitshandle();
kSZInput->open(kSZFile);
kSZInput->goto_hdu(2);
//kSZ template, Dl=l(l+1)/2pi Cl, units uK
unsigned int sizekSZ = kSZInput->nelems(2);
double * ellkSZ = new double[sizekSZ];
double * ClkSZ = new double[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];
kSZInput->close();
delete [] ellkSZ;
delete [] ClkSZ;
delete kSZInput;
_n.push_back("Aksz");
}
// End : kSZ
......@@ -676,32 +676,32 @@ void HiLLiPOP::ProcessNuisanceVariables(const string SZFile,const string pathToC
for(unsigned int f = 0; f < _nXFreq; f++) {
//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";
tmpFile=Parser::CheckPath(os.str());
fitshandle * SZxCIBInput = new fitshandle();
SZxCIBInput->open(tmpFile);
SZxCIBInput->goto_hdu(2);
unsigned int sizeSZxCIB = SZxCIBInput->nelems(2);
double * ellSZxCIB = new double[sizeSZxCIB];
double * ClSZxCIB = new double[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];
SZxCIBInput->close();
delete [] ellSZxCIB;
delete [] ClSZxCIB;
delete SZxCIBInput;
}
_n.push_back("Aszxcib");
}
// End : SZxCIB
......@@ -954,8 +954,8 @@ void HiLLiPOP::computeResiduals( vector<double> ClCMB, vector<double>& nuisance,
unsigned int modeOffset = 0;
char tmpPar[32];
vector<double> Aps(_nXFreq);
double Asz, Acib, Aksz, Aszxcib;
double AdustTT, AdustTP, AdustPP;
double Asz = 0.0, Acib = 0.0, Aksz = 0.0, Aszxcib = 0.0;
double AdustTT = 0.0, AdustTP = 0.0, AdustPP = 0.0;
_residual.resize(5);
_instrumental.resize(5);
......@@ -1041,14 +1041,14 @@ void HiLLiPOP::computeResiduals( vector<double> ClCMB, vector<double>& nuisance,
lminXFreq = _lminXFreq[m][f];
lmaxXFreq = _lmaxXFreq[m][f];
}
for(unsigned int l = lminXFreq; l <= lmaxXFreq; l++) {
if( m == 0 ) foregrounds = Aps[f]+Asz*_ClSZ[f][l]+Acib*_ClCIB[f][l]+AdustTT*_ClDust[f][modeOffset+l]+Aksz*_ClkSZ[l]+Aszxcib*_ClSZxCIB[f][l];
if( m == 1 || m == 2 ) foregrounds = AdustPP*_ClDust[f][modeOffset+l];
if( m == 3 || m == 4 ) foregrounds = AdustTP*_ClDust[f][modeOffset+l];
for(unsigned int c = 0; c < _XFreq2XSpectra[f].size(); c++) {
if( !fullRange) {
if( !fullRange) {
lminXSpectra = _lminXSpectra[m][_XFreq2XSpectra[f][c]];
lmaxXSpectra = _lmaxXSpectra[m][_XFreq2XSpectra[f][c]];
}
......@@ -1072,7 +1072,7 @@ void HiLLiPOP::computeResiduals( vector<double> ClCMB, vector<double>& nuisance,
if( m == 3 ) {
tmpPolEffTP = 1+epsilon[_XSpectra2Maps[_XFreq2XSpectra[f][c]][1]];
_residual[m][f][l] += tmpWeight*(_ClData[_XFreq2XSpectra[f][c]][modeOffset+l]-tmpCal*tmpPolEffTP*tmpBeamEigenmodes*(ClCMB[modeOffset+l]+foregrounds));
_instrumental[m][f][l] += tmpWeight*tmpCal*tmpPolEffTP*tmpBeamEigenmodes;
_instrumental[m][f][l] += tmpWeight*tmpCal*tmpPolEffTP*tmpBeamEigenmodes;
}
if( m == 4 ) {
tmpPolEffPT = 1+epsilon[_XSpectra2Maps[_XFreq2XSpectra[f][c]][0]];
......@@ -1227,7 +1227,7 @@ void HiLLiPOP::WriteOutput(const vector<double>& ClCMB,const vector<double>& nui
tmpWeight = _ClWeightData[_XFreq2XSpectra[f][c]][modeOffset+l];
tmpBeamEigenmodes = pow(1+beta[_XFreq2XSpectra[f][c]]*_beamEigenmodes[_XFreq2XSpectra[f][c]][l],2);
tmpCal = 1+cal[_XSpectra2Maps[_XFreq2XSpectra[f][c]][0]]+cal[_XSpectra2Maps[_XFreq2XSpectra[f][c]][1]];
#ifdef PolarizationEfficiency
#ifdef PolarizationEfficiency
if( m == 0 ) {
residual[m][f][l] += tmpWeight*(_ClData[_XFreq2XSpectra[f][c]][modeOffset+l]-tmpCal*tmpBeamEigenmodes*(ClCMB[modeOffset+l]+foregrounds));
instrumental[m][f][l] += tmpWeight*tmpCal*tmpBeamEigenmodes;
......@@ -1240,7 +1240,7 @@ void HiLLiPOP::WriteOutput(const vector<double>& ClCMB,const vector<double>& nui
if( m == 3 ) {
tmpPolEffTP = 1+epsilon[_XSpectra2Maps[_XFreq2XSpectra[f][c]][1]];
residual[m][f][l] += tmpWeight*(_ClData[_XFreq2XSpectra[f][c]][modeOffset+l]-tmpCal*tmpPolEffTP*tmpBeamEigenmodes*(ClCMB[modeOffset+l]+foregrounds));
instrumental[m][f][l] += tmpWeight*tmpCal*tmpPolEffTP*tmpBeamEigenmodes;
instrumental[m][f][l] += tmpWeight*tmpCal*tmpPolEffTP*tmpBeamEigenmodes;
}
if( m == 4 ) {
tmpPolEffPT = 1+epsilon[_XSpectra2Maps[_XFreq2XSpectra[f][c]][0]];
......@@ -1305,7 +1305,7 @@ void HiLLiPOP::WriteOutput(const vector<double>& ClCMB,const vector<double>& nui
of << "#Chi2 = " << _chi2 << endl;
of << "#Multipole";
for(unsigned int m = 0; m < 4; m++) of << "\t CMB" << CMBmode[m];
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])))
for(unsigned int f = 0; f < _nXFreq; f++) {
of << "\t Residual" << CMBmode[m] << _XFreq[f];
......@@ -1318,7 +1318,7 @@ void HiLLiPOP::WriteOutput(const vector<double>& ClCMB,const vector<double>& nui
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])))
for(unsigned int f = 0; f < _nXFreq; f++) of << "\t Dust" << CMBmode[m] << _XFreq[f];
of << endl;
......@@ -1330,7 +1330,7 @@ void HiLLiPOP::WriteOutput(const vector<double>& ClCMB,const vector<double>& nui
for(unsigned int m = 0; m < 4; m++) {
if( _modeStatus[m] || (m==3 && (_modeStatus[3] || _modeStatus[4])))
for(unsigned int f = 0; f < _nXFreq; f++) {
if( _XFreqStatus[m][f] ) of << " " << _residual[m][f][l]; // << " " << _instrumental[m][f][l];
if( _XFreqStatus[m][f] ) of << " " << _residual[m][f][l]; // << " " << _instrumental[m][f][l];
else of << " " << 0;// << " " << 0;
}
}
......@@ -1359,4 +1359,3 @@ void HiLLiPOP::WriteOutput(const vector<double>& ClCMB,const vector<double>& nui
of.close();
// End: output
}
src/camel/CMB/HighEll.cc
View file @ 57b8f7d1
This diff is collapsed.
src/camel/CMB/Hillipop.cc
View file @ 57b8f7d1
......@@ -37,7 +37,7 @@ void Hillipop::Init(const string fileName)
cout << "Using: " << fileName << endl;
paramfile parameters(fileName);
string multipolesFile = Parser::CheckPath(parameters.find<string>("MultipolesRange"));
string SZFile = Parser::CheckPath(parameters.find<string>("SZ"));
string kSZFile = Parser::CheckPath(parameters.find<string>("kSZ"));
......@@ -66,7 +66,7 @@ void Hillipop::Init(const string fileName)
_typeStatus[0] = _modeStatus[0];
_typeStatus[1] = _modeStatus[1];
_typeStatus[2] = _modeStatus[2];
if( _modeStatus[3] == 1 || _modeStatus[4] == 1) _typeStatus[3] = 1; else _typeStatus[3] = 0;
if( _modeStatus[3] == 1 || _modeStatus[4] == 1) _typeStatus[3] = 1; else _typeStatus[3] = 0;
ProduceList();
......@@ -141,7 +141,7 @@ void Hillipop::ProcessMultipoles(const string multipolesFile)
for(unsigned int m = 0; m < 5; m++) {
if( m < 4) type = m; else type = 3;
input.goto_hdu(type+2);
input.goto_hdu(type+2);
_lminXSpectra[m].resize(_nXSpectra);
_lmaxXSpectra[m].resize(_nXSpectra);
......@@ -258,7 +258,7 @@ void Hillipop::ProcessXSpectra(const string pathToXSpectra)
// TT - EE - BB - TE
//sprintf(XSpectrumFile,"%s_%d_%d.fits",pathToXSpectra.c_str(),_XSpectra2Maps[c][0],_XSpectra2Maps[c][1]);
//planck_assert(file_present(XSpectrumFile),string("missing file : ")+XSpectrumFile);
ostringstream os ;
os << pathToXSpectra << "_" << _XSpectra2Maps[c][0] << "_" << _XSpectra2Maps[c][1] << ".fits";
string XSpectrumFile=Parser::CheckPath(os.str());
......@@ -323,7 +323,7 @@ void Hillipop::ProcessXSpectraErrors(const string pathToXSpectraErrors)
//sprintf(XSpectrumFile,"%s_%d_%d.fits",pathToXSpectraErrors.c_str(),_XSpectra2Maps[c][0],_XSpectra2Maps[c][1]);
//planck_assert(file_present(XSpectrumFile),string("missing file : ")+XSpectrumFile);
ostringstream os ;
os << pathToXSpectraErrors << "_" << _XSpectra2Maps[c][0] << "_" << _XSpectra2Maps[c][1] << ".fits";
string XSpectrumFile=Parser::CheckPath(os.str());
......@@ -335,7 +335,7 @@ void Hillipop::ProcessXSpectraErrors(const string pathToXSpectraErrors)
// ET
//sprintf(XSpectrumFile,"%s_%d_%d.fits",pathToXSpectraErrors.c_str(),_XSpectra2Maps[c][1],_XSpectra2Maps[c][0]);
//planck_assert(file_present(XSpectrumFile),string("missing file : ")+XSpectrumFile);
os.clear();os.str("");
os << pathToXSpectraErrors << "_" << _XSpectra2Maps[c][1] << "_" << _XSpectra2Maps[c][0] << ".fits";
XSpectrumFile=Parser::CheckPath(os.str());
......@@ -364,7 +364,7 @@ void Hillipop::ProcessXSpectraErrors(const string pathToXSpectraErrors)
if( ell[size-1]-_lmaxXSpectra[type][c] < 0 ) status = false;
planck_assert(status,string("multipoles mismatch"));
for(unsigned int l = 0; l <= _maxOflmax[0]-ell[0]; l++) _ClWeightData[c][m*(_maxOflmax[0]+1)+ell[l]] = pow(2*M_PI/ell[l]/(ell[l]+1)*ClErr[l]*1e12,-2);
for(unsigned int l = 0; l <= _maxOflmax[0]-ell[0]; l++) _ClWeightData[c][m*(_maxOflmax[0]+1)+ell[l]] = pow(2*M_PI/ell[l]/(ell[l]+1)*ClErr[l]*1e12,-2);
delete [] ell;
delete [] ClErr;
......@@ -525,14 +525,14 @@ void Hillipop::ProcessNuisanceVariables(const string SZFile,const string pathToC
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
planck_assert(file_present(SZFile),string("missing file : ")+SZFile);
fitshandle * SZInput = new fitshandle();
SZInput->open(SZFile);
SZInput->goto_hdu(2);
unsigned int sizeSZ = SZInput->nelems(2);
double * ellSZ = new double[sizeSZ];
double * Cl1haloSZ = new double[sizeSZ];
......@@ -540,7 +540,7 @@ void Hillipop::ProcessNuisanceVariables(const string SZFile,const string pathToC
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);
......@@ -548,7 +548,7 @@ void Hillipop::ProcessNuisanceVariables(const string SZFile,const string pathToC
}
SZInput->close();
delete [] ellSZ;
delete [] Cl1haloSZ;
delete [] Cl2haloSZ;
......@@ -565,15 +565,15 @@ void Hillipop::ProcessNuisanceVariables(const string SZFile,const string pathToC
for(unsigned int f = 0; f < _nXFreq; f++) {
//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";
tmpFile=Parser::CheckPath(os.str());
fitshandle * CIBInput = new fitshandle();
CIBInput->open(tmpFile);
CIBInput->goto_hdu(2);
unsigned int sizeCIB = CIBInput->nelems(2);
double * ellCIB = new double[sizeCIB];
double * Cl1haloCIB = new double[sizeCIB];
......@@ -581,18 +581,18 @@ void Hillipop::ProcessNuisanceVariables(const string SZFile,const string pathToC
CIBInput->read_column_raw(1,ellCIB,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];
CIBInput->close();
delete [] ellCIB;
delete [] Cl1haloCIB;
delete [] Cl2haloCIB;
delete CIBInput;
}
_n.push_back("Acib");
}
// End : CIB
......@@ -620,7 +620,7 @@ void Hillipop::ProcessNuisanceVariables(const string SZFile,const string pathToC
for(unsigned int f = 0; f < _nXFreq; f++) {
//sprintf(tmpFile,"%s_%s.fits",pathToDustSpectra.c_str(),_XFreq[f].c_str());
//planck_assert(file_present(tmpFile),string("missing file : ")+tmpFile);
ostringstream os ;
os << pathToDustSpectra << "_" << _XFreq[f] << ".fits";
tmpFile=Parser::CheckPath(os.str());
......@@ -669,12 +669,12 @@ void Hillipop::ProcessNuisanceVariables(const string SZFile,const string pathToC
double * ClkSZ = new double[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];
kSZInput->close();
delete [] ellkSZ;
delete [] ClkSZ;
delete kSZInput;
......@@ -690,27 +690,27 @@ void Hillipop::ProcessNuisanceVariables(const string SZFile,const string pathToC
for(unsigned int f = 0; f < _nXFreq; f++) {
//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";
tmpFile=Parser::CheckPath(os.str());
fitshandle * SZxCIBInput = new fitshandle();
SZxCIBInput->open(tmpFile);
SZxCIBInput->goto_hdu(2);
unsigned int sizeSZxCIB = SZxCIBInput->nelems(2);
double * ellSZxCIB = new double[sizeSZxCIB];
double * ClSZxCIB = new double[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];
SZxCIBInput->close();
delete [] ellSZxCIB;
delete [] ClSZxCIB;
delete SZxCIBInput;
......@@ -862,7 +862,7 @@ void Hillipop::WriteOutput(const vector<double>& ClCMB,const vector<double>& nui