Commit 07994ce5 authored by Plaszczynski Stephane's avatar Plaszczynski Stephane

add Olivier test

parent 3cdd6261
#include "Parser.hh"
#include "Class/MnClassEngine.hh"
#ifdef PICO
#include "pico/MnPicoEngine.hh"
#endif
#include"arr.h"
#include"fitshandle.h"
#include <iostream>
#include<fstream>
#include<cassert>
#include<numeric>
#include<string>
using namespace std;
int main(int argc,char** argv){
planck_assert(argc==4,"Usage: spectra.exe parfile(in) lmax(in) file.fits(out)");
planck_assert(file_present(argv[1]),"mssing par file");
//decodage arguments
Parser parser(argv[1]);
int lmax=atoi(argv[2]);
//choose engine
string engine= parser.params.find<string>("engine","CLASS");
Engine *e;
if (engine=="pico") {
throw Message_error("pico cannot be used for Pk !");
#ifndef PICO
throw Message_error("pico support not compiled: see reqiurement");
#else
e=new MnPicoEngine(parser.vars(),lmax,PICODATA,true);
#endif
}
else{
//CLASS
ClassParams classparms(parser.classparms);
classparms.add("l_max_scalars",lmax);
classparms.add("output","tCl,pCl,lCl,mPk"); //polar +lens+clphi
classparms.add("z_pk","0,0.32,0.57"); //polar +lens+clphi
int do_nl =0;
try{
size_t test = classparms.findIndex("non linear");
cout << " >>> CLASS nl par found !! "<< test<<endl;
do_nl = 1 ;
}catch (Message_error& msg){
cout << " >>> CLASS nl par NOT found !! " <<endl;
}
for (size_t i=0;i<classparms.size();i++)
cout << "CLASS \t--> "<< classparms.key(i) << "\t" << classparms.value(i) <<endl;
//precision file
string pre="";
if (parser.params.param_present("precisionFile")){
pre=Parser::getParFile(parser.params.find<string>("precisionFile",""));
}
cout << "CLASS \t--> precision parameters taken from file=" << pre << endl;
e=new MnClassEngine(parser.vars(),
classparms,
pre);
//derived
//expansion
double H0=e->get_H0();
double theta=e->theta();
cout << " expansion: H0=" << H0 << " 100*theta_s=" << 100*theta << endl << endl;;
cout << " Omega_m=" << e->Omega_m() << endl ;
//YHE
cout << " YHe=" << e->YHe() << endl << endl;;
//reio
double tau=e->get("tau_reio");
double zreio=e->get("z_reio");
cout << " reionization : tau_reio=" << tau << " z_reio=" << zreio << endl << endl;;
//BAO
double rfid=149.28 ; // from BOSS Anderson 2014
double rd=e->rs_drag();
double r=rd/rfid;
cout << " Dv*(rfid/rd) @z=0.32: " << e->get_Dv(0.32)/r << " @z=0.57: " << e->get_Dv(0.57)/r << endl;
cout << " Da*(rfid/rd): @z=0.57: " << e->get_Da(0.57)/r << " H=" << e->get_H(0.57) << endl << endl;;
//sigma8
double s8=e->get_sigma8(0.);
cout << " sigma8=" << s8 << endl;
cout << " s8*(Omega_m/0.27)^0.3=" << s8*pow(e->Omega_m()/0.27,0.3) << endl;
cout << " f sigma8(@z=0.57)=" << e->get_f(0.57)*e->get_sigma8(0.57) << endl;
//ecriture output fits
const std::string fileName(argv[3]);
cout << "writing FITS file for Cells =" << fileName << endl;
if (file_present(fileName)) {
std::cout << "the file " + fileName + " already exists: removing"<<endl;
remove_file(fileName);
}
fitshandle fout(fileName,fitshandle::CREATE, READWRITE);
std::vector<fitscolumn> cols;
//definition des colonnes
cols.push_back(fitscolumn("TT", "K",1,TDOUBLE));
cols.push_back(fitscolumn("EE", "K",1,TDOUBLE));
cols.push_back(fitscolumn("BB", "K",1,TDOUBLE));
cols.push_back(fitscolumn("TE", "K",1,TDOUBLE));
fout.insert_bintab (cols);
//lvec>2
vector<unsigned> lvec(lmax-1,1);
lvec[0]=2;
partial_sum(lvec.begin(),lvec.end(),lvec.begin());
vector<double> cltt,clte,clee,clbb;
e->getCls(lvec,cltt,clte,clee,clbb);
{
arr<double> 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];
fout.write_column(1,cl);
//EE
for (size_t i=2;i<=lmax;i++) cl[i]=1e-12*clee[i-2];
fout.write_column(2,cl);
//BB
for (size_t i=2;i<=lmax;i++) cl[i]=1e-12*clbb[i-2];
fout.write_column(3,cl);
//TE
for (size_t i=2;i<=lmax;i++) cl[i]=1e-12*clte[i-2];
fout.write_column(4,cl);
}
//add lensing if available
vector<double> clpp,cltp,clep;
if (e->getLensing(lvec,clpp,cltp,clep)) {
std::vector<fitscolumn> cols2;
//definition des colonnes
cols2.push_back(fitscolumn("PP", "K",1,TDOUBLE));
cols2.push_back(fitscolumn("TP", "K",1,TDOUBLE));
cols2.push_back(fitscolumn("EP", "K",1,TDOUBLE));
fout.insert_bintab(cols2);
{
arr<double> cl(lmax+1);
cl[0]=0;
cl[1]=0;
for (size_t i=2;i<=lmax;i++) cl[i]=clpp[i-2];
fout.write_column(1,cl);
for (size_t i=2;i<=lmax;i++) cl[i]=cltp[i-2];
fout.write_column(2,cl);
for (size_t i=2;i<=lmax;i++) cl[i]=clep[i-2];
fout.write_column(3,cl);
}
}
// pk at z=0
//ecriture output fits
std::vector<fitscolumn> pkcols;
pkcols.push_back(fitscolumn("k", "h/Mpc",1,TDOUBLE));
pkcols.push_back(fitscolumn("Pklin", "(h/Mpc)^-3" ,1,TDOUBLE));
pkcols.push_back(fitscolumn("Pknl", "(h/Mpc)^-3",1,TDOUBLE));
fout.insert_bintab(pkcols);
int npks = 1000 ;
double kmax = 0.25;
vector<double> pklin;
vector<double> pknl ;
vector<double> kval ;
double z = 0.;
for (int k =0; k<npks ; k++){
kval.push_back(kmax *((double (k)+1.)/ npks )) ;
double myk = kval[k];
//double toto = e->get_test(myk, z);
pklin.push_back( e->get_Pk(myk, z) );
if ( do_nl == 1) pknl.push_back( e->get_PkNL(myk, z) );
}
arr<double> apk(npks);
for (size_t i=0; i<npks ; i++) apk[i]=kval[i];
fout.write_column(1,apk);
for (size_t i=0; i<npks ; i++) apk[i]=pklin[i];
fout.write_column(2,apk);
if ( do_nl == 1){
for (size_t i=0; i<npks ; i++) apk[i]=pknl[i];
fout.write_column(3,apk);
}
fout.close();
}
delete e;
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment