source: Sophya/trunk/Cosmo/SimLSS/cmvconchprof.cc@ 3115

#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include "timing.h"
#include "hisprof.h"

void usage(void);
void usage(void)
{
 cout<<"cmvconchprof -n namehprof file1.ppf file2.ppf ..."<<endl;
}

int main(int narg,char *arg[])
{
 string nameh = ""; // "hpkgen" "hpkgenf" "hpkrec"
 vector<string> ppfname;

 // --- Decodage arguments
 char c;
 while((c = getopt(narg,arg,"hn:")) != -1) {
  switch (c) {
  case 'n' :
    nameh = optarg;
    break;
  case 'h' :
  default :
    usage(); return -1;
  }
 }
 if(nameh.size()<=0 || optind>=narg) {usage(); return -1;}
 for (int i=optind;i<narg;i++) ppfname.push_back(arg[i]);
 cout<<"nameh = "<<nameh<<" , number of file to be read = "<<ppfname.size()<<endl;

 // --- lecture et moyenne des HProf
 HProf *hpconc = NULL;
 int nhp = 0;
 int itest = 0;
 double sum=0., sum2=0., nsum=0;
 for (int ifile=0;ifile<ppfname.size();ifile++) {
   PInPersist pis(ppfname[ifile].c_str());
   HProf hp;
   pis.GetObject(hp,nameh);
   int n = hp.NBins();
   cout<<ifile<<" : "<<ppfname[ifile]<<"  nbin="<<n<<endl;
   if(ifile==0) { 
     hpconc = new HProf(hp.XMin(),hp.XMax(),hp.NBins());
     hpconc->SetErrOpt(true);
     nhp = n;
     itest = int(hpconc->NBins()/5.+0.5);
   }
   if(n!=nhp) {cout<<"BAD NUMBER OF BINS"<<endl; continue;}
   for(int i=0;i<nhp;i++) hpconc->AddBin(i,hp(i));
   sum+=hp(itest); sum2+=hp(itest)*hp(itest), nsum++;
 }
 hpconc->UpdateHisto(true);
 if(nsum>0) {sum/=nsum; sum2/=nsum; sum2-=sum*sum;}
 cout<<"Test bin "<<itest<<" mean="<<sum<<" sigma^2="<<sum2<<"  nsum="<<nsum<<endl;
 cout<<"         "<<" mean="<<(*hpconc)(itest)<<" sigma^2="<<hpconc->Error2(itest)<<endl;

 // --- ecriture ascii
 if(1) {
   FILE * fdata = fopen("cmvconchprof.data","w");
   fprintf(fdata,"# k(Mpc^-1) <pkz_reconstruit> variance_pkz\n");
   for(int i=0;i<hpconc->NBins();i++) {
     double k = hpconc->BinCenter(i);
     fprintf(fdata,"%e %e %e\n",k,(*hpconc)(i),hpconc->Error2(i));
   }
   fclose(fdata);
 }

 // --- ecriture ppf
 string tagobs = "cmvconchprof.ppf";
 POutPersist posobs(tagobs);
 tagobs = "hpconc"; posobs.PutObject(*hpconc,tagobs);

 delete hpconc;

 return 0;
}

/*
openppf cmvconchprof.ppf

disp hpconc
n/plot hpconc.val%log10(x) x>0 ! "connectpoints"
*/