diff --git a/src/AbsRand/OofNoise.cc b/src/AbsRand/OofNoise.cc index c18aebb0f7da7d16490b4369e11a60da433f2a21..ada73c92c506418f6a6f13a116ac77e76f814be3 100644 --- a/src/AbsRand/OofNoise.cc +++ b/src/AbsRand/OofNoise.cc @@ -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 p(Nproc); - vector z(Nproc); - - double dp=(Wmax-Wmin)/Nproc; - p[0]=Wmin+(1-a/2)/2*dp; - z[0]=p[0]+a/2*dp; + vector p(Nproc); + vector z(Nproc); - for(int i=1;i pole(Nproc); - vector zero(Nproc); - for (int i=0;i pole(Nproc); + vector zero(Nproc); + for (int i=0;inormal(); for (unsigned i=0;ifilter(x2); return x2*_sigma; -}; - +} diff --git a/src/camel/CMB/HiLLiPOP.cc b/src/camel/CMB/HiLLiPOP.cc index e79dc8bc32f3320bd4ea08d5cca2fa33955dc852..03be755f1a3822367329e896d798a7ec57c06bab 100644 --- a/src/camel/CMB/HiLLiPOP.cc +++ b/src/camel/CMB/HiLLiPOP.cc @@ -38,7 +38,7 @@ void HiLLiPOP::Init(const string fileName) cout << "Using: " << fileName << endl; paramfile parameters(fileName); - + string multipolesFile = Parser::CheckPath(parameters.find("MultipolesRange")); string SZFile = Parser::CheckPath(parameters.find("SZ")); string kSZFile = Parser::CheckPath(parameters.find("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 ClCMB, vector& nuisance, unsigned int modeOffset = 0; char tmpPar[32]; vector 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 ClCMB, vector& 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 ClCMB, vector& 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& ClCMB,const vector& 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& ClCMB,const vector& 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& ClCMB,const vector& 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& ClCMB,const vector& 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& ClCMB,const vector& 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& ClCMB,const vector& nui of.close(); // End: output } - diff --git a/src/camel/CMB/HighEll.cc b/src/camel/CMB/HighEll.cc index 8b6a032123698c1859b6d2073d22e094012c1942..03d892e6968d4185cb972d59a18f6b3cb0b5e39c 100644 --- a/src/camel/CMB/HighEll.cc +++ b/src/camel/CMB/HighEll.cc @@ -17,25 +17,25 @@ HighEll::HighEll(const string parfile ) { _plik_mode = false; read_par_file( parfile ); - // nuisance names for all + // nuisance names for all if (_plik_mode){ - // params for plik + // params for plik _n.push_back("A_sz"); // tSZ - _n.push_back("ksz_norm"); // kSZ - _n.push_back("A_cib_217"); // CIB - _n.push_back("xi_sz_cib"); + _n.push_back("ksz_norm"); // kSZ + _n.push_back("A_cib_217"); // CIB + _n.push_back("xi_sz_cib"); }else{ - _n.push_back("Asz");// tSZ + _n.push_back("Asz");// tSZ _n.push_back("Aksz"); - _n.push_back("Acib"); - _n.push_back("Aszxcib"); + _n.push_back("Acib"); + _n.push_back("Aszxcib"); } _numNuisAll = 4; } - - + + HighEll::~HighEll(){ } @@ -45,7 +45,7 @@ void HighEll:: read_par_file(const string filename ) paramfile parameters(filename); _parfile = filename ; - //(Parser::CheckPath(.find("DataSet",""))).c_str(); + //(Parser::CheckPath(.find("DataSet",""))).c_str(); //_verbose = 0; _verbose = parameters.find("verbose"); @@ -58,17 +58,17 @@ void HighEll:: read_par_file(const string filename ) _nfreqs = parameters.find("nb_freqs"); for(unsigned int i = 0; i < _nfreqs ; i++) { - + sprintf(tmp,"freq%d",i); _freq.push_back(parameters.find(tmp)); } _name = parameters.find("Name") ; - _nBins = parameters.find("nBins"); // nbre de bins ds la 'theorie - // min max en ell pour les theories (fg et CMB) + _nBins = parameters.find("nBins"); // nbre de bins ds la 'theorie + // min max en ell pour les theories (fg et CMB) _lmin = parameters.find("ell_min"); _lmax = parameters.find("ell_max"); _lmaxCMB = parameters.find("ell_max_CMB"); - + _size = 0; _tSZ_file = Parser::CheckPath(parameters.find("tSZ_file")); @@ -91,13 +91,13 @@ void HighEll:: read_par_file(const string filename ) _chi2_prec = 0.; cout << " verbose "<< _verbose<0) { for (int ja=0; ja< _nparts ; ja++){ - int jbmin = ja; + int jbmin = ja; if (ib != ia) jbmin = 0; for (int jb=jbmin; jb<_nparts ; jb++){ - + name = Parser::CheckPath(get_blname( _freq[ia], _freq[ib], _part[ja], _part[jb])); - + cout << " name of file "<< name<<" "<< _lmax << endl; - + int size = _nBins; if (ib ==0) size = _ibmax[0]; - + HepMatrix winf( size , _lmax-2); //TOCHECK ifstream file(name.c_str(),ios::in); double ell,x; - + for (int l=2 ; l<_lmax ; l++){ file >> ell; //cout <<" at "<> x; //cout << x ; - winf[k][l-2] = x; - } + winf[k][l-2] = x; + } // cout << endl; } _winFct.push_back( winf ); - + file.close(); - + } } }else{ @@ -153,18 +153,18 @@ void HighEll::read_window_functions(){ HepMatrix winf( size , _lmax-2); //TOCHECK ifstream file(name.c_str(),ios::in); double ell,x; - + for (int l=2 ; l<_lmax ; l++){ file >> ell; - // cout << " read "<< ell<> x; //cout << " read "<0 ) if ( _select_freq[ib] !=0) continue; if (_nparts >0) { - + for (int ja=0; ja< _nparts ; ja++){ - int jbmin = ja; + int jbmin = ja; if (ib != ia) jbmin = 0; for (int jb=jbmin; jb<_nparts ; jb++){ - - string name; - + + string name; + name = Parser::CheckPath(get_spname( _freq[ia], _freq[ib], _part[ja], _part[jb])); cout << " reading "< tempv; vector tempe; vector tempev; int ityp = getType( _freq[ia], _freq[ib]); - + for (int k=0 ; k< _ibmax[ityp]- _ibmin[ityp] ; k++){ file >> a >> b >> c; - tempv.push_back(b); - tempe.push_back(a); - tempev.push_back(c); + tempv.push_back(b); + tempe.push_back(a); + tempev.push_back(c); } _dataLength += _ibmax[ityp]- _ibmin[ityp]; - _size += _ibmax[ityp]- _ibmin[ityp]; // doubleon ?? + _size += _ibmax[ityp]- _ibmin[ityp]; // doubleon ?? file.close(); _inSpectra.push_back(tempv); _inSpectraErrors.push_back(tempev); @@ -222,7 +222,7 @@ void HighEll::read_spectra_data() } } } else { - string name; + string name; string toto; name = Parser::CheckPath(get_spname( _freq[ia], _freq[ib],toto,toto)); cout << " reading "<> a >> b >> c; - tempv.push_back(b); - tempe.push_back(a); - tempev.push_back(c); + tempv.push_back(b); + tempe.push_back(a); + tempev.push_back(c); } _dataLength += _ibmax[ityp]- _ibmin[ityp]; - _size += _ibmax[ityp]- _ibmin[ityp]; // doubleon ?? - + _size += _ibmax[ityp]- _ibmin[ityp]; // doubleon ?? + file.close(); cout <<" read "<< tempv.size()<<" "<& l,vector& cltt,vector& clte,vector& clee,vector& clbb,vector& nuisance){ @@ -277,7 +277,7 @@ void HighEll::computeResiduals(const vector& l,vector& clt vector < vector > myDust = getDust( ); cout << " got dust template"<& l,vector& clt fileout = new ofstream (fname.c_str(),ios::out | ios::trunc); for (int k=0 ; k< _nSpectra ; k++) { - HepVector clth(_lmax-2); + HepVector clth(_lmax-2); for (int i=0; i<_lmax-2 ; i++) { clth[i] = 0.; if ( i+2< _lmaxCMB){ @@ -300,8 +300,8 @@ void HighEll::computeResiduals(const vector& l,vector& clt //cout<<_winFct[k]<1)cout << "spec no "< calib data and covmat + + clbin = _winFct[k]*clth / _calibs[_types[k]]; // calibrates model <=> calib data and covmat for (int ii = 0 ; ii< clbin.num_row() ; ii++){ if (_verbose >1) cout << " CMB spect "<& l,vector& clt fileout = new ofstream (fname.c_str(),ios::out | ios::trunc); for (int k=0 ; k< _nSpectra ; k++) { - HepVector clth(_lmax-2); + HepVector clth(_lmax-2); for (int i=0; i<_lmax-2 ; i++) { clth[i] = 0.; - clth[i] = myCIB[ _types[k]][i] ;// + clth[i] = myCIB[ _types[k]][i] ;// if ( _ellFactor != 0) clth[i] *= 2.*M_PI/(i+2.)/(i+3.) ; //cout << " components : cib "<< myCIB[ _types[k]][i] << " tsz "<< mytSZ[ _types[k]][i] <<" dust "<< myDust[ _types[k]][i] << " "<< _aPS[_types[k]] << " tot "<< clth[i] << endl; } HepVector clbin; cout << "spec no "< calib data and covmat + clbin = _winFct[k]*clth / _calibs[_types[k]]; // calibrates model <=> calib data and covmat for ( int ii = 0 ; ii< clbin.num_row() ; ii++) { if (_verbose >1)cout << " CIB spect "<close(); delete fileout; @@ -351,21 +351,21 @@ void HighEll::computeResiduals(const vector& l,vector& clt fileout = new ofstream (fname.c_str(),ios::out | ios::trunc); for (int k=0 ; k< _nSpectra ; k++) { - HepVector clth(_lmax-2); + HepVector clth(_lmax-2); for (int i=0; i<_lmax-2 ; i++) { clth[i] = 0.; - clth[i] = mytSZxCIB[ _types[k]][i] ;// + clth[i] = mytSZxCIB[ _types[k]][i] ;// if ( _ellFactor != 0) clth[i] *= 2.*M_PI/(i+2.)/(i+3.) ; //cout << " components : cib "<< myCIB[ _types[k]][i] << " tsz "<< mytSZ[ _types[k]][i] <<" dust "<< myDust[ _types[k]][i] << " "<< _aPS[_types[k]] << " tot "<< clth[i] << endl; } HepVector clbin; cout << "spec no "< calib data and covmat + clbin = _winFct[k]*clth / _calibs[_types[k]]; // calibrates model <=> calib data and covmat for ( int ii = 0 ; ii< clbin.num_row() ; ii++) { if (_verbose >1)cout << " tSZxCIB spect "<close(); delete fileout; @@ -378,7 +378,7 @@ void HighEll::computeResiduals(const vector& l,vector& clt fileout = new ofstream (fname.c_str(),ios::out | ios::trunc); for (int k=0 ; k< _nSpectra ; k++) { - HepVector clth(_lmax-2); + HepVector clth(_lmax-2); for (int i=0; i<_lmax-2 ; i++) { clth[i] = 0.; clth[i] = myDust[ _types[k]][i] ; @@ -387,7 +387,7 @@ void HighEll::computeResiduals(const vector& l,vector& clt } HepVector clbin; if (_verbose >1)cout << "spec no "< calib data and covmat + clbin = _winFct[k]*clth / _calibs[_types[k]]; // calibrates model <=> calib data and covmat for ( int ii = 0 ; ii< clbin.num_row() ; ii++) { cout << " Dust spect "<& l,vector& clt delete fileout; if (_verbose >1)cout << " end loop on spectra for Dust "<1)cout << " tSZ: begin loop on spectra "<& l,vector& clt fileout = new ofstream (fname.c_str(),ios::out | ios::trunc); for (int k=0 ; k< _nSpectra ; k++) { - HepVector clth(_lmax-2); - for (int i=0; i<_lmax-2 ; i++) { + HepVector clth(_lmax-2); + for (int i=0; i<_lmax-2 ; i++) { clth[i] = mytSZ[ _types[k]][i] ; //cout << " components : cib "<< myCIB[ _types[k]][i] << " tsz "<< mytSZ[ _types[k]][i] <<" dust "<< myDust[ _types[k]][i] << " "<< _aPS[_types[k]] << " tot "<< clth[i] << endl; if ( _ellFactor != 0) clth[i] *= 2.*M_PI/(i+2.)/(i+3.) ; @@ -415,7 +415,7 @@ void HighEll::computeResiduals(const vector& l,vector& clt //cout<<_winFct[k]<1)cout << "spec no "< calib data and covmat + clbin = _winFct[k]*clth / _calibs[_types[k]]; // calibrates model <=> calib data and covmat for (int ii = 0 ; ii< clbin.num_row() ; ii++){ cout << " SZ spect "<& l,vector& clt fileout = new ofstream (fname.c_str(),ios::out | ios::trunc); for (int k=0 ; k< _nSpectra ; k++) { - HepVector clth(_lmax-2); + HepVector clth(_lmax-2); for (int i=0; i<_lmax-2 ; i++) { - + clth[i] = mykSZ[ _types[k]][i] ; //cout << " components : cib "<< myCIB[ _types[k]][i] << " tsz "<< mytSZ[ _types[k]][i] <<" dust "<< myDust[ _types[k]][i] << " "<< _aPS[_types[k]] << " tot "<< clth[i] << endl; if ( _ellFactor != 0) clth[i] *= 2.*M_PI/(i+2.)/(i+3.) ; @@ -444,7 +444,7 @@ void HighEll::computeResiduals(const vector& l,vector& clt //cout<<_winFct[k]< calib data and covmat + clbin = _winFct[k]*clth / _calibs[_types[k]]; // calibrates model <=> calib data and covmat for (int ii = 0 ; ii< clbin.num_row() ; ii++){ if (_verbose >1) cout << " kSZ spect "<& l,vector& clt delete fileout; if (_verbose >1)cout << " end loop on spectra for kSZ "<& l,vector& clt fileout = new ofstream (fname.c_str(),ios::out | ios::trunc); for (int k=0 ; k< _nSpectra ; k++) { - HepVector clth(_lmax-2); + HepVector clth(_lmax-2); for (int i=0; i<_lmax-2 ; i++) { - + clth[i] = _aPS[_types[k]] * pow(((i+2.)/3000.),2.) ; //cout << " components : cib "<< myCIB[ _types[k]][i] << " tsz "<< mytSZ[ _types[k]][i] <<" dust "<< myDust[ _types[k]][i] << " "<< _aPS[_types[k]] << " tot "<< clth[i] << endl; if ( _ellFactor != 0) clth[i] *= 2.*M_PI/(i+2.)/(i+3.) ; @@ -476,7 +476,7 @@ void HighEll::computeResiduals(const vector& l,vector& clt //cout<<_winFct[k]<1)cout << "spec no "< calib data and covmat + clbin = _winFct[k]*clth / _calibs[_types[k]]; // calibrates model <=> calib data and covmat for (int ii = 0 ; ii< clbin.num_row() ; ii++){ cout << "PS spect "<& l,vector& clt delete fileout; if (_verbose >1) cout << " end loop on spectra for PS "<& l,vector& if (_verbose >2) cout << " high ell likelihood compute for"<< _name<<" nuisance size "<< nuisance.size()<2) cout << " handled nuisance params for "<< _name<2) { string fname = "outLike_"; fname = fname + _name; @@ -520,7 +520,7 @@ double HighEll::computeLikelihood(const vector& l,vector& if (_verbose >2)cout <<" win fct sizew"<< _winFct.size() << endl; for (int k=0 ; k< _nSpectra ; k++) { //cout << " spectrum "<& l,vector& clth[i] = clth[i] + myDust[ _types[k]][i] ; clth[i] = clth[i] + myCIB[ _types[k]][i] ; clth[i] = clth[i] + mytSZxCIB[ _types[k]][i] ; - clth[i] = clth[i] + _aPS[_types[k]] * pow(((i+2.)/3000.),2.); // PS + clth[i] = clth[i] + _aPS[_types[k]] * pow(((i+2.)/3000.),2.); // PS if (i+2< _lmaxCMB) clth[i] = clth[i] + cltt[i+2]*(i+2.)*(i+3.)/2./M_PI; // CMB [uK^2 ?] if ( clth[i] < min) min = clth[i]; if ( clth[i] > max) max = clth[i]; - if ( _ellFactor != 0) clth[i] *= 2.*M_PI/(i+2.)/(i+3.) ; // for ACT + if ( _ellFactor != 0) clth[i] *= 2.*M_PI/(i+2.)/(i+3.) ; // for ACT } HepVector clbin; - clbin = _winFct[k]*clth / cali ; // calibrates model <=> calib data and covmat + clbin = _winFct[k]*clth / cali ; // calibrates model <=> calib data and covmat clSky_bin.push_back(clbin); } @@ -552,33 +552,35 @@ double HighEll::computeLikelihood(const vector& l,vector& if ( _verbose > 2 ) (*fileout) << "# type num bin binTh model data ell err clchi"< 2 ) (*fileout) << _types[k] <<" "< 2 ) cout <<" bin "<< k<<" clchi "<< clChi[k]< 2 ) cout <<" bin "<< k<<" clchi "<< clChi[k]<2)cout << "chi2 computation ndof = "<< clChi.num_row()<< " matrix dims "<< _cov_matrix->num_col() <<"x"<< _cov_matrix->num_row()<< " "<< _size<<" "<< _dataLength<< endl; mychi2 = _cov_matrix->similarity(clChi); - + _chi2_prec = mychi2; - if (_verbose >1)cout << " High ell likelihood for "<< _name<<" chi2 " << mychi2 <<" ndof "<< _size < 2 ) fileout->close(); - if ( _verbose > 2 ) delete fileout; - if ( _verbose > 2 ) computeResiduals(l, cltt, clte, clee, clbb, nuisance); + if (_verbose >1) cout << " High ell likelihood for "<< _name<<" chi2 " << mychi2 <<" ndof "<< _size < 2) { + fileout->close(); + delete fileout; + computeResiduals(l, cltt, clte, clee, clbb, nuisance); + } _lnlkl= mychi2/2.; - + return(mychi2/2.); } @@ -600,7 +602,7 @@ double HighEll::planck_func_ratio(double freq,double fref,double t_eff){ double nu0=fref*1.e9; double xx=h_pl*nu/(k_b*t_eff); double xx0=h_pl*nu0/(k_b*t_eff); - + corr = ((nu/nu0)*(nu/nu0)*(nu/nu0))*(exp(xx0)-1.)/(exp(xx)-1.); return(corr); @@ -612,7 +614,7 @@ double HighEll::flxtotemp(double fa,double fref){ double nu0=fref*1.e9; double xx=h_pl*nu/(k_b*t_cmb); double xx0=h_pl*nu0/(k_b*t_cmb); - + corr = pow(nu0/nu,4.)*(exp(xx0)/exp(xx))*pow(((exp(xx)-1.)/(exp(xx0)-1.)),2) ; return(corr); @@ -641,7 +643,7 @@ void HighEll::readAllCIBTemplates(int lmin, int lmax) { vector vbuf; int nmax_read = 13000; for (int i= 0; (i> a >> b1 >> b2 >> b3 >> b4 >> b5 >> b6>>b7>>b8>>b9>> b10>> b11 >> b12 >> b13 >> b14 >> b15 >> b16>>b17>>b18>>b19>>b20>>b21 ; double spe = 0.; @@ -653,16 +655,16 @@ void HighEll::readAllCIBTemplates(int lmin, int lmax) { if ( _idfre[k] == 2 && _idfre[j] == 2) spe = b12; //if ( a < 100.) cout << _idfre[k]<< " "<< _idfre[j] <<" a= "<=lmin && ( (int) a <=lmax ) ) { - vbuf.push_back(spe* a*(a+1.)/2./M_PI/facta/factb); + vbuf.push_back(spe* a*(a+1.)/2./M_PI/facta/factb); // CONVERSION en MJy Puis les facteurs de la note de G. Lagache (entres dans les classes derivees) - + } } cout << " read "<< vbuf.size() << " values "<0 ) @@ -701,7 +703,7 @@ void HighEll::readAlltSZxCIBTemplates(int lmin, int lmax) { vector vbuf; int nmax_read = 13000; for (int i= 0; (i> a >> b1 >> b2 >> b3 >> b4 >> b5 >> b6>>b7>>b8>>b9>> b10>> b11 >> b12 >> b13 >> b14 >> b15 >> b16>>b17>>b18>>b19>>b20>>b21 ; double spe = 0.; @@ -717,14 +719,14 @@ void HighEll::readAlltSZxCIBTemplates(int lmin, int lmax) { // fact/fact_ref -> extrapole de hfi -> act/spt // puis CONVERSION en MJy avec facteurs de la note de G. Lagache (entres dans les classes derivees) - + } } cout << " read "<< vbuf.size() << " values "<=lmin && ( (int) a <=lmax ) ){ vbuf.push_back(b); //i++ ; @@ -764,7 +766,7 @@ void HighEll::readtSZxCIBTemplate(int lmin, int lmax) { if ( _select_freq[k] !=0) continue; double facta = _cibcc_to_p217[k]; - + for (int j=k; j<_nfreqs ; j++) { if (_select_freq.size()>0 ) @@ -802,7 +804,7 @@ void HighEll::readCIBTemplates(int lmin, int lmax) { int nreadmax = 13000; for (int i=0; i<=nreadmax&& !myfile.eof() ; i++){ myfile >> a >> b ; - //cout << i<< " read " <=lmin && ( (int) a <=lmax ) ){ vbuf.push_back(b); //i++ ; @@ -811,7 +813,7 @@ void HighEll::readCIBTemplates(int lmin, int lmax) { myfile.close(); cout << "readCIBTemplates file closed "< wbuf(vbuf); - + std::transform(wbuf.begin(), wbuf.end(), wbuf.begin(), std::bind2nd(std::multiplies(),sqrt(facta*factb))); - - // std::bind2nd(std::multiplies(), 1./facta/factb )); + + // std::bind2nd(std::multiplies(), 1./facta/factb )); _templ_CIB.push_back( wbuf ); } } cout << " after read cib template "<< _templ_CIB.size()< > HighEll::getCIB(){// std::vector& nuisance vector < vector > temp ( _templ_CIB ); // constructeur par copie ? int ityp=0; - for (int k=0 ; k< _nfreqs ; k++) + for (int k=0 ; k< _nfreqs ; k++) for (int l=k ; l< _nfreqs ; l++) { if (_select_freq.size()>0 ) - if ( ( _select_freq[k] !=0) || ( _select_freq[l] !=0) ) continue ; // ## TOCHECK + if ( ( _select_freq[k] !=0) || ( _select_freq[l] !=0) ) continue ; // ## TOCHECK ityp = getType( _freq[k ], _freq[l ]); std::transform(temp[ityp].begin(), temp[ityp].end(), temp[ityp].begin(), std::bind2nd(std::multiplies(), _A_CIB )); - + } return(temp); @@ -862,21 +864,21 @@ return(temp); vector < vector > HighEll::gettSZxCIB(){// std::vector& nuisance ){ if (_verbose >2) cout << " gettSZxCIB for "<<_name<< " ampli "<<_AtSZxCIB < > temp ( _templ_tSZxCIB ); // constructeur par copie ? - + double fact = _AtSZxCIB ; if (_plik_mode) fact = sqrt(_AtSZ*_A_CIB ) * _AtSZxCIB ; int ityp=0; - for (int k=0 ; k< _nfreqs ; k++) + for (int k=0 ; k< _nfreqs ; k++) for (int l=k ; l< _nfreqs ; l++) { if (_select_freq.size()>0 ) - if ( ( _select_freq[k] !=0) || ( _select_freq[l] !=0) ) continue ; // ## TOCHECK + if ( ( _select_freq[k] !=0) || ( _select_freq[l] !=0) ) continue ; // ## TOCHECK ityp = getType( _freq[k ], _freq[l ]); std::transform(temp[ityp].begin(), temp[ityp].end(), temp[ityp].begin(), std::bind2nd(std::multiplies(), fact )); - - + + } @@ -893,7 +895,7 @@ return(temp); void HighEll::readtSZTemplates(int lmin, int lmax) { cout << " reading tSZ template in "<< _tSZ_file <<" lmin "<< lmin<<" lmax "< vbuf; - // reference = 143 GHz + // reference = 143 GHz double nuref = 143.; double ref_szt = sz_func(nuref); @@ -906,7 +908,7 @@ void HighEll::readtSZTemplates(int lmin, int lmax) { int nreadmax = 13000; for (int i=0; i<=nreadmax&& !myfile.eof() ; i++){ myfile >> a >> b ; - //cout << i<< " read " <=lmin && ( (int) a <=lmax ) ){ vbuf.push_back(b); //i++ ; @@ -918,7 +920,7 @@ void HighEll::readtSZTemplates(int lmin, int lmax) { if (_plik_mode) renorm = vbuf[3000-lmin]; for (int k=0; k<_nfreqs ; k++) - { + { for (int j=k; j<_nfreqs ; j++) { if (_verbose >2) cout << " for : "<(), sz_func( _frtsz[k])* sz_func( _frtsz[j])/(ref_szt*ref_szt*renorm ))); if (_verbose >2) cout << " for : "<> a >> b ; - //cout << i<< " read " <=lmin && ( (int) a <=lmax ) ){ vbuf.push_back(b); // i++; @@ -962,11 +964,11 @@ void HighEll::readkSZTemplates(int lmin, int lmax) { for (int k=0; k<_nfreqs ; k++) { - + for (int j=k; j<_nfreqs ; j++) { if (_verbose >2) cout << " for : "< wbuf(vbuf); + vector wbuf(vbuf); if (_plik_mode) std::transform(wbuf.begin(), wbuf.end(), wbuf.begin(), std::bind2nd(std::multiplies(),1./renorm)); @@ -976,14 +978,14 @@ void HighEll::readkSZTemplates(int lmin, int lmax) { } } cout << " kSZ templates read ! " < > HighEll::gettSZ(){// std::vector& nuisance ){ vector < vector > temp ( _templ_tSZ ); // constructeur par copie ? int ityp=0; - for (int k=0 ; k< _nfreqs ; k++) + for (int k=0 ; k< _nfreqs ; k++) for (int l=k ; l< _nfreqs ; l++) { ityp = getType( _freq[k], _freq[l]); std::transform(temp[ityp].begin(), temp[ityp].end(), temp[ityp].begin(), @@ -1000,7 +1002,7 @@ vector < vector > HighEll::getkSZ(){// std::vector& nuisance ) vector < vector > temp ( _templ_kSZ ); // constructeur par copie ? int ityp=0; - for (int k=0 ; k< _nfreqs ; k++) + for (int k=0 ; k< _nfreqs ; k++) for (int l=k ; l< _nfreqs ; l++) { ityp = getType( _freq[k], _freq[l]); std::transform(temp[ityp].begin(), temp[ityp].end(), temp[ityp].begin(), diff --git a/src/camel/CMB/Hillipop.cc b/src/camel/CMB/Hillipop.cc index 53d5655f27c9481f43addd500927e51ad8075f13..6daab265c095c11b0a44fab9a491eb489a189074 100644 --- a/src/camel/CMB/Hillipop.cc +++ b/src/camel/CMB/Hillipop.cc @@ -37,7 +37,7 @@ void Hillipop::Init(const string fileName) cout << "Using: " << fileName << endl; paramfile parameters(fileName); - + string multipolesFile = Parser::CheckPath(parameters.find("MultipolesRange")); string SZFile = Parser::CheckPath(parameters.find("SZ")); string kSZFile = Parser::CheckPath(parameters.find("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& ClCMB,const vector& nui of << "#Chi2 = " << _chi2 << endl; of << "#Multipole \t"; 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]; @@ -876,7 +876,7 @@ void Hillipop::WriteOutput(const vector& ClCMB,const vector& 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; @@ -888,7 +888,7 @@ void Hillipop::WriteOutput(const vector& ClCMB,const vector& 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; } } @@ -925,8 +925,8 @@ void Hillipop::computeResiduals( vector ClCMB, vector& nuisance, double tmpBeamEigenmodes = 1.; double foregrounds = 0; unsigned int modeOffset = 0; - double ApsRadio, ApsDusty, Asz, Acib, Aksz, Aszxcib; - double AdustTT, AdustTP, AdustPP; + double ApsRadio = 0.0, ApsDusty = 0.0, 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); @@ -1015,7 +1015,7 @@ void Hillipop::computeResiduals( vector ClCMB, vector& nuisance, 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]]; } @@ -1039,7 +1039,7 @@ void Hillipop::computeResiduals( vector ClCMB, vector& 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]]; diff --git a/src/camel/CMB/SPT_pol.cc b/src/camel/CMB/SPT_pol.cc index ce7ba48c7c627875c5a87233c5d4db2c5c502732..2fcf82e61885297aabe8e85bebc6cf56a8f54f3b 100644 --- a/src/camel/CMB/SPT_pol.cc +++ b/src/camel/CMB/SPT_pol.cc @@ -1,6 +1,6 @@ #include "CMB/SPT_pol.hh" #include "Parser.hh" -#include "paramfile.h" +#include "paramfile.h" #include #include #include "float.h" @@ -10,7 +10,9 @@ SPTpol::~SPTpol(){ cout << "Chi2 from SPTpol: " << SPTpolLnLike << endl; - delete _cov_data, _cov_sims, _windows; + delete _cov_data; + delete _cov_sims; + delete _windows; } @@ -18,16 +20,16 @@ SPTpol::SPTpol(const std::string parfile) { paramfile parameters(parfile); _name="SPT_pol"; - + cout << " Likelihood: " << _name << endl; cout << " param file "<< parfile<( "sptpol_bp_file", "")); string cov_data_file = Parser::CheckPath(parameters.find( "sptpol_cov_data_file", "")); string cov_sims_file = Parser::CheckPath(parameters.find( "sptpol_cov_sims_file", "")); @@ -39,11 +41,11 @@ SPTpol::SPTpol(const std::string parfile) spt_windows_lmax = parameters.find( "spt_windows_lmax", 6998); doEE = parameters.find( "EE", true); - if( spt_windows_lmin < 2 || spt_windows_lmin >= spt_windows_lmax) + if( spt_windows_lmin < 2 || spt_windows_lmin >= spt_windows_lmax) throw Message_error( "Invalid lranges for sptpol"); InitSPTpolData( bp_file, cov_data_file, cov_sims_file, win_folder); - + //def nuisance _n.push_back( "SPTpol_Tcal"); _n.push_back( "SPTpol_Pcal"); @@ -66,7 +68,7 @@ void SPTpol::InitSPTpolData( string bp_file, string cov_data_file, string cov_si double ell, cell, val, winval; ifstream fp; - + cout << "nband=" << nband << endl; cout << "nall=" << nall << endl; cout << "bands_per_freq=" << bands_per_freq << endl; @@ -126,7 +128,7 @@ void SPTpol::InitSPTpolData( string bp_file, string cov_data_file, string cov_si } fp.close(); } - + } @@ -138,7 +140,7 @@ double SPTpol::computeLikelihood(const vector& l, vector& clbb, vector& nuisance) { - + // Reshaping the theoretical CMB power spectra Cl (TT,TE,EE,BB) HepMatrix raw_spectra(4,spt_windows_lmax+2,0); int i=0; @@ -150,13 +152,13 @@ double SPTpol::computeLikelihood(const vector& l, i++; } - //The derivative at point n is roughly (raw_spectra[n+1]*l^3 - raw_spectra[n-1]*l^3)/2. + //The derivative at point n is roughly (raw_spectra[n+1]*l^3 - raw_spectra[n-1]*l^3)/2. //Then we devide by l^2 for the proper scaling for the kappa parameter as described in Manzotti, et al. 2014, equation (32). HepMatrix cl_derivative(4,spt_windows_lmax+1,0); for( int j=0; j<3; j++) for(unsigned int l=spt_windows_lmin; l<=spt_windows_lmax; l++) cl_derivative[j][l] = (raw_spectra[j][l+1]*double(l+1)*double(l+1)*double(l+1) - raw_spectra[j][l-1]*double(l-1)*double(l-1)*double(l-1))/2./(double)l/(double)l; - + //DataParams for SPTpol likelihood are: [MapTcal, MapPcal, Czero_EE, CzeroTE, TTepsilon, kappa] const double MapTcal = nuisance[getIndex("SPTpol_Tcal" )]; const double MapPcal = nuisance[getIndex("SPTpol_Pcal" )]; @@ -166,12 +168,12 @@ double SPTpol::computeLikelihood(const vector& l, const double kappa = nuisance[getIndex("SPTpol_kappa" )]; i=0; - double epsilon, calib; + double epsilon = 0.0, calib = 0.0; HepVector deltacb(nall,0); for( int j=0; j& l, //EE (poisson - kappa*deriv) if( k==2) for( int l=spt_windows_lmin; l<=spt_windows_lmax; l++) cl_fgs[l] = raw_spectra[2][l] + ApsEE/3000. - kappa*cl_derivative[2][l]; - + //convert Cl to Dl HepVector dl_th(spt_windows_lmax+1,0); for( int l=spt_windows_lmin; l<=spt_windows_lmax; l++) dl_th[l] = cl_fgs[l] * (double)l*(double)(l+1)/2./M_PI; - for( int b=0; b<=spt_windows_lmax; b++) + for( int b=0; b<=spt_windows_lmax; b++) if( isnan(cl_fgs[b])) { if( i==1) cout << "cl_fgs[TE] "<< b << " " << cl_fgs[b] << " " << raw_spectra[1][b] << " " << cl_derivative[1][b] << endl; if( i==2) cout << "cl_fgs[EE] "<< b << " " << cl_fgs[b] << " " << raw_spectra[2][b] << " " << cl_derivative[2][b] << endl; } - + //Bin dl_th HepVector clbin(nbin,0); // for( int b=0; bsub((i-1)*nbin+1,i*nbin,1,spt_windows_lmax+1)*dl_th; - //If i==1, spec[0] = SPTpol_TE. Scale by Tcal and Pcal and add the TT leakage (epsilon*TT) to the theory bandpowers. + //If i==1, spec[0] = SPTpol_TE. Scale by Tcal and Pcal and add the TT leakage (epsilon*TT) to the theory bandpowers. if( i==1) { epsilon = TTepsilon*MapTcal*MapTcal; calib = MapTcal*MapTcal*MapPcal; @@ -209,24 +211,24 @@ double SPTpol::computeLikelihood(const vector& l, epsilon = TTepsilon*TTepsilon*MapTcal*MapTcal; calib = MapTcal*MapTcal*MapPcal*MapPcal; } - + //add TT leakage clbin += epsilon * _spec[2]; //Compute (model - c*data) deltacb.sub((i-1)*nbin+1, clbin - calib*_spec[i-1]); - + } } - + //We have to scale cov_data by the appropriate calibration parameters. vector cal_vector(nall); - + //If this is the TE portion, only hit with Tcal^2 Pcal. for( int i=0; i<=nbin; i++) cal_vector[i] = MapTcal*MapTcal * MapPcal; //If this is the EE portion, hit with Tcal^2 Pcal^2. for( int i=nbin+1; i& l, cov = cov.sub(1,nbin); } SPTpolLnLike = GaussianLogLike(cov,deltacb); - + _lnlkl = SPTpolLnLike/2.; //return -ln(L) = chi2/2 @@ -254,7 +256,7 @@ double SPTpol::computeLikelihood(const vector& l, // -// Returns -2Log Likelihood for Gaussian: (d^T Cov^{-1} d + log|Cov|) +// Returns -2Log Likelihood for Gaussian: (d^T Cov^{-1} d + log|Cov|) // double SPTpol::GaussianLogLike( HepSymMatrix cov, HepVector deltacb) { @@ -342,7 +344,7 @@ int SPTpol::CholDec(HepSymMatrix* A) *====================================================================*/ { int n = A->num_col(); - + double sum; //Not positive definite @@ -356,7 +358,7 @@ int SPTpol::CholDec(HepSymMatrix* A) //Not positive definite if( sum < DBL_EPSILON) return( -1); - + (*A)[i][i] = sqrt(sum); for( int j=i+1; j std::string str(const T &x){ std::ostringstream os; os << x; return os.str(); -}; +} //specilization template<> std::string str (const float &x){ std::ostringstream os; @@ -86,7 +86,7 @@ void ClassParams::updateVal(const unsigned& i,const string& newval) {pars[i].sec // Constructors -- //---------------- ClassEngine::ClassEngine(const ClassParams& pars): cl(0),dofree(true),pvecback(0){ - + _lmax=-1; //default cout << "Running CLASS version " << _VERSION_ << endl; @@ -95,7 +95,7 @@ ClassEngine::ClassEngine(const ClassParams& pars): cl(0),dofree(true),pvecback(0 size_t n=pars.size(); // parser_init(&fc,n,"pipo",_errmsg); - + //config for (size_t i=0;i%s\n",errmsg); @@ -358,8 +358,8 @@ int ClassEngine::class_main( int ClassEngine::freeStructs(){ - - + + if (lensing_free(&le) == _FAILURE_) { printf("\n\nError in spectra_free \n=>%s\n",le.error_message); return _FAILURE_; @@ -369,23 +369,23 @@ ClassEngine::freeStructs(){ printf("\n\nError in spectra_free \n=>%s\n",sp.error_message); return _FAILURE_; } - + if (transfer_free(&tr) == _FAILURE_) { printf("\n\nError in transfer_free \n=>%s\n",tr.error_message); return _FAILURE_; - } - + } + if (nonlinear_free(&nl) == _FAILURE_) { printf("\n\nError in nonlinear_free \n=>%s\n",nl.error_message); return _FAILURE_; } - + if (primordial_free(&pm) == _FAILURE_) { printf("\n\nError in primordial_free \n=>%s\n",pm.error_message); return _FAILURE_; } - + if (perturb_free(&pt) == _FAILURE_) { printf("\n\nError in perturb_free \n=>%s\n",pt.error_message); @@ -411,7 +411,7 @@ ClassEngine::getCl(Engine::cltype t,const long &l){ if (!dofree) throw out_of_range("no Cl available because CLASS failed"); if (output_total_cl_at_l(&sp,&le,&op,static_cast(l),cl) == _FAILURE_){ - cerr << ">>>fail getting Cl type=" << (int)t << " @l=" << l <>>fail getting Cl type=" << (int)t << " @l=" << l <& lvec, //input - std::vector& cltt, - std::vector& clte, - std::vector& clee, +void +ClassEngine::getCls(const std::vector& lvec, //input + std::vector& cltt, + std::vector& clte, + std::vector& clee, std::vector& clbb) { cltt.resize(lvec.size()); clte.resize(lvec.size()); clee.resize(lvec.size()); clbb.resize(lvec.size()); - + for (size_t i=0;i& lvec, //input } } - -bool -ClassEngine::getLensing(const std::vector& lvec, //input - std::vector& clpp , - std::vector& cltp , + +bool +ClassEngine::getLensing(const std::vector& lvec, //input + std::vector& clpp , + std::vector& cltp , std::vector& clep ){ - + clpp.resize(lvec.size()); cltp.resize(lvec.size()); clep.resize(lvec.size()); - + for (size_t i=0;i unsigned inline + template unsigned inline add(const string& key,const T& val){ pars.push_back(make_pair(key,str(val))); return pars.size(); @@ -97,7 +98,7 @@ public: ClassEngine(const ClassParams& pars); //with a class .pre file ClassEngine(const ClassParams& pars,const string & precision_file); - + // destructor ~ClassEngine(); @@ -123,27 +124,27 @@ public: //throws std::execption if pb virtual std::string name() const {return "ClassEngine";} - double getCl(Engine::cltype t,const long &l); - void getCls(const std::vector& lVec, //input - std::vector& cltt, - std::vector& clte, - std::vector& clee, + double getCl(Engine::cltype t,const long &l); + void getCls(const std::vector& lVec, //input + std::vector& cltt, + std::vector& clte, + std::vector& clee, std::vector& clbb); - - bool getLensing(const std::vector& lVec, //input - std::vector& clphiphi, - std::vector& cltphi, + + bool getLensing(const std::vector& lVec, //input + std::vector& clphiphi, + std::vector& cltphi, std::vector& clephi); //recombination inline double z_rec() const {return th.z_rec;} - inline double rs_rec() const {return th.rs_rec;} + inline double rs_rec() const {return th.rs_rec;} - //baryon drag + //baryon drag inline double z_drag() const {return th.z_d;} - inline double rs_drag() const {return th.rs_d;} + inline double rs_drag() const {return th.rs_d;} double get_Dv(double z); double get_Da(double z); @@ -164,13 +165,13 @@ public: bool get_PklinNodes(std::vector& knodes,std::vector& Pknodes); bool get_PkNLNodes(std::vector& knodes,std::vector& Pknodes); // the followingdepends on switch (prefer previous methods) - inline bool get_PkNodes(std::vector& knodes,std::vector& Pknodes) - {return (has_PkNL()? get_PkNLNodes(knodes,Pknodes) : get_PklinNodes(knodes,Pknodes));} + inline bool get_PkNodes(std::vector& knodes,std::vector& Pknodes) + {return (has_PkNL()? get_PkNLNodes(knodes,Pknodes) : get_PklinNodes(knodes,Pknodes));} //vectorized access: faster (uses class_extra from camel) // outout vector pks must be already sized (see Engine.hh) - bool get_Pkvec(const std::vector& kin, const std::vector& z,std::vector& pks); - + bool get_Pkvec(const std::vector& kin, const std::vector& z,std::vector& pks); + inline bool get_Pkvec(const std::vector& kvec,const double& z, std::vector& pks){ std::vector zvec(1,z); return this->get_Pkvec(kvec,zvec,pks); @@ -251,7 +252,7 @@ private: protected: - + std::string _output; double* pvecback; // helper to avoid too much allocations @@ -259,4 +260,3 @@ protected: }; #endif - diff --git a/src/camel/Engine.cc b/src/camel/Engine.cc index 18288e116590c24354039c6ee697084c32bbae4c..665cb8eb665b6ac14e09974ed74ef3dc5a8c39d0 100644 --- a/src/camel/Engine.cc +++ b/src/camel/Engine.cc @@ -104,7 +104,8 @@ Engine::get(const std::string& key) planck_assert(key.substr(key.size()-1,1)==")",key+": keyword not finishing with )"); var=key.substr(0,index); istringstream isNum(key.substr(index+1,key.size()-index-1)); - planck_assert(isNum>>z,key+": Engine canot decode z value in keyword"); + isNum >> z; + planck_assert(isNum.fail(),key+": Engine canot decode z value in keyword"); } if (var=="z_drag") return z_drag(); else if (var=="rs_drag") return rs_drag(); diff --git a/src/camel/JLA/JLA_chi2.cc b/src/camel/JLA/JLA_chi2.cc index 592c6f35c653efe4c5449935ebe5a7b073562efa..bedaeaf5cced39e962387c7aa89235b2bd0bb4e6 100644 --- a/src/camel/JLA/JLA_chi2.cc +++ b/src/camel/JLA/JLA_chi2.cc @@ -50,37 +50,37 @@ std::vector JLA_chi2::requires() const{ double JLA_chi2::chi2_eff(const std::vector& par) const { double _jla_nuisance[4]; - // fills nuisance param array + // fills nuisance param array _jla_nuisance[0] = par[_alpha_index]; _jla_nuisance[1] = par[_beta_index]; - _jla_nuisance[2] = par[_M_index]; - _jla_nuisance[3] = par[_DeltaM_index]; + _jla_nuisance[2] = par[_M_index]; + _jla_nuisance[3] = par[_DeltaM_index]; // gets distance moduli from class - if (_verbosity>3) + if (_verbosity>3) cout << name() << " chi2_eff nuisance params : "<< _jla_nuisance[0] <<" "<< _jla_nuisance[1] <<" "<< _jla_nuisance[2] << " "<< _jla_nuisance[3]<size()]; + std::vector mu( _jlalike->size()); for (int i=0; i < _jlalike->size(); i++) { - mu[i] = engine->get_DMod (_Z[i] ); + mu[i] = engine->get_DMod (_Z[i] ); } // JLA likelihood evaluation - double retchi2 = _jlalike->computeLikelihood(mu,_jla_nuisance ); + double retchi2 = _jlalike->computeLikelihood(&mu[0],_jla_nuisance ); if (_verbosity>4) { - // puts distance modulus and standardized residuals in an ascii file for later plotting + // puts distance modulus and standardized residuals in an ascii file for later plotting ofstream myfile; - double residuals[ _jlalike->size() ] ; - int ret = _jlalike->computeResiduals(mu, _jla_nuisance, residuals); + std::vector residuals( _jlalike->size()); + int ret = _jlalike->computeResiduals(&mu[0], _jla_nuisance, &residuals[0]); if (ret !=0) cout<<" PROBLEM WITH COMPUTE RESIDUALS !!" <size(); i++){ myfile<< i << " " << _Z[i]<<" "<< mu[i]<<" "<< residuals[i] <<"\n"; } myfile.close(); - - // + + // } return(retchi2); @@ -89,16 +89,16 @@ double JLA_chi2::chi2_eff(const std::vector& par) const void JLA_chi2::read_par_file(const std::string fileName) { - + // string pathToDustSpectra = parameters.find("Dust"); // _modeStatus[0] = parameters.find("TT"); - + cout << "using:"<< fileName << endl; paramfile parameters(fileName); // _dataset_name = (Parser::CheckPath(parameters.find("DataSet",""))).c_str(); // cout << " dataset "<< _dataset_name << endl; //Parser::CheckPath(parser.params.find("BAO2DFile","")); - _verbosity = 3 ; + _verbosity = 3 ; if (parameters.param_present("verbose")) _verbosity = parameters.find("verbose"); _jlalike = new JLALikelihood (_verbosity); @@ -108,7 +108,7 @@ void JLA_chi2::read_par_file(const std::string fileName) _jlalike->read( (Parser::CheckPath(fileName)).c_str() ); _Z = _jlalike->getZ(); - // + // } diff --git a/src/camel/JLA/jla_likelihood_v3/src/jla.cc b/src/camel/JLA/jla_likelihood_v3/src/jla.cc index bb5598db6fabebe2153f696afe0232a9185f9f9c..2938d3d1267211fa5f34a11919b2b96f13deda27 100644 --- a/src/camel/JLA/jla_likelihood_v3/src/jla.cc +++ b/src/camel/JLA/jla_likelihood_v3/src/jla.cc @@ -147,11 +147,11 @@ void JLALikelihood::read(const char * datasetfile) cout << "Reading config from " << datasetfile << endl; - // TEST + // TEST if (ini_parse(datasetfile, configuration_handler, &config) < 0) { - if (verbosity > 0) + if (verbosity > 0) cerr << "Can't load '" << datasetfile << "'\n"; exit(-1); } @@ -178,7 +178,7 @@ void JLALikelihood::read(const char * datasetfile) cerr << "Can't load '" << config.data_file << "'\n"; exit(-1); } - + while(fid.get() == '#') fid.getline(buffer, 1024); fid.unget(); @@ -189,7 +189,7 @@ void JLALikelihood::read(const char * datasetfile) if (fid) lcpars.push_back(SN); } - + if (verbosity > 1) cout << "Read " << lcpars.size() << " SNe in file " << config.data_file << endl; } @@ -218,10 +218,11 @@ void JLALikelihood::read(const char * datasetfile) double JLALikelihood::computeLikelihood(double * distance_modulii, double * nuisance_parameters) { - double residuals[size()]; + + std::vector residuals(size()); int status; - status = computeResiduals(distance_modulii, nuisance_parameters, residuals); - + status = computeResiduals(distance_modulii, nuisance_parameters, &residuals[0]); + double chi2 = 0; if (status == 0) { @@ -258,7 +259,7 @@ double JLALikelihood::computeLikelihood(double * distance_modulii, * standardized residuals r_i at the end of the execution. * The minization criterion is $\chi^2 = \sum_i r_i^2$. * - * Return: + * Return: * ------- * 0 if the computation is successful, -1 otherwise. */ @@ -275,15 +276,15 @@ int JLALikelihood::computeResiduals(double * distance_modulii, double * nuisance // Covariance matrix computation int nsq = n*n; - double cov[nsq]; - cblas_dcopy(nsq, C00, incx, cov, incy); - cblas_daxpy(nsq, sq(alpha), C11, 1, cov, 1); - cblas_daxpy(nsq, sq(beta), C22, 1, cov, 1); - cblas_daxpy(nsq, 2.0 * alpha, C01, 1, cov, 1); - cblas_daxpy(nsq, -2.0 * beta, C02, 1, cov, 1); - cblas_daxpy(nsq, -2.0 * alpha * beta, C12, 1, cov, 1); - - + std::vector cov(nsq); + cblas_dcopy(nsq, C00, incx, &cov[0], incy); + cblas_daxpy(nsq, sq(alpha), C11, 1, &cov[0], 1); + cblas_daxpy(nsq, sq(beta), C22, 1, &cov[0], 1); + cblas_daxpy(nsq, 2.0 * alpha, C01, 1, &cov[0], 1); + cblas_daxpy(nsq, -2.0 * beta, C02, 1, &cov[0], 1); + cblas_daxpy(nsq, -2.0 * alpha * beta, C12, 1, &cov[0], 1); + + for (int i = 0; i < n; i++) { LCPar sn = lcpars[i]; @@ -298,17 +299,17 @@ int JLALikelihood::computeResiduals(double * distance_modulii, double * nuisance - 2.0 * alpha * beta * sn.cov_s_c; } - + // Whiten the residuals int nhrs = 1, info = 0; - dpotrf("U", &n, cov, &n, &info); - + dpotrf("U", &n, &cov[0], &n, &info); + if (info != 0){ if (verbosity > 0) cerr << "Cholesky Error: " << info << endl; return -1; } - dtrtrs("U", "T", "N", &n, &nhrs, cov, &n, residuals, &n, &info); + dtrtrs("U", "T", "N", &n, &nhrs, &cov[0], &n, residuals, &n, &info); if (info != 0){ if (verbosity > 0) cerr << "Solve Error: " << info << endl; @@ -357,4 +358,3 @@ JLALikelihood::~JLALikelihood() free(C12); } } - diff --git a/src/camel/MCMC/exec/mainMCMC.cc b/src/camel/MCMC/exec/mainMCMC.cc index 4f257e0cc421973ed4ed682ab0d12bd2a23b59a6..d3b85d2b223681c88b8c44f570f27cb0b0bad21b 100644 --- a/src/camel/MCMC/exec/mainMCMC.cc +++ b/src/camel/MCMC/exec/mainMCMC.cc @@ -15,7 +15,7 @@ #include "CLHEP/Matrix/SymMatrix.h" #include "CLHEP/Matrix/Matrix.h" #include -//planck lkh +//planck lkh #include "Chi2Factory.hh" #include "cxxsupport/paramfile.h" @@ -27,20 +27,20 @@ #include //#include -#include +#include #include -using namespace std; +using namespace std; int main(int argc, char *argv[]) { - + planck_assert(argc==3,"Usage: mcmc parfile(in) txtfile(out)"); planck_assert(file_present(argv[1]),"missing par file"); - - Parser parser(argv[1]); + + Parser parser(argv[1]); const char* chain_run_name=argv[2]; const string fileName(chain_run_name); @@ -67,17 +67,17 @@ int main(int argc, char *argv[]) bool use_ext_cor=(proposal_cor.size()!=0); bool use_ext_cov=(proposal_cov.size()!=0); - + //--------------------------------------------------------- - + Chi2Combiner* chi2=Chi2Factory::gimeChi2(parser); chi2->setVerbose(false); planck_lkh* lkh=new planck_lkh(chi2,parser); - const int dim=lkh->Get_par_num(); + const int dim=lkh->Get_par_num(); cout <<"dim=" <SetBunchSize(bunchSize); //directly from Parser - chain->set_par(parser); - + chain->set_par(parser); + if (verbose){ cout << chain->min_par() << endl; cout << chain->max_par() << endl; cout << chain->starting_point() << endl; } - //or by hand using set/get functions - + //or by hand using set/get functions + vector par_err(dim); for(int i=0;ipar_err())[i]; if (verbose) cout << par_err[i] <starting_point())[i]; } - + HepSymMatrix cov(dim,1); HepSymMatrix identity(dim,1); - + //-----define the lkh---------------------------------------------- pdf_gaussian* proposal= new pdf_gaussian(dim,mean_pro,identity); @@ -157,7 +157,7 @@ int main(int argc, char *argv[]) } else { if(use_ext_cor) - { + { cout << "Read proposal cor matrix.." << endl; cout << proposal_cor.c_str() << endl; proposal->set_cov(proposal_cor.c_str()); @@ -166,7 +166,7 @@ int main(int argc, char *argv[]) cov = proposal->Get_cov_matrix(); proposal->calc_cov_matrix( cov, par_err); } - + if (verbose) cout << "cov proposal= " << proposal->Get_cov_matrix() << endl; //scale @@ -175,7 +175,7 @@ int main(int argc, char *argv[]) vector scale_vs_length; proposal->init_generator(seed); - + //if we want to displace the starting point @@ -192,7 +192,7 @@ int main(int argc, char *argv[]) chain->Set_starting_point(chain->starting_point()+start); if (verbose) cout << "start_new=" << chain->starting_point() << endl; } - + chain->Set_burn_in(0); @@ -204,8 +204,8 @@ int main(int argc, char *argv[]) //pdf* p=new prior_beta(chain->_min_par[i],chain->_max_par[i], 1,1); prior_list[i]=p; } - - + + prior* flat_p=new prior(prior_list); pdf* target=new post_planck(lkh,flat_p); @@ -222,26 +222,26 @@ int main(int argc, char *argv[]) //proposal->write_cov("cov_proposal.txt"); - + //acceptance rate simple double rate=chain->acc_rate(0,length); if (verbose) cout << "acc_rate=" << rate << endl ; - + //---autocorr----and acc rate------------- vector ar_vs_length; - int max_lag; - int min_lag; + // int max_lag; + // int min_lag; vector autocorr_vs_length; for(int step=100; step<=length; step+=100) { - max_lag=50; - min_lag=max_lag; + // max_lag=50; + // min_lag=max_lag; if((stepacc_rate(0,step)); } else if(step>=t0 && algo=="ada" ) {ar_vs_length.push_back(chain->acc_rate(step-100,step)); } - else + else {ar_vs_length.push_back(chain->acc_rate(0,step));} } @@ -251,18 +251,17 @@ int main(int argc, char *argv[]) { MCMC_adaptive* ada=dynamic_cast(chain); chain->write_if_useful(proposal->_scale_vs_length,"scale_vs_length.txt"); - - ofstream f("corr.txt"); + + ofstream f("corr.txt"); if(!f) { cout << "Error in file creation!"; } f << ada->mean() << endl //mean << proposal->Get_cov_matrix() << endl; //cov matrix - f.close(); + f.close(); cout<<"output written in:corr.txt"< void printVec(const vector& v,unsigned n,const string & s){ - cout << s << " (size=" << v.size() <<") :" ; - if (v.size()(cout,"\t")); - cout << "..." << endl; - } -; + cout << s << " (size=" << v.size() <<") :" ; + if (v.size()(cout,"\t")); + cout << "..." << endl; +} /*********************/ /* global variables */ /*********************/ @@ -59,7 +58,7 @@ void dumpMin(ostream& o,const FunctionMinimum& min1){ o << "IsAboveMaxEdm " << min1.IsAboveMaxEdm() < firstpar=migrad.Params(); if(rm_cosmo_lim) { for (size_t i=0;i lastpar=migrad.Params(); double chi2=(*theChi2)(lastpar); bool valid=min1.IsValid(); - + cout << "Profile" << " " ; for(int ii=0;iigetEngine()->get(parser.derivedVars()[i]); }catch (Message_error &e){ std::cerr << "could not get: " << e.what() << std::endl; - + } of << val <<"\t"; } @@ -177,12 +176,12 @@ int main(int argc,char** argv){ sprintf(DeChi2,"%5.7f",chi2); of << DeChi2 << " " << " " << valid << " " << min1.Edm() ; of << " " << min1.HasReachedCallLimit() - << " " << min1.HasValidParameters() - << " " << min1.HasCovariance() - << " " << min1.HasValidCovariance() - << " " << min1.HasAccurateCovar() - << " " << min1.HasPosDefCovar() - << " " << min1.HesseFailed() + << " " << min1.HasValidParameters() + << " " << min1.HasCovariance() + << " " << min1.HasValidCovariance() + << " " << min1.HasAccurateCovar() + << " " << min1.HasPosDefCovar() + << " " << min1.HesseFailed() << " " << min1.IsAboveMaxEdm() << endl; } @@ -190,7 +189,6 @@ int main(int argc,char** argv){ if (argc ==4){ ofstream of(argv[3]); const MnUserCovariance& cov=min1.UserCovariance(); - int pr = of.precision(6); unsigned int n = cov.Nrow(); for (unsigned int i = 0; i < n; i++) { for (unsigned int j = 0; j < n; j++) { diff --git a/src/camel/exec/writeSpectraPk.cc b/src/camel/exec/writeSpectraPk.cc index 320c1b2516a575b00b0c916f20d5dd9116a10528..8921ee280cc9a6fda6889536b668c8e6764f8a3b 100644 --- a/src/camel/exec/writeSpectraPk.cc +++ b/src/camel/exec/writeSpectraPk.cc @@ -36,7 +36,7 @@ int main(int argc,char** argv){ planck_assert(file_present(argv[1]),"mssing par file"); - float z0 = 0.; + float z0 = 0.; if (argc==5) z0=atof(argv[4]); @@ -59,7 +59,7 @@ int main(int argc,char** argv){ classparms.add("z_pk",z0); //polar +lens+clphi //classparms.add("z_max_pk",z0); //polar +lens+clphi - + //increase kmax if higher value requested const double kmax = parser.params.find("k_end",1); try{ @@ -78,7 +78,7 @@ int main(int argc,char** argv){ if (parser.params.param_present("precisionFile")){ pre=Parser::CheckPath(parser.params.find("precisionFile","")); } - + cout << "CLASS \t--> precision parameters taken from file=" << pre << endl; e=new MnClassEngine(parser.vars(), @@ -136,7 +136,7 @@ int main(int argc,char** argv){ cols.push_back(fitscolumn("TE", "K",1,TDOUBLE)); fout.insert_bintab (cols); - + //lvec>2 vector lvec(lmax-1,1); lvec[0]=2; @@ -149,7 +149,7 @@ int main(int argc,char** argv){ arr cl(lmax+1); cl[0]=0; cl[1]=0; - + //en K //TT for (size_t i=2;i<=lmax;i++) cl[i]=1e-12*cltt[i-2]; @@ -189,7 +189,7 @@ int main(int argc,char** argv){ } - // pk + // pk cout << "Computing P(k) at z=" << z0 << endl; cout << "NL spectrum will" << (!do_nl? " NOT ":" ") << "be computed" <("n_pks",1000); const double kmin = parser.params.find("k_beg",8.e-4); - logspace(kmin,kmax,npks,kval) ; + logspace(kmin,kmax,npks,kval) ; } //test fast Pk @@ -238,11 +238,11 @@ int main(int argc,char** argv){ for (size_t i=0;i<100;i++){ for (int k =0; kget_Pklin(kval[k], 0); + e->get_Pklin(kval[k], 0); } } cout << "std=" << timer.partial() < pklin; vector pknl ; for (int k =0; k pkfast(kval.size()); try{ e->get_Pkvec(kval,z0,pkfast); - } + } catch (Message_error &e){}; int a(0); @@ -273,6 +273,6 @@ int main(int argc,char** argv){ fout.close(); } delete e; - + } diff --git a/src/cxxsupport/cxxutils.cc b/src/cxxsupport/cxxutils.cc index e2618a559f9907bee56bbcba4439f36103878e22..e5ee4e9975de44fb63f8a22473e1c5d1c5c95b0f 100644 --- a/src/cxxsupport/cxxutils.cc +++ b/src/cxxsupport/cxxutils.cc @@ -31,7 +31,7 @@ using namespace std; bool file_present (const string &filename) { ifstream dummy(filename.c_str()); - return dummy; + return dummy.good(); } void assert_present (const string &filename) @@ -214,15 +214,6 @@ void IO_backend_status() << iohandle_current::Type() << " I/O." << endl; } -void CVS_status() - { -#ifndef CVSTAG - cout << "CVS checkout status: unknown" << endl; -#else - cout << "CVS checkout tag: " CVSTAG " (" CVSSTATUS ")" << endl; -#endif - } - } //unnamed namespace void announce (const string &name) diff --git a/src/cxxsupport/datatypes.h b/src/cxxsupport/datatypes.h index 7a5f5c4216866e275bfba8c735c5abc1af5d3a82..b43012ca500cbff07adfcbd3a541f958eebd3923 100644 --- a/src/cxxsupport/datatypes.h +++ b/src/cxxsupport/datatypes.h @@ -18,54 +18,73 @@ // should not be used outside this file. template struct sizeChooserHelper - { typedef void TYPE; }; +{ typedef void TYPE; }; template struct sizeChooserHelper - { typedef T TYPE; }; +{ typedef T TYPE; }; template struct sizeChooserHelper2 - { typedef T1 TYPE; }; +{ typedef T1 TYPE; }; template struct sizeChooserHelper2 - { typedef T2 TYPE; }; +{ typedef T2 TYPE; }; template struct sizeChooserHelper2 - { typedef T3 TYPE; }; +{ typedef T3 TYPE; }; template <> struct sizeChooserHelper2 - { }; +{ }; template - struct sizeChooser - { +struct sizeChooser +{ typedef typename sizeChooserHelper2 - ::TYPE, - typename sizeChooserHelper::TYPE, - typename sizeChooserHelper::TYPE >::TYPE TYPE; - }; + ::TYPE, + typename sizeChooserHelper::TYPE, + typename sizeChooserHelper::TYPE >::TYPE TYPE; +}; +#if (__cplusplus>=201103L) + +#include + +typedef int8_t int8; +typedef uint8_t uint8; + +typedef int16_t int16; +typedef uint16_t uint16; + +typedef int32_t int32; +typedef uint32_t uint32; + +typedef int64_t int64; +typedef uint64_t uint64; + +#else typedef signed char int8; typedef unsigned char uint8; typedef sizeChooser<2, short, int>::TYPE - int16; +int16; typedef sizeChooser<2, unsigned short, unsigned int>::TYPE - uint16; +uint16; typedef sizeChooser<4, int, long, short>::TYPE - int32; +int32; typedef sizeChooser<4, unsigned int, unsigned long, unsigned short>::TYPE - uint32; +uint32; typedef sizeChooser<8, long, long long>::TYPE - int64; +int64; typedef sizeChooser<8, unsigned long, unsigned long long>::TYPE - uint64; +uint64; + +#endif typedef sizeChooser<4, float, double>::TYPE - float32; +float32; typedef sizeChooser<8, double, long double>::TYPE - float64; +float64; //! unsigned integer type which should be used for array sizes typedef std::size_t tsize; @@ -89,32 +108,32 @@ enum { PLANCK_INT8 = 0, template struct typehelper {}; template<> struct typehelper - { enum { id=PLANCK_INT8 }; }; +{ enum { id=PLANCK_INT8 }; }; template<> struct typehelper - { enum { id=PLANCK_UINT8 }; }; +{ enum { id=PLANCK_UINT8 }; }; template<> struct typehelper - { enum { id=PLANCK_INT16 }; }; +{ enum { id=PLANCK_INT16 }; }; template<> struct typehelper - { enum { id=PLANCK_UINT16 }; }; +{ enum { id=PLANCK_UINT16 }; }; template<> struct typehelper - { enum { id=PLANCK_INT32 }; }; +{ enum { id=PLANCK_INT32 }; }; template<> struct typehelper - { enum { id=PLANCK_UINT32 }; }; +{ enum { id=PLANCK_UINT32 }; }; template<> struct typehelper - { enum { id=PLANCK_INT64 }; }; +{ enum { id=PLANCK_INT64 }; }; template<> struct typehelper - { enum { id=PLANCK_UINT64 }; }; +{ enum { id=PLANCK_UINT64 }; }; template<> struct typehelper - { enum { id=PLANCK_FLOAT32 }; }; +{ enum { id=PLANCK_FLOAT32 }; }; template<> struct typehelper - { enum { id=PLANCK_FLOAT64 }; }; +{ enum { id=PLANCK_FLOAT64 }; }; template<> struct typehelper - { enum { id=PLANCK_BOOL }; }; +{ enum { id=PLANCK_BOOL }; }; template<> struct typehelper - { enum { id=PLANCK_STRING }; }; +{ enum { id=PLANCK_STRING }; }; inline int type2size (int type) - { +{ switch (type) { case PLANCK_INT8 : return 1; @@ -132,10 +151,10 @@ inline int type2size (int type) default: planck_fail ("type2size: unsupported data type"); } - } +} inline int string2type(const std::string &type) - { +{ if (type=="FLOAT64") return PLANCK_FLOAT64; if (type=="FLOAT32") return PLANCK_FLOAT32; if (type=="INT8") return PLANCK_INT8; @@ -149,10 +168,10 @@ inline int string2type(const std::string &type) if (type=="BOOL") return PLANCK_BOOL; if (type=="STRING") return PLANCK_STRING; planck_fail ("string2type: unknown data type '"+type+"'"); - } +} inline const char *type2string (int type) - { +{ switch (type) { case PLANCK_INT8 : return "INT8"; @@ -170,37 +189,37 @@ inline const char *type2string (int type) default: planck_fail ("type2string: unsupported data type"); } - } +} template inline const char *type2typename () - { return "unknown type"; } +{ return "unknown type"; } template<> inline const char *type2typename () - { return "signed char"; } +{ return "signed char"; } template<> inline const char *type2typename () - { return "unsigned char"; } +{ return "unsigned char"; } template<> inline const char *type2typename () - { return "short"; } +{ return "short"; } template<> inline const char *type2typename () - { return "unsigned short"; } +{ return "unsigned short"; } template<> inline const char *type2typename () - { return "int"; } +{ return "int"; } template<> inline const char *type2typename () - { return "unsigned int"; } +{ return "unsigned int"; } template<> inline const char *type2typename () - { return "long"; } +{ return "long"; } template<> inline const char *type2typename () - { return "unsigned long"; } +{ return "unsigned long"; } template<> inline const char *type2typename () - { return "long long"; } +{ return "long long"; } template<> inline const char *type2typename () - { return "unsigned long long"; } +{ return "unsigned long long"; } template<> inline const char *type2typename () - { return "float"; } +{ return "float"; } template<> inline const char *type2typename () - { return "double"; } +{ return "double"; } template<> inline const char *type2typename () - { return "bool"; } +{ return "bool"; } template<> inline const char *type2typename () - { return "std::string"; } +{ return "std::string"; } #endif