Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Open sidebar
Xavier Garrido
CAMEL
Commits
28ed2caa
Commit
28ed2caa
authored
Jul 19, 2018
by
Xavier Garrido
Committed by
Matthieu Tristram
Jul 19, 2018
Browse files
Fix warnings
parent
2cef189c
Changes
15
Expand all
Hide whitespace changes
Inline
Side-by-side
Showing
15 changed files
with
547 additions
and
538 deletions
+547
-538
src/AbsRand/OofNoise.cc
src/AbsRand/OofNoise.cc
+41
-42
src/camel/CMB/HiLLiPOP.cc
src/camel/CMB/HiLLiPOP.cc
+46
-47
src/camel/CMB/HighEll.cc
src/camel/CMB/HighEll.cc
+138
-136
src/camel/CMB/Hillipop.cc
src/camel/CMB/Hillipop.cc
+36
-36
src/camel/CMB/SPT_pol.cc
src/camel/CMB/SPT_pol.cc
+32
-30
src/camel/Class/ClassEngine.cc
src/camel/Class/ClassEngine.cc
+48
-49
src/camel/Class/ClassEngine.hh
src/camel/Class/ClassEngine.hh
+21
-21
src/camel/Engine.cc
src/camel/Engine.cc
+2
-1
src/camel/JLA/JLA_chi2.cc
src/camel/JLA/JLA_chi2.cc
+16
-16
src/camel/JLA/jla_likelihood_v3/src/jla.cc
src/camel/JLA/jla_likelihood_v3/src/jla.cc
+22
-22
src/camel/MCMC/exec/mainMCMC.cc
src/camel/MCMC/exec/mainMCMC.cc
+36
-37
src/camel/exec/Minimize.cc
src/camel/exec/Minimize.cc
+25
-27
src/camel/exec/writeSpectraPk.cc
src/camel/exec/writeSpectraPk.cc
+11
-11
src/cxxsupport/cxxutils.cc
src/cxxsupport/cxxutils.cc
+1
-10
src/cxxsupport/datatypes.h
src/cxxsupport/datatypes.h
+72
-53
No files found.
src/AbsRand/OofNoise.cc
View file @
28ed2caa
...
...
@@ -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 @
28ed2caa
...
...
@@ -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 @
28ed2caa
This diff is collapsed.
Click to expand it.
src/camel/CMB/Hillipop.cc
View file @
28ed2caa
...
...
@@ -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
());