Changeset 3282 in Sophya for trunk

Timestamp:

Jul 19, 2007, 10:13:51 AM (18 years ago)

Author:

cmv

Message:

calcul covariance pour une ligne Kz ou Kt 2D cmv 19/07/2007

File:

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

trunk/Cosmo/SimLSS/cmvconcherr.cc

@@ -13 +13 @@
 
 int ObjectType(string nameh,string nameppf);
-int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp,bool do_cov);
-int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp);
+int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp=2,bool do_cov=false);
+int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp=2,int ibinkt=-1,int ibinkz=-1);
 void usage(void);
 void usage(void)
@@ -22 +22 @@
 <<"  name_histoerr: object name"<<endl
 <<"  wtyp: bit0 write ASCII, bit1 write PPF, bit2 write FITS (all=7)"<<endl
-<<"  -C : compute covariance for 1D spectrum"<<endl;
+<<"  -C : compute covariance for 1D spectrum"<<endl
+<<"  -Z ibinKt : compute covariance of line Pk(ibinKt,.) along Kz for 2D"<<endl
+<<"  -T ibinKz : compute covariance of line Pk(.,ibinKz) along Kt for 2D"<<endl
+<<endl;
 }
 
@@ -32 +35 @@
  string nameh = "";
  vector<string> ppfname;
- int Write_Type = 3;
+ int Write_Type = 2;
  bool Do_Cov = false;
+ int ibinKt=-1, ibinKz=-1;
 
  // --- Decodage arguments
  char c;
- while((c = getopt(narg,arg,"hCn:w:")) != -1) {
+ while((c = getopt(narg,arg,"hCn:w:Z:T:")) != -1) {
   switch (c) {
   case 'n' :
@@ -47 +51 @@
   case 'C' :
     Do_Cov = true;
+    break;
+  case 'Z' :
+    sscanf(optarg,"%d",&ibinKt);
+    break;
+  case 'T' :
+    sscanf(optarg,"%d",&ibinKz);
     break;
   case 'h' :
@@ -65 +75 @@
  int nread=0;
  if(obtyp==1) nread = ConcatHistoErr(nameh,ppfname,Write_Type,Do_Cov);
- else if(obtyp==2) nread = ConcatHisto2DErr(nameh,ppfname,Write_Type);
+ else if(obtyp==2) nread = ConcatHisto2DErr(nameh,ppfname,Write_Type,ibinKt,ibinKz);
  else {
    cout<<"Type of object "<<nameh<<" not decoded"<<endl;
@@ -108 +118 @@
  int nherr=0, nread=0, itest=0;
  double sum=0., sum2=0., nsum=0;
- for (int ifile=0;ifile<ppfname.size();ifile++) {
+ for (unsigned 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());
@@ -122 +133 @@
      continue;
    }
+
    if(herr.NMean()!=1) cout<<"ATTENTION: file="<<ppfname[ifile]
-                             <<" nmean="<<herr.NMean()<<endl;
+                           <<" nmean="<<herr.NMean()<<endl;
    nread++;
    for(int i=0;i<nherr;i++) herrconc->AddBin(i,herr(i));
@@ -175 +187 @@
 
 //---------------------------------------------------------------
-int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp)
+int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp,int ibinkt,int ibinkz)
 {
  if(ppfname.size()<=0) return 0;
 
  Histo2DErr *herrconc = NULL;
+ TVector<r_8> Tsum_kz; TMatrix<r_8> Tsum2_kz;
+ TVector<r_8> Tsum_kt; TMatrix<r_8> Tsum2_kt;
  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++) {
+ for (unsigned int ifile=0;ifile<ppfname.size();ifile++) {
    PInPersist pis(ppfname[ifile].c_str());
    Histo2DErr herr;
@@ -188 +202 @@
    cout<<ifile<<" : "<<ppfname[ifile]<<"  nbin="
        <<herr.NBinX()<<","<<herr.NBinY()<<endl;
+
    if(ifile==0) { 
      herrconc = new Histo2DErr(herr.XMin(),herr.XMax(),herr.NBinX()
@@ -194 +209 @@
      itestx = int(herrconc->NBinX()/5.+0.5);
      itesty = int(herrconc->NBinY()/5.+0.5);
+     if(ibinkt>=0) { // matrix de covariance pour une ligne Kz
+       if(ibinkt>=nherrx) ibinkt = nherrx-1;
+       Tsum_kz.ReSize(nherry); Tsum_kz=0.;
+       Tsum2_kz.ReSize(nherry,nherry); Tsum2_kz=0.;
+     }
+     if(ibinkz>=0) { // matrix de covariance pour une ligne Kt
+       if(ibinkz>=nherry) ibinkz = nherry-1;
+       Tsum_kt.ReSize(nherrx); Tsum_kt=0.;
+       Tsum2_kt.ReSize(nherrx,nherrx); Tsum2_kt=0.;
+     }
    } 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;
+                           <<" 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++;
+
+   if(ibinkt>=0) {
+     for(int i=0;i<Tsum_kz.Size();i++) {
+       Tsum_kz(i) += herr(ibinkt,i);
+       for(int j=0;j<Tsum_kz.Size();j++) Tsum2_kz(i,j) += herr(ibinkt,i)*herr(ibinkt,j);
+     }
+   }
+
+   if(ibinkz>=0) {
+     for(int i=0;i<Tsum_kt.Size();i++) {
+       Tsum_kt(i) += herr(i,ibinkz);
+       for(int j=0;j<Tsum_kt.Size();j++) Tsum2_kt(i,j) += herr(i,ibinkz)*herr(j,ibinkz);
+     }
+   }
+
  }
  herrconc->ToVariance();
@@ -212 +253 @@
                   <<" sigma^2="<<herrconc->Error2(itestx,itesty)<<endl;
 
+ // --- Covariance
+ if(ibinkt>=0 && nread>1) {
+   Tsum_kz /= nread;
+   Tsum2_kz /= nread;
+   for(int i=0;i<Tsum_kz.Size();i++)
+     for(int j=0;j<Tsum_kz.Size();j++) Tsum2_kz(i,j) -= Tsum_kz(i)*Tsum_kz(j);
+ }
+ if(ibinkz>=0 && nread>1) {
+   Tsum_kt /= nread;
+   Tsum2_kt /= nread;
+   for(int i=0;i<Tsum_kt.Size();i++)
+     for(int j=0;j<Tsum_kt.Size();j++) Tsum2_kt(i,j) -= Tsum_kt(i)*Tsum_kt(j);
+ }
+
  // --- ecriture
  if(wrtyp&1) {   // ecriture ASCII
@@ -221 +276 @@
    POutPersist posobs(tagobs);
    tagobs = "herrconc2"; posobs.PutObject(*herrconc,tagobs);
+   if(ibinkt>=0) {
+     char str[32];
+     sprintf(str,"kzm%d",ibinkt);
+     posobs.PutObject(Tsum_kz,str);
+     sprintf(str,"kzc%d",ibinkt);
+     posobs.PutObject(Tsum2_kz,str);
+   }
+   if(ibinkz>=0) {
+     char str[32];
+     sprintf(str,"ktm%d",ibinkz);
+     posobs.PutObject(Tsum_kt,str);
+     sprintf(str,"ktc%d",ibinkz);
+     posobs.PutObject(Tsum2_kt,str);
+   }
  }
  if(wrtyp&4) {   // ecriture FITS
@@ -233 +302 @@
 
 /*
+####
 #### HistoErr 1D
+####
 openppf cmvconcherr.ppf
 
@@ -264 +335 @@
 imag cor
 
+####
 #### Histo2DErr 2D
+####
 openppf cmvconcherr2.ppf
 
@@ -271 +344 @@
 imag herrconc2 "hbinerr"
 imag herrconc2 "hbinent"
+
+zone
+
+set mean kzm140
+set cov kzc140
+
+set mean ktm32
+set cov ktc32
+
+n/plot $mean.val%n ! ! "connectpoints"
+disp $cov
+n/plot $cov.val%c r==32 ! "connectpoints"
+
+set cor ${cov}cor
+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;
+
+disp $cor
+n/plot $cor.val%c r==32 ! "connectpoints"
 */