Commit 38e61f3b authored by Matthieu Tristram's avatar Matthieu Tristram
Browse files

Add Billipop for Bmodes recombination scales

parent b4e89a2e
#include "Billipop.hh"
#include "Parser.hh"
#include "paramfile.h"
#include "fitshandle.h"
#include "cxxutils.h"
namespace hlpBB {
Billipop::Billipop(const string fileName)
{
// Variables
planck_assert(file_present(fileName),string("missing file : ")+fileName);
Init(fileName);
_name = "Billipop";
}
Billipop::~Billipop()
{
std::cout << "Chi2 from Billipop: " << chi2 << std::endl;
}
void Billipop::Init(const string fileName)
{
std::stringstream file;
std::cout << "Using: " << fileName << std::endl;
paramfile parameters(fileName);
Lmin = parameters.find< int>( "lmin", 30 );
Lmax = parameters.find< int>( "lmax", 300);
_nell = Lmax-Lmin+1;
string DataFile = Parser::CheckPath(parameters.find<string>("Datafile"));
string CovFile = Parser::CheckPath(parameters.find<string>("Clcov"));
std::cout << "Read Dataset" << std::endl;
Billipop::ReadDataSet( DataFile);
std::cout << "Read Cl-covariance" << std::endl;
Billipop::ReadCov( Lmin, Lmax, CovFile);
// Planck global calibration factor
_n.push_back("A_planck");
}
void Billipop::ReadDataSet( string DataFile)
//read (ell,EE,BB,EB) in muK^2
{
unsigned il=0;
double var_l, var_cl;
std::ifstream fp;
_data = new HepVector(_nell);
//read data file
fp.open( DataFile.c_str());
while( fp >> var_l >> var_cl) {
if( var_l < Lmin || var_l > Lmax) continue;
(*_data)[il] = var_cl;
il++;
// std::cout << "l=" << var_l << std::endl;
}
fp.close();
if( il != _nell) {
printf( "Wrong number of l in datafile (%d)\n", (long)il);
exit(-1);
}
}
void Billipop::ReadCov( int lmin, int lmax, string CovFile)
//read covariance in muK^2 and invert
{
fitshandle * input = new fitshandle();
input->open(CovFile);
int nel = (lmax+1);
vector<long> nax = input->axes();
if( nax[0] < nel || nax[1] < nel) {
printf( "Wrong Matrix size !");
exit(-1);
}
std::cout << "init cov size: " << nax[0] << "," << nax[1] << std::endl;
arr2<double> data(nax[0],nax[1]);
input->read_image(data);
int ndata = nax[0]/2.;
_invcov = new HepSymMatrix(_nell);
for( int i=0; i<_nell; i++)
for( int j=0; j<_nell; j++)
(*_invcov)[i][j] = data[lmin+i][lmin+j];
// std::cout << *_invcov << std::endl;
std::cout << "Condition Number: " << condition(*_invcov) << std::endl;
// std::cout << "Determinant: " << _invcov->determinant() << std::endl;
//invert
int ifail;
_invcov->invert(ifail);
if (ifail !=0) {
std::cout << "Billipop cov matrix inversion fails (err="<< ifail << ")" << std::endl;
exit(-1);
}
// std::cout << *_invcov << std::endl;
}
double Billipop::computeLikelihood( const vector<unsigned int>& l,
vector<double>& cltt,
vector<double>& clte,
vector<double>& clee,
vector<double>& clbb,
vector<double>& nuisance)
//cl in muK2
{
int il=0;
HepVector model(_nell);
//Absolute Calibration
const double Aplanck=nuisance[getIndex("A_planck")];
for( int i=0; i<l.size(); i++) {
if( l[i] < Lmin || l[i] > Lmax) continue;
model[il] = clbb[i]*Aplanck*Aplanck;
il++;
}
if( il != _nell) exit(-1);
// std::ofstream fout("model_class.dat");
// for( int i=0; i<model.size(); i++)
// fout << model[i].l << " " << model[i].ee << std::endl;
// fout.close();
//compute x.C^-1.x
chi2 = (double) _invcov->similarity((*_data)-model);
_lnlkl = chi2/2.;
//return -ln(L) = chi2/2
return( chi2/2.);
}
}
//--------------------------------------------------------------------------
//
// Description:
// class Lollipop :
// base class to the LAL low-ell likelihood
//
//
// Author List:
// Matthieu Tristram (tristram@lal.in2p3.fr)
//
// History (add to end):
// creation: Avr 11 2014
// Modified: Avr 13 2015 change to oHL
//
//------------------------------------------------------------------------
#ifndef Billipop_hh
#define Billipop_hh
#include "CMB/ClLikelihood.hh"
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <sstream>
using std::vector;
namespace hlpBB {
class Billipop : public ClLikelihood
{
public:
enum cltype {TT=0,EE,BB,TE};
// Constructor
Billipop(const std::string fileName);
// Get TT maximum multipole. This will be the one use by CLASS.
int getTTmax() const {return ( (Lmax > 300) ? Lmax : 300);}
// Overidden
double computeLikelihood(const vector<unsigned int>& l,
std::vector<double>& cltt,
std::vector<double>& clte,
std::vector<double>& clee,
std::vector<double>& clbb,
std::vector<double>& nuisance);
// Use base definition
double computeLikelihood(const std::vector<unsigned int>& l,
std::vector<double>& cltt,
std::vector<double>& clte,
std::vector<double>& clee,
std::vector<double>& clbb) {
std::vector<double> nuisance;
return Billipop::computeLikelihood(l,cltt,clte,clee,clbb,nuisance);
}
// destructor
~Billipop();
private:
// Process the parameters file
void Init(const std::string fileName);
unsigned int Lmax;
unsigned int Lmin;
unsigned int _nell;
double chi2;
CLHEP::HepVector *_data;
CLHEP::HepSymMatrix *_invcov;
void ReadDataSet( string DataFile);
void ReadCov( int lmin, int lmax, string CovFile);
};
}
#endif
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