#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 "histerr.h"

void usage(void);
void usage(void)
{
 cout<<"cmvconchprof -n name_histoerr 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 HistoErr
 HistoErr *herrconc = NULL;
 int nherr = 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());
   HistoErr herr;
   pis.GetObject(herr,nameh);
   int n = herr.NBins();
   cout<<ifile<<" : "<<ppfname[ifile]<<"  nbin="<<n<<endl;
   if(ifile==0) { 
     herrconc = new HistoErr(herr.XMin(),herr.XMax(),herr.NBins());
     nherr = n;
     itest = int(herrconc->NBins()/5.+0.5);
   }
   if(n!=nherr) {cout<<"BAD NUMBER OF BINS"<<endl; continue;}
   for(int i=0;i<nherr;i++) herrconc->AddBin(i,herr(i));
   sum+=herr(itest); sum2+=herr(itest)*herr(itest), nsum++;
 }
 herrconc->ToCorrel();
 if(nsum>0) {sum/=nsum; sum2/=nsum; sum2-=sum*sum;}
 cout<<"Test bin "<<itest<<" mean="<<sum<<" sigma^2="<<sum2<<"  nsum="<<nsum<<endl;
 cout<<"         "<<" mean="<<(*herrconc)(itest)<<" sigma^2="<<herrconc->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<herrconc->NBins();i++) {
     double k = herrconc->BinCenter(i);
     fprintf(fdata,"%e %e %e\n",k,(*herrconc)(i),herrconc->Error2(i));
   }
   fclose(fdata);
 }

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

 delete herrconc;

 return 0;
}

/*
openppf cmvconchprof.ppf

disp herrconc

n/plot herrconc.val%log10(x) x>0 ! "connectpoints"
n/plot herrconc.err%log10(x) x>0 ! "connectpoints same red"

n/plot herrconc.err/val%log10(x) x>0&&val>0 ! "connectpoints"
*/