Merge branch 'dev' into 'master'

Dev

See merge request !13

cmt/camel_setup.sh

@@ -84,6 +84,22 @@ cmt show uses
 echo "your CAMELROOT is $CAMELROOT"
 echo "CAMEL RELATIVE PATH=$CAMELROOT/lik"
 
+# put 1 at the end of next line if you want a closer look at compilaton
+export VERBOSE=
+#echo "for more verbose compilation set VERBOSE env variable"
+
+
+export PATH=${CAMELROOT}/${CAMELCONFIG}/:${PATH}
+
+#horrible hack to find if the class version is above v2.7
+export CLASS_VERSION=$(awk '/_VERSION_/{print substr($3,3,5)}' $CLASSROOT/include/common.h)
+echo "CLASS version=$CLASS_VERSION"
+if [ $(echo $CLASS_VERSION | awk '{print(substr($0,1,3)>=2.7)}') == 1 ] ; then
+echo "CLASSV_ABOVE_2_7 defined"
+export CMTEXTRATAGS=CLASSV_ABOVE_2_7
+fi
+
+
 
 echo "------------------"
 echo "CMTSITE=$CMTSITE"
@@ -92,10 +108,3 @@ echo "C++ compiler:"
 $(cmt show macro_value cpp) --version  | head -1
 echo "with flags"
 cmt show macro_value cppflags
-
-# put 1 at the end of next line if you want a closer look at compilaton
-export VERBOSE=
-echo "for more verbose compilation set VERBOSE env variable"
-
-
-export PATH=${CAMELROOT}/${CAMELCONFIG}/:${PATH}

cmt/fix4_thetas.txt

+
+*/* SP fix 06/09/18 see https://github.com/lesgourg/class_public/issues/213
+
+if minimization fails with CLASS 
+
+1/try using H0 instead of theta_s
+2/
+
+A simple thing which could work is to try changing the two occurrences
+of the class_call_try() macro to a normal class_call(), i.e. change
+
+in source/input.c :
+
+class_call_try(input_find_root(&xzero,
+                                     &fevals,
+                                     &fzw,
+                                     errmsg),
+                     errmsg,
+                     pba->shooting_error,
+                     shooting_failed=_TRUE_);
+
+to:
+
+class_call(input_find_root(&xzero,
+                                     &fevals,
+                                     &fzw,
+                                     errmsg),
+                     errmsg,
+                     errmsg);
+
+
+and similarly for the other class_call_try().

cmt/requirements-dbug

@@ -117,3 +117,9 @@ macro_append testBAO_dependencies " MinuitFit "
 macro_append testHiLLiPOP_dependencies " MinuitFit "
 macro_append test_jla_dependencies " MinuitFit "
 
+
+#SP 7/6/19:
+#Pk signature in CLASS changed in v2.7 : need to propagate a macro
+# defined in camel_setup
+macro_append cppflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"
+macro_append cflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"

cmt/requirements-jla.generic

@@ -115,3 +115,9 @@ macro_append testBAO_dependencies " MinuitFit "
 macro_append testHiLLiPOP_dependencies " MinuitFit "
 macro_append test_jla_dependencies " MinuitFit "
 
+
+#SP 7/6/19:
+#Pk signature in CLASS changed in v2.7 : need to propagate a macro
+# defined in camel_setup
+macro_append cppflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"
+macro_append cflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"

cmt/requirements-lite.generic

@@ -60,8 +60,6 @@ macro_append  CAMEL_linkopts " -lcfitsio"
 #macro_append cppflags ' -DPICODATA=\"$(PICO_DATA)\" '
 #application testPico -group=test -s=$(CAMELROOT)/src/camel/pico/exec testPico.cc
 
-
-
 ###FROM THERE YOU DONT NEED TO TOUCH
 ###########################################################################################
 #CAMEL LIBRARIES
@@ -117,3 +115,9 @@ macro_append testBAO_dependencies " MinuitFit "
 macro_append testHiLLiPOP_dependencies " MinuitFit "
 macro_append test_jla_dependencies " MinuitFit "
 
+
+#SP 7/6/19:
+#Pk signature in CLASS changed in v2.7 : need to propagate a macro
+# defined in camel_setup
+macro_append cppflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"
+macro_append cflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"

cmt/requirements-lsst.cc

@@ -53,8 +53,6 @@ macro_append  CAMEL_linkopts " -L$(CFITSIODIR)/lib -lcfitsio "
 #application FitClassFromPico -group=exec -s=../src/camel exec/FitClassFromPico.cc
 
 
-#CLASS FastPk access
-library class_extra -no_share $(CAMELROOT)/src/class_extra/*.c
 macro_append CAMEL_linkopts " -lclass_extra "
 
 #CAMEL LIBRARIES
@@ -112,3 +110,10 @@ macro_append testHiLLiPOP_dependencies " MinuitFit "
 macro_append test_jla_dependencies " MinuitFit "
 macro_append testPico_dependencies " MinuitFit "
 
+#SP 7/6/19:
+#Pk signature in CLASS changed in v2.7 : need to propagate a macro
+# defined in camel_setup
+macro_append cppflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"
+macro_append cflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"
+library class_extra -no_share $(CAMELROOT)/src/class_extra/*.c
+

cmt/requirements.cc

@@ -102,3 +102,9 @@ macro_append testHiLLiPOP_dependencies " MinuitFit "
 macro_append test_jla_dependencies " MinuitFit "
 macro_append testPico_dependencies " MinuitFit "
 
+
+#SP 7/6/19:
+#Pk signature in CLASS changed in v2.7 : need to propagate a macro
+# defined in camel_setup
+macro_append cppflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"
+macro_append cflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"

cmt/requirements.dec

@@ -115,3 +115,9 @@ macro_append testBAO_dependencies " MinuitFit "
 macro_append testHiLLiPOP_dependencies " MinuitFit "
 macro_append test_jla_dependencies " MinuitFit "
 
+
+#SP 7/6/19:
+#Pk signature in CLASS changed in v2.7 : need to propagate a macro
+# defined in camel_setup
+macro_append cppflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"
+macro_append cflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"

cmt/requirements.nersc

@@ -98,3 +98,9 @@ macro_append testHiLLiPOP_dependencies " MinuitFit "
 macro_append test_jla_dependencies " MinuitFit "
 macro_append testPico_dependencies " MinuitFit "
 
+
+#SP 7/6/19:
+#Pk signature in CLASS changed in v2.7 : need to propagate a macro
+# defined in camel_setup
+macro_append cppflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"
+macro_append cflags "" CLASSV_ABOVE_2_7 " -DCLASSV_ABOVE_2_7"

par/examples/hlpTT_LCDM/hlpTT_LCDM.par

@@ -22,8 +22,6 @@ par   Aps217x217          nui     7.6E-05	6E-06     0.0     0.1
 par   Asz                 nui     1		0.1       0.0     10
 par   Acib                nui     1.		0.1       0.0     10
 par   AdustTT		  nui	  1	        0.1	  0.0     2
-fix   AdustPP		  nui	  0.00        0.1	  0.0     2
-fix   AdustTP		  nui	  0.00        0.1	  0.0     2
 par   Aksz                nui     0.00        1.0	  0.0     10
 par   Aszxcib             nui     0.00        1.0	  0.0     10
 

par/examples/hlpTT_LCDM/hlpTT_PS.par

@@ -37,8 +37,6 @@ par   Adusty		  nui     1.		1E-01     0.0     2
 par   Asz                 nui     1.		0.1       0.0     10
 par   Acib                nui     1.		0.1       0.0     10
 par   AdustTT		  nui	  1.	        0.1	  0.0     2
-fix   AdustPP		  nui	  0.		0.1	  0.0     2
-fix   AdustTP		  nui	  0.		0.1	  0.0     2
 par   Aksz                nui     0.		1.0	  0.0     10
 par   Aszxcib             nui     0.		1.0	  0.0     10
 

src/camel/Class/ClassEngine.cc

@@ -565,7 +565,14 @@ double ClassEngine::get_Pklin(double k, double z){
   background_tau_of_z(&ba,z,&tau);
   background_at_tau(&ba,tau,ba.long_info,ba.inter_normal, &index, pvecback);
 
+#ifdef CLASSV_ABOVE_2_7
+  double pk_cb=0.;
+  double pk_cb_ic=0.;
+  int ret=spectra_pk_at_k_and_z(&ba,&pm,&sp,k,z,&mypk,pk_ic,&pk_cb,&pk_cb_ic);
+#else
   int ret=spectra_pk_at_k_and_z(&ba,&pm,&sp,k,z,&mypk,pk_ic);
+#endif
+
   return mypk;
 }
 
@@ -577,11 +584,44 @@ double ClassEngine::get_PkNL(double k, double z){
   //cout << "pknl  after bkg "<<tau<< endl;
   background_at_tau(&ba,tau,ba.long_info,ba.inter_normal, &index, pvecback);
   //cout << "pknl  after bkg @tau "<<tau<< " z "<<z<< " k "<< k << endl;
-  int ret = spectra_pk_nl_at_k_and_z(&ba,&pm,&sp,k,z,&mypk);
 
+#ifdef CLASSV_ABOVE_2_7
+  double pk_cb_tot=0;
+  int ret = spectra_pk_nl_at_k_and_z(&ba,&pm,&sp,k,z,&mypk,&pk_cb_tot);
+#else
+  int ret = spectra_pk_nl_at_k_and_z(&ba,&pm,&sp,k,z,&mypk);
+#endif
   return mypk;
 }
 
+
+bool
+ClassEngine::get_Pkvec(const std::vector<double>& knodes, const std::vector<double>& zvec, std::vector<double>& pks)
+{
+
+if (has_PkNL()){
+//non linear case
+#ifdef CLASSV_ABOVE_2_7
+  vector<double> pk_cb_tot(pks.size());
+  return (spectra_fast_pk_at_kvec_and_zvec(&ba,&sp,const_cast<double*>(&knodes[0]),knodes.size(),const_cast<double*>(&zvec[0]),zvec.size(),&pks[0],&pk_cb_tot[0],1) == _TRUE_);
+#else
+  return (camel_spectra_fast_pk_nl_at_kvec_and_zvec(&ba,&sp,const_cast<double*>(&knodes[0]),knodes.size(),const_cast<double*>(&zvec[0]), zvec.size(),&pks[0])==_TRUE_);
+#endif
+}
+else{
+//linear case
+#ifdef CLASSV_ABOVE_2_7
+  vector<double> pk_cb_tot(pks.size());
+    return (spectra_fast_pk_at_kvec_and_zvec(&ba,&sp,const_cast<double*>(&knodes[0]),knodes.size(),const_cast<double*>(&zvec[0]),zvec.size(),&pks[0],&pk_cb_tot[0],0) == _TRUE_);
+#else
+    return (camel_spectra_fast_pk_at_kvec_and_zvec(&ba,&sp,const_cast<double*>(&knodes[0]),knodes.size(),const_cast<double*>(&zvec[0]), zvec.size(),&pks[0])==_TRUE_);
+#endif
+
+}
+
+}
+
+
 double ClassEngine::get_sigma8(double z)
 {
   double tau;
@@ -719,15 +759,6 @@ double  ClassEngine::get_DMod(double z){
 
 }
 
-bool
-ClassEngine::get_Pkvec(const std::vector<double>& knodes, const std::vector<double>& zvec, std::vector<double>& pks)
-{
-  return ( has_PkNL() ? 
-(spectra_fast_pk_nl_at_kvec_and_zvec(&ba,&sp,const_cast<double*>(&knodes[0]),knodes.size(),const_cast<double*>(&zvec[0]), zvec.size(),&pks[0])==_TRUE_) 
-: (spectra_fast_pk_at_kvec_and_zvec(&ba,&sp,const_cast<double*>(&knodes[0]),knodes.size(),const_cast<double*>(&zvec[0]), zvec.size(),&pks[0])==_TRUE_) );
-
-}
-
 bool
 ClassEngine::get_PklinNodes(std::vector<double>& knodes,std::vector<double>& pknodes){
 

src/camel/exec/writeSpectraPk.cc

@@ -228,21 +228,7 @@ int main(int argc,char** argv){
   logspace(kmin,kmax,npks,kval) ;
   }
 
-  //test fast Pk
-  Timer timer;
-  cout <<"Pk timer test..." << endl;
-  vector<double> pkvec(kval.size());
-  for (size_t i=0;i<100;i++)
-    e->get_Pkvec(kval,z0,pkvec);
-  cout << "fast=" << timer.partial() <<endl;
-
-  for (size_t i=0;i<100;i++){
-    for (int k =0; k<kval.size() ; k++){
-      e->get_Pklin(kval[k], 0);
-    }
-  }
-  cout << "std=" << timer.partial() <<endl;
-
+ 
   vector<double>  pklin;
   vector<double>  pknl ;
   for (int k =0; k<kval.size() ; k++){

src/camel/test/testKlass.cc

 
 //KLASS
 #include"Class/ClassEngine.hh"
-
+#include "Timer.hh"
 #include"Parser.hh"
 #include <iostream>
 #include<string>
@@ -58,6 +58,9 @@ int main(int argc,char** argv){
   pars.add("z_pk","0,0.57");
   pars.add("z_max_pk",1);
 
+  //test halofit
+  //pars.add("non linear","halofit");
+
 
   //the class calculator (engine) 
   ClassEngine* engine(0);
@@ -70,10 +73,9 @@ int main(int argc,char** argv){
     engine=new ClassEngine(pars);
   }
   cout << "----Engine creating with parameters: " << endl;
-  //engine->printFC();
-
-  cout << "----Printing some keyword values: " << endl;
+  engine->printFC();
 
+  cout << "----Check some keywords: " << endl;
 
   //test keyword access
   const char* derived_par[]={ "H0","H(0)","z_reio","tau_reio","z_drag","rs_drag","sigma8","Dv(0.57)","Da(0.57)","H(0.57)","sigma8(0.57)", "f(0.57)","growthD(0.57)"};
@@ -89,21 +91,43 @@ int main(int argc,char** argv){
 
   //interp
   for (size_t i=0;i<Nz;i++){
-    cout << "z=" << zi[i] << " growthD=" << engine->get_growthD(zi[i]) <<endl;
+    cout << "z=" << zi[i] << "\tchi=" << engine->com_distance(zi[i]) << " Mpc\tgrowthD=" << engine->get_growthD(zi[i]) <<endl;
   }
    
   cout << "----Testing Pk output: " << endl;
+  
+  cout << "is there a Pk spectrum? : " <<engine->has_Pklin()<<endl;
+  cout << "is it non linear? : " << engine->has_PkNL()<<endl;
+ 
+
   //test pk
   double kval[]={1e-4,3e-4,1e-3,2e-3,1e-2,3e-2,0.1};
   const vector<double> kvec(kval,kval+sizeof(kval)/sizeof(double));
   vector<double> pkvec(kvec.size());
 
-  engine->get_Pkvec(kvec,0.,pkvec);
+
+   //test fast Pk is fast
+  Timer timer;
+  cout <<"Pk timing test..." << endl;
+
+  for (size_t i=0;i<100000;i++){
+    for (size_t k =0; k<kvec.size() ; k++){
+      engine->get_Pklin(kval[k], 0);
+    }
+  }
+  cout << "std=" << timer.partial() <<endl;
+  for (size_t i=0;i<100000;i++)
+    engine->get_Pkvec(kvec,0.,pkvec);
+  cout << "fast=" << timer.partial() <<endl;
+  cout << "----" << endl;
+
+  //compare values
   for (size_t i=0;i<kvec.size();i++) 
-    cout << "k=" << kvec[i] << "\t p(k)=" << pkvec[i] << endl;
+    cout << "k=" << kvec[i] << "\t p(k)=" << engine->get_Pklin(kvec[i],0.) << "\t fast=" << pkvec[i] << endl;
   //clean
 
-  cout << "----Testing update: " << endl;
+
+  cout << "----Testing params update: " << endl;
   const int i_ns=pars.findIndex("n_s");
   cout <<" n_s index=" << i_ns <<endl;
   pars.updateVal(i_ns,"1");

src/class_extra/spectra_extra.c

 
 #include "spectra.h"
 
+/*   stuf for class below v2.7*/
 
-/*   extra stuf */
-
-int spectra_fast_pk_at_kvec_and_zvec(
+#ifndef CLASSV_ABOVE_2_7
+int camel_spectra_fast_pk_at_kvec_and_zvec(
                     struct background * pba,
                     struct spectra * psp,
 		    double * kvec,
@@ -112,7 +112,7 @@ int spectra_fast_pk_at_kvec_and_zvec(
 
 
 
-int spectra_fast_pk_nl_at_kvec_and_zvec(
+int camel_spectra_fast_pk_nl_at_kvec_and_zvec(
                     struct background * pba,
                     struct spectra * psp,
 		    double * kvec,
@@ -217,3 +217,5 @@ int spectra_fast_pk_nl_at_kvec_and_zvec(
 	      
   return _SUCCESS_;
 }
+
+#endif

src/class_extra/spectra_extra.h

@@ -11,7 +11,7 @@ extern "C" {
 #endif
 /*   extra stuf */
 
-  int spectra_fast_pk_at_kvec_and_zvec(
+  int camel_spectra_fast_pk_at_kvec_and_zvec(
                     struct background * pba,
                     struct spectra * psp,
 		    double * kvec,
@@ -20,7 +20,7 @@ extern "C" {
 		    int zvec_size,
                     double * pk_tot_out);
 
- int spectra_fast_pk_nl_at_kvec_and_zvec(
+ int camel_spectra_fast_pk_nl_at_kvec_and_zvec(
                     struct background * pba,
                     struct spectra * psp,
 		    double * kvec,

work/batch/cc/runMCMC_class.sh

@@ -64,6 +64,8 @@ fi
 cat > mc_setup <<EOF
 length=$NSAMPLES
 proposal_cov=$cov
+derived=sigma8 
+do_mPK=true
 EOF
 
 echo "**************************************"

work/ref/README_tests

+
+#check CLASS versions:
+writeSpectraPk planck_nl.cml 2000 tt.fits
++ compspec.py
+
+#check camep imp of CLASS
+testKlass
+
+test_jla : JLA likelihood evaluation: 703.525
+
+
+writeChi2 hlpTT_PS_H0.par:      chi2=11055.1 (sBBN_2017_marcucci.dat)
+
+batch/cc/Minimize5 hlpTT_PS_H0.par 1-10
+
+
+v2.7.2:
+
+||0||    H0        	||     66.816	||    0.61935	||  limited 	||
+||1||    omega_b   	||   0.022112	|| 0.00081095	||  limited 	||
+||2||    omega_cdm 	||    0.12108	|| 0.00039075	||  limited 	||
+||3||    ln10^{10}A_s	||     3.0602	||  0.0045976	||  limited 	||
+||4||    n_s       	||    0.96445	||   0.010725	||  limited 	||
+||5||    tau_reio  	||    0.06295	|| 0.00083529	||  limited 	||
+||6||    A_planck  	||          1	||  0.0015371	||  limited 	||
+||7||    c0        	|| 0.00194	|| 0.00116	||  limited 	||
+||8||    c1        	|| 0.0021	|| 0.00116	||  limited 	||
+||9||    c2        	||     0	||     0	||   const  	||
+||10||    c3        	|| -0.00218	|| 0.00101	||  limited 	||
+||11||    c4        	|| 0.00207	|| 0.00133	||  limited 	||
+||12||    c5        	|| 0.0022	|| 0.00132	||  limited 	||
+||13||    Aradio    	||  1.63	|| 0.0744	||  limited 	||
+||14||    Adusty    	|| 0.789	|| 0.0427	||  limited 	||
+||15||    Asz       	||  1.12	|| 0.143	||  limited 	||
+||16||    Acib      	|| 0.745	|| 0.051	||  limited 	||
+||17||    AdustTT   	|| 0.956	|| 0.0486	||  limited 	||
+||18||    Aksz      	|| 8.33e-05	|| 0.499	||  limited 	||
+||19||    Aszxcib   	||  1.98	|| 0.255	||  limited 	||
+WARNING: FunctionMinimum is invalid.
+===== Hillipop:DX11dHM_superExt_CO_TT.lik: chi2=9997.42
+===== CMB_all: chi2=9997.42
+===== Gauss1 constraint on tau_reio(0.058+-0.012): chi2=0.170144
+===== Gauss1 constraint on tau_reio(0.058+-0.012): chi2=0.170144
+===== Gauss1 constraint on AdustTT(1+-0.2): chi2=0.0485947
+===== Gauss1 constraint on AdustTT(1+-0.2): chi2=0.0485947
+===== Gauss1 constraint on Asz(1+-0.2): chi2=0.377042
+===== Gauss1 constraint on Asz(1+-0.2): chi2=0.377042
+===== Gauss1 constraint on Acib(1+-0.2): chi2=1.63119
+===== Gauss1 constraint on Acib(1+-0.2): chi2=1.63119
+===== Gauss1 constraint on Aradio(1+-0.2): chi2=9.8846
+===== Gauss1 constraint on Aradio(1+-0.2): chi2=9.8846
+===== Gauss1 constraint on Adusty(1+-0.2): chi2=1.10816
+===== Gauss1 constraint on Adusty(1+-0.2): chi2=1.10816
+===== Gauss1 constraint on A_planck(1+-0.0025): chi2=0.000122272
+===== Gauss1 constraint on A_planck(1+-0.0025): chi2=0.000122272
+===== Gauss1 constraint on c0(0+-0.002): chi2=0.939762
+===== Gauss1 constraint on c0(0+-0.002): chi2=0.939762
+===== Gauss1 constraint on c1(0+-0.002): chi2=1.10743
+===== Gauss1 constraint on c1(0+-0.002): chi2=1.10743
+===== Gauss1 constraint on c3(0+-0.002): chi2=1.18895
+===== Gauss1 constraint on c3(0+-0.002): chi2=1.18895
+===== Gauss1 constraint on c4(0.002+-0.002): chi2=0.00123482
+===== Gauss1 constraint on c4(0.002+-0.002): chi2=0.00123482
+===== Gauss1 constraint on c5(0.002+-0.002): chi2=0.00958915
+===== Gauss1 constraint on c5(0.002+-0.002): chi2=0.00958915
+TIMER TOTAL TIME=19676.4 s	=327.94 min	=5.46567h
+===== Chi2Combiner: chi2=10013.9
+chi2min= 10013.88797181
+#iter= 3290 	tot time= 5.46567  h	time/iter= 5.98067 s
+#

work/ref/class_spectra/compspec.py

+#########################################################################
+#define some CI tests for class spectra 
+# run : ../../amd64_sl7/writeSpectraPk planck_nl.cml 2000 tt.fits
+#
+# S.Plaszczynski 07/06/19
+#########################################################################
+
+from pylab import *
+from astropy.io import fits
+
+def mrdfits(fn,hdu,header=False,silent=False):
+    hdulist = fits.open(fn)
+    #hdulist.info()
+    
+    h=hdulist[hdu].header
+    if not silent:
+        print(str(h).strip())
+    data=hdulist[hdu].data
+    hdulist.close()
+    if header:
+        return (data,h)
+    else:
+        return data
+
+
+f1="spectra_class_2.6.3.fits"
+f2="tt.fits"
+
+t1=mrdfits(f1,1)
+t2=mrdfits(f2,1)
+
+l=arange(len(t1.TT))
+llp=l*(l+1)/(2*pi)*1e12
+
+figure()
+vars=['TT','TE','EE','BB']
+i=1
+for var in vars:
+    subplot(2,2,i)
+    plot(l,llp*t1[var],'r')
+    plot(l,llp*t2[var],'k',label=var)
+    legend()
+    axhline(0,c='k',lw=0.5)
+    i+=1
+show()
+
+#lensing
+t1=mrdfits(f1,2)
+t2=mrdfits(f2,2)
+
+figure()
+vars=['PP','TP','EP']
+i=1
+for var in vars:
+    subplot(2,2,i)
+    plot(l,llp**2*t1[var],'r')
+    plot(l,llp**2*t2[var],'k',label=var)
+    legend()
+    axhline(0,c='k',lw=0.5)
+    i+=1
+show()
+
+#pk(z=0)
+t1=mrdfits(f1,3)
+t2=mrdfits(f2,3)
+figure()
+semilogx(t1.k,t1.Pklin,'r',label='lin')
+plot(t1.k,t1.Pknl,'r--',label='NL')
+
+plot(t2.k,t2.Pklin,'k')
+plot(t2.k,t2.Pknl,'k--')
+
+show()
+
+

work/ref/classcomp.scr

-Tcmb= 2.7255
-
-file1='cl_classv2.4.3.fits'
-t1=mrdfits(file1,1)
-
-file2='cl_classv2.5.0_nl.fits'
-t2=mrdfits(file2,1)
-
-file3='plikTT_bflike_pico.fits'
-;file3='cl_classv2.4.3_nl.fits'
-t3=mrdfits(file3,1)
-
-;camb
-readcol,'base_plikHM_TT_lowTEB.minimum.theory_cl',ell,dl,skip=1
-
-;class
-;file3='output/plikbf_cl_lensed.dat'
-;readcol,file3,ell,dltt,skip=7
-
-
-;restart @lmin=2
-l=dindgen(n_elements(t1.tt))
-llp=l*(l+1)/(2*!dpi)
-
-dl1=llp*t1.tt*1e12
-dl1=dl1[2:*]
-
-dl2=llp*t2.tt*1e12
-dl2=dl2[2:*]
-
-dl3=llp*t3.tt*1e12