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

Timestamp:

Mar 21, 2007, 7:18:25 PM (19 years ago)

Author:

cmv

Message:

petis changements cmv 21/03/2007

File:

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

trunk/Cosmo/SimLSS/cmvconcherr.cc

@@ -13 +13 @@
 
 int ObjectType(string nameh,string nameppf);
-int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp);
+int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp,bool do_cov);
 int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp);
 void usage(void);
@@ -21 +21 @@
 <<"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;
+<<"  wtyp: bit0 write ASCII, bit1 write PPF, bit2 write FITS (all=7)"<<endl
+<<"  -C : compute covariance for 1D spectrum"<<endl;
 }
 
@@ -32 +33 @@
  vector<string> ppfname;
  int Write_Type = 3;
+ bool Do_Cov = false;
 
  // --- Decodage arguments
  char c;
- while((c = getopt(narg,arg,"hn:w:")) != -1) {
+ while((c = getopt(narg,arg,"hCn:w:")) != -1) {
   switch (c) {
   case 'n' :
@@ -43 +45 @@
     sscanf(optarg,"%d",&Write_Type);
     break;
+  case 'C' :
+    Do_Cov = true;
+    break;
   case 'h' :
   default :
@@ -59 +64 @@
  // --- lecture et moyenne des HistoErr ou Histo2DErr
  int nread=0;
- if(obtyp==1) nread = ConcatHistoErr(nameh,ppfname,Write_Type);
+ if(obtyp==1) nread = ConcatHistoErr(nameh,ppfname,Write_Type,Do_Cov);
  else if(obtyp==2) nread = ConcatHisto2DErr(nameh,ppfname,Write_Type);
  else {
@@ -95 +100 @@
 
 //---------------------------------------------------------------
-int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp)
+int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp,bool do_cov)
 {
  if(ppfname.size()<=0) return 0;
 
  HistoErr *herrconc = NULL;
+ TVector<r_8> Tsum; TMatrix<r_8> Tsum2;
  int nherr=0, nread=0, itest=0;
  double sum=0., sum2=0., nsum=0;
@@ -111 +117 @@
      nherr = herr.NBins();
      itest = int(herrconc->NBins()/5.+0.5);
+     if(do_cov) {Tsum.ReSize(nherr); Tsum=0.; Tsum2.ReSize(nherr,nherr); Tsum2=0.;}
    } else if(nherr!=herr.NBins()) {
      cout<<"BAD NUMBER OF BINS"<<endl;
@@ -120 +127 @@
    for(int i=0;i<nherr;i++) herrconc->AddBin(i,herr(i));
    sum+=herr(itest); sum2+=herr(itest)*herr(itest); nsum++;
+
+   if(do_cov) {
+     for(int i=0;i<nherr;i++) {
+       Tsum(i) += herr(i);
+       for(int j=0;j<nherr;j++) Tsum2(i,j) += herr(i)*herr(j);
+     }
+   }
+
  }
  herrconc->ToVariance();
@@ -127 +142 @@
                   <<" sigma^2="<<herrconc->Error2(itest)<<endl;
 
+ // --- Covariance
+ if(do_cov && nread>1) {
+   Tsum /= nread;
+   Tsum2 /= nread;
+   for(int i=0;i<nherr;i++)
+     for(int j=0;j<nherr;j++) Tsum2(i,j) -= Tsum(i)*Tsum(j);
+ }
+
  // --- ecriture
  if(wrtyp&1) {   // ecriture ASCII
@@ -132 +155 @@
    herrconc->WriteASCII(asname);
  }
- if(wrtyp&2) {   // ecriture PPF
+ if(wrtyp&2 || do_cov) {   // ecriture PPF
    string tagobs = "cmvconcherr.ppf";
    POutPersist posobs(tagobs);
    tagobs = "herrconc"; posobs.PutObject(*herrconc,tagobs);
+   if(do_cov) {
+     tagobs = "mean"; posobs.PutObject(Tsum,tagobs);
+     tagobs = "cov"; posobs.PutObject(Tsum2,tagobs);
+   }
  }
  if(wrtyp&4) {   // ecriture FITS
@@ -220 +247 @@
 n/plot herrconc.sqrt(err2)/val%log10(x) x>0&&err2>0&&val>0 ! "connectpoints"
 
+disp mean
+imag cov
+
+del cor
+c++exec \
+TMatrix<r_8> cor(cov,false); cor = 0.; \
+for(int i=0;i<cor.NRows();i++) { \
+  for(int j=0;j<cor.NCols();j++) { \
+    double v = cov(i,i)*cov(j,j); \
+    if(v<=0.) continue; \
+    cor(i,j) = cov(i,j)/sqrt(v); \
+} } \
+KeepObj(cor); \
+cout<<"Matrice cor cree "<<endl;
+
+imag cor
+
 #### Histo2DErr 2D
 openppf cmvconcherr2.ppf