Changeset 3150 in Sophya for trunk/Cosmo/SimLSS/cmvconcherr.cc

Timestamp:

Jan 18, 2007, 7:23:56 PM (19 years ago)

Author:

cmv

Message:

suite modifs structure code + bug HistoErr/2D dans ComputeSpectrum cmv 18/01/2007

File:

• trunk/Cosmo/SimLSS/cmvconcherr.cc (modified) (3 diffs)

trunk/Cosmo/SimLSS/cmvconcherr.cc

@@ -9 +9 @@
 #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<<"cmvconchprof -n name_histoerr file1.ppf file2.ppf ..."<<endl;
-}
-
+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[])
 {
- string nameh = ""; // "hpkgen" "hpkgenf" "hpkrec"
+ // "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:")) != -1) {
+ 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' :
@@ -35 +50 @@
  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
+ 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;
- int itest = 0;
+ int nherr=0, nread=0, itest=0;
  double sum=0., sum2=0., nsum=0;
  for (int ifile=0;ifile<ppfname.size();ifile++) {
@@ -46 +106 @@
    HistoErr herr;
    pis.GetObject(herr,nameh);
-   int n = herr.NBins();
-   cout<<ifile<<" : "<<ppfname[ifile]<<"  nbin="<<n<<endl;
+   cout<<ifile<<" : "<<ppfname[ifile]<<"  nbin="<<herr.NBins()<<endl;
    if(ifile==0) { 
      herrconc = new HistoErr(herr.XMin(),herr.XMax(),herr.NBins());
-     nherr = n;
+     nherr = herr.NBins();
      itest = int(herrconc->NBins()/5.+0.5);
+   } else if(nherr!=herr.NBins()) {
+     cout<<"BAD NUMBER OF BINS"<<endl;
+     continue;
    }
-   if(n!=nherr) {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->ToCorrel();
+ 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 ascii
- if(1) {
-   FILE * fdata = fopen("cmvconchprof.data","w");
-   fprintf(fdata,"# k(Mpc^-1) <pkz_reconstruit> variance_pkz\n");
+ 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 k = herrconc->BinCenter(i);
-     fprintf(fdata,"%e %e %e\n",k,(*herrconc)(i),herrconc->Error2(i));
+     double x = herrconc->BinCenter(i);
+     fprintf(fdata,"%e %e %e\n",x,(*herrconc)(i),herrconc->Error2(i));
    }
    fclose(fdata);
  }
-
- // --- ecriture ppf
- string tagobs = "cmvconchprof.ppf";
- POutPersist posobs(tagobs);
- tagobs = "herrconc"; posobs.PutObject(*herrconc,tagobs);
+ 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 0;
+ 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<nherrx;j++) herrconc->AddBin(i,j,herr(i,j));
+   sum+=herr(itestx,itesty); sum2+=herr(itestx,itesty)*herr(itestx,itesty), nsum++;
+ }
+ herrconc->ToVariance();
+ 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;
 }
 
 /*
-openppf cmvconchprof.ppf
-
-disp herrconc
+#### HistoErr 1D
+openppf cmvconcherr.ppf
+
+disp herrconc "hbincont err"
+disp herrconc "hbinerr"
+disp herrconc "hbinent"
 
 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"
+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
+
+disp herrconc2 "hbincont"
+disp herrconc2 "hbinerr"
+disp herrconc2 "hbinent"
 */