source: Sophya/trunk/Cosmo/SimLSS/cmvconcherr.cc@ 3157

#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"
#include "hist2err.h"
#include "fitshisterr.h"

int ObjectType(string nameh,string nameppf);
int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp);
int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp);
void usage(void);
void usage(void)
{
cout
<<"cmvconcherr -w wtyp -n name_histoerr file1.ppf file2.ppf ..."<<endl
<<"  name_histoerr: object name"<<endl
<<"  wtyp: bit0 write ASCII, bit1 write PPF, bit2 write FITS (all=7)"<<endl;
}


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

 // --- Decodage arguments
 char c;
 while((c = getopt(narg,arg,"hn:w:")) != -1) {
  switch (c) {
  case 'n' :
    nameh = optarg;
    break;
  case 'w' :
    sscanf(optarg,"%d",&Write_Type);
    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()
     <<" , write_type="<<Write_Type<<endl;

 // --- Quel type d'objet ?
 int obtyp = ObjectType(nameh,ppfname[0]);

 // --- lecture et moyenne des HistoErr ou Histo2DErr
 int nread=0;
 if(obtyp==1) nread = ConcatHistoErr(nameh,ppfname,Write_Type);
 else if(obtyp==2) nread = ConcatHisto2DErr(nameh,ppfname,Write_Type);
 else {
   cout<<"Type of object "<<nameh<<" not decoded"<<endl;
   return -1;
 }
 cout<<"Number of file read: "<<nread<<endl;
 return 0;
}

//---------------------------------------------------------------
int ObjectType(string nameh,string nameppf)
// retourne: 1 si HistoErr
//           2 si Histo2DErr,
//           0 si autre objet
//          -1 si il n'y a pas d'objet de nom "nameh"
{
 PInPersist pis(nameppf.c_str());

 bool found_tag = pis.GotoNameTag(nameh.c_str());
 if(!found_tag) return -1;

 PPersist *pps = pis.ReadObject();
 AnyDataObj *obj = pps->DataObj();
 cout<<"Object type read from input PPF file : "<<typeid(*obj).name()<< endl;

 HistoErr *herr = dynamic_cast<HistoErr *>(obj);
 if(herr) return 1;

 Histo2DErr *herr2 = dynamic_cast<Histo2DErr *>(obj);
 if(herr2) return 2;

 return 0;
}

//---------------------------------------------------------------
int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp)
{
 if(ppfname.size()<=0) return 0;

 HistoErr *herrconc = NULL;
 int nherr=0, nread=0, 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);
   cout<<ifile<<" : "<<ppfname[ifile]<<"  nbin="<<herr.NBins()<<endl;
   if(ifile==0) { 
     herrconc = new HistoErr(herr.XMin(),herr.XMax(),herr.NBins());
     nherr = herr.NBins();
     itest = int(herrconc->NBins()/5.+0.5);
   } else if(nherr!=herr.NBins()) {
     cout<<"BAD NUMBER OF BINS"<<endl;
     continue;
   }
   if(herr.NMean()!=1) cout<<"ATTENTION: file="<<ppfname[ifile]
                             <<" nmean="<<herr.NMean()<<endl;
   nread++;
   for(int i=0;i<nherr;i++) herrconc->AddBin(i,herr(i));
   sum+=herr(itest); sum2+=herr(itest)*herr(itest); nsum++;
 }
 herrconc->ToVariance();
 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
 if(wrtyp&1) {   // ecriture ASCII
   FILE * fdata = fopen("cmvconcherr.data","w");
   fprintf(fdata,"# x <mean> variance\n");
   fprintf(fdata,"%.17e %.17e %d %d\n",herrconc->XMin(),herrconc->XMax()
          ,herrconc->NBins(),herrconc->NMean());
   for(int i=0;i<herrconc->NBins();i++) {
     double x = herrconc->BinCenter(i);
     fprintf(fdata,"%e %e %e\n",x,(*herrconc)(i),herrconc->Error2(i));
   }
   fclose(fdata);
 }
 if(wrtyp&2) {   // ecriture PPF
   string tagobs = "cmvconcherr.ppf";
   POutPersist posobs(tagobs);
   tagobs = "herrconc"; posobs.PutObject(*herrconc,tagobs);
 }
 if(wrtyp&4) {   // ecriture FITS
   FitsInOutFile fio("!cmvconcherr.fits",FitsInOutFile::Fits_Create);
   fio << *herrconc;
 }

 delete herrconc;

 return nread;
}

//---------------------------------------------------------------
int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp)
{
 if(ppfname.size()<=0) return 0;

 Histo2DErr *herrconc = NULL;
 int nherrx=0, nherry=0, nread=0, itestx=0, itesty=0;
 double sum=0., sum2=0., nsum=0;
 for (int ifile=0;ifile<ppfname.size();ifile++) {
   PInPersist pis(ppfname[ifile].c_str());
   Histo2DErr herr;
   pis.GetObject(herr,nameh);
   cout<<ifile<<" : "<<ppfname[ifile]<<"  nbin="
       <<herr.NBinX()<<","<<herr.NBinY()<<endl;
   if(ifile==0) { 
     herrconc = new Histo2DErr(herr.XMin(),herr.XMax(),herr.NBinX()
                              ,herr.YMin(),herr.YMax(),herr.NBinY());
     nherrx = herr.NBinX(); nherry = herr.NBinY();
     itestx = int(herrconc->NBinX()/5.+0.5);
     itesty = int(herrconc->NBinY()/5.+0.5);
   } else if(nherrx!=herr.NBinX() || nherry!=herr.NBinY()) {
     cout<<"BAD NUMBER OF BINS"<<endl;
     continue;
   }
   if(herr.NMean()!=1) cout<<"ATTENTION: file="<<ppfname[ifile]
                             <<" nmean="<<herr.NMean()<<endl;
   nread++;
   for(int i=0;i<nherrx;i++)
     for(int j=0;j<nherry;j++) herrconc->AddBin(i,j,herr(i,j));
   sum+=herr(itestx,itesty); sum2+=herr(itestx,itesty)*herr(itestx,itesty); nsum++;
 }
 herrconc->ToVariance();
 if(nsum>0) {sum/=nsum; sum2/=nsum; sum2-=sum*sum;}
 cout<<"Test bin "<<itestx<<","<<itesty
     <<" mean="<<sum<<" sigma^2="<<sum2<<"  nsum="<<nsum<<endl;
 cout<<"         "<<" mean="<<(*herrconc)(itestx,itesty)
                  <<" sigma^2="<<herrconc->Error2(itestx,itesty)<<endl;

 // --- ecriture
 if(wrtyp&1) {   // ecriture ASCII
   cout<<"ERROR: No ASCII writing implemented"<<endl;
 }
 if(wrtyp&2) {   // ecriture PPF
   string tagobs = "cmvconcherr2.ppf";
   POutPersist posobs(tagobs);
   tagobs = "herrconc2"; posobs.PutObject(*herrconc,tagobs);
 }
 if(wrtyp&4) {   // ecriture FITS
   FitsInOutFile fio("!cmvconcherr2.fits",FitsInOutFile::Fits_Create);
   fio << *herrconc;
 }

 delete herrconc;

 return nread;
}

/*
#### HistoErr 1D
openppf cmvconcherr.ppf

zone 2 2
disp herrconc "hbincont err"
disp herrconc "hbinerr"
disp herrconc "hbinent"

zone
n/plot herrconc.val%log10(x) x>0 ! "connectpoints"
n/plot herrconc.sqrt(err2)%log10(x) x>0&&err2>0 ! "connectpoints same red"

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

#### Histo2DErr 2D
openppf cmvconcherr2.ppf

zone 2 2
imag herrconc2 "hbincont"
imag herrconc2 "hbinerr"
imag herrconc2 "hbinent"
*/