| 1 | // Utilisation de SOPHYA pour faciliter les tests ...
 | 
|---|
| 2 | #include "sopnamsp.h"
 | 
|---|
| 3 | #include "machdefs.h"
 | 
|---|
| 4 | 
 | 
|---|
| 5 | /* ---------------------------------------------------------- 
 | 
|---|
| 6 |    Projet BAORadio - (C) LAL/IRFU  2008-2010
 | 
|---|
| 7 | 
 | 
|---|
| 8 |    Programme de lecture des fichiers vecteurs/matrices de 
 | 
|---|
| 9 |    visibilites produits par mcrd / vismfib 
 | 
|---|
| 10 |    R. Ansari, C. Magneville   -  LAL/Irfu
 | 
|---|
| 11 |    V : Mai 2009
 | 
|---|
| 12 |    ---------------------------------------------------------- */
 | 
|---|
| 13 | 
 | 
|---|
| 14 | // include standard c/c++
 | 
|---|
| 15 | #include <math.h>
 | 
|---|
| 16 | #include <stdio.h>
 | 
|---|
| 17 | #include <stdlib.h>
 | 
|---|
| 18 | #include <string.h>
 | 
|---|
| 19 | 
 | 
|---|
| 20 | #include <iostream>
 | 
|---|
| 21 | #include <string>
 | 
|---|
| 22 | 
 | 
|---|
| 23 | #include "pexceptions.h"
 | 
|---|
| 24 | #include "tvector.h"
 | 
|---|
| 25 | #include "fioarr.h"
 | 
|---|
| 26 | // #include "tarrinit.h"
 | 
|---|
| 27 | #include "ntuple.h" 
 | 
|---|
| 28 | #include "datatable.h" 
 | 
|---|
| 29 | #include "histinit.h" 
 | 
|---|
| 30 | #include "matharr.h" 
 | 
|---|
| 31 | #include "timestamp.h"
 | 
|---|
| 32 | #include <utilarr.h>
 | 
|---|
| 33 | 
 | 
|---|
| 34 | // include sophya mesure ressource CPU/memoire ...
 | 
|---|
| 35 | #include "resusage.h"
 | 
|---|
| 36 | #include "ctimer.h"
 | 
|---|
| 37 | #include "timing.h"
 | 
|---|
| 38 | 
 | 
|---|
| 39 | 
 | 
|---|
| 40 | //--------------------------- Fonctions de ce fichier   ------------------- 
 | 
|---|
| 41 | int Usage(bool fgshort=true);
 | 
|---|
| 42 | 
 | 
|---|
| 43 | int DecodeProc(int narg, char* arg[]);
 | 
|---|
| 44 | int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, 
 | 
|---|
| 45 |                 vector<sa_size_t> tfrlist);
 | 
|---|
| 46 | int FillVisDTable(BaseDataTable& visdt_, TMatrix< complex<r_4> > vismtx_, TVector< uint_4> chanum_,
 | 
|---|
| 47 |                   sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_) ;
 | 
|---|
| 48 | 
 | 
|---|
| 49 | int DecodeProcVJun09(int narg, char* arg[]);
 | 
|---|
| 50 | int ProcSVFilesVJun09(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card);
 | 
|---|
| 51 | //------------------------------------------------------------------------------------------------------------
 | 
|---|
| 52 | 
 | 
|---|
| 53 | /* --Fonction-- */
 | 
|---|
| 54 | int Usage(bool fgshort)
 | 
|---|
| 55 | {
 | 
|---|
| 56 |   cout << " --- svv2mtx.cc : Read PPF files produced by mcrd/visfmib to make time-frequency\n"
 | 
|---|
| 57 |        << "     matrices and Visibilites=f(time) datatable " << endl;
 | 
|---|
| 58 |   cout << " Usage: svv2mtx -mcrd InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq [card=1] \n" 
 | 
|---|
| 59 |        << "   OR \n" 
 | 
|---|
| 60 |        << "        svv2mtx InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq VisMtxRowList" << endl; 
 | 
|---|
| 61 |   if (fgshort) {
 | 
|---|
| 62 |     cout << " svv2mtx -h for detailed instructions" << endl;
 | 
|---|
| 63 |     return 1;
 | 
|---|
| 64 |   }
 | 
|---|
| 65 |   cout << " -mcrd : Read mcrd output files -> ProcSVFilesVJun09() \n" 
 | 
|---|
| 66 |        << " InOutPath : Input/Output directory name \n" 
 | 
|---|
| 67 |        << " Imin,Imax,IStep: Input PPF files sequence number \n"
 | 
|---|
| 68 |        << "    FileNames=InOutPath/Ch12_II.fits Imin<=II<=Imax II+=IStep \n" 
 | 
|---|
| 69 |        << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for ntuple\n"
 | 
|---|
| 70 |        << " VisMtxRowList : List of visibiliy matrix rows (0,2,...) -> time-freq map \n"
 | 
|---|
| 71 |        << " card=1 Ch12 , card=2 Ch34    (mcrd output files) " << endl;
 | 
|---|
| 72 |   return 1;
 | 
|---|
| 73 | }
 | 
|---|
| 74 | 
 | 
|---|
| 75 | //----------------------------------------------------
 | 
|---|
| 76 | //----------------------------------------------------
 | 
|---|
| 77 | int main(int narg, char* arg[])
 | 
|---|
| 78 | {
 | 
|---|
| 79 |   if ((narg>1)&&(strcmp(arg[1],"-h")==0))  return Usage(false);
 | 
|---|
| 80 |   if (narg<4) return Usage(true);
 | 
|---|
| 81 |   HiStatsInitiator _inia;
 | 
|---|
| 82 |   //   TArrayInitiator  _inia;
 | 
|---|
| 83 | 
 | 
|---|
| 84 |   int rc = 0;
 | 
|---|
| 85 |   try {
 | 
|---|
| 86 |     bool fgvjun09=false;
 | 
|---|
| 87 |     if ((strcmp(arg[1],"-mcrd")==0))  fgvjun09=true;
 | 
|---|
| 88 |     ResourceUsage resu;
 | 
|---|
| 89 |     if (fgvjun09)  DecodeProcVJun09(narg-2, arg+2);
 | 
|---|
| 90 |     else  DecodeProc(narg-1, arg+1);
 | 
|---|
| 91 |     resu.Update();
 | 
|---|
| 92 |     cout << resu;
 | 
|---|
| 93 |   }
 | 
|---|
| 94 |   catch (PException& exc) {
 | 
|---|
| 95 |     cerr << " svv2mtx.cc catched PException " << exc.Msg() << endl;
 | 
|---|
| 96 |     rc = 77;
 | 
|---|
| 97 |   }  
 | 
|---|
| 98 |   catch (std::exception& sex) {
 | 
|---|
| 99 |     cerr << "\n svv2mtx.cc std::exception :" 
 | 
|---|
| 100 |          << (string)typeid(sex).name() << "\n msg= " 
 | 
|---|
| 101 |          << sex.what() << endl;
 | 
|---|
| 102 |     rc = 78;
 | 
|---|
| 103 |   }
 | 
|---|
| 104 |   catch (...) {
 | 
|---|
| 105 |     cerr << " svv2mtx.cc catched unknown (...) exception  " << endl; 
 | 
|---|
| 106 |     rc = 79; 
 | 
|---|
| 107 |   } 
 | 
|---|
| 108 | 
 | 
|---|
| 109 |   cout << ">>>> svv2mtx.cc ------- END ----------- RC=" << rc << endl;
 | 
|---|
| 110 |   return rc;
 | 
|---|
| 111 | 
 | 
|---|
| 112 | }
 | 
|---|
| 113 | 
 | 
|---|
| 114 | //-------------------------------------------------------------------- 
 | 
|---|
| 115 | //       Traitement fichiers produits par vismfib  (V Nov09)
 | 
|---|
| 116 | /* --Fonction-- */
 | 
|---|
| 117 | int DecodeProc(int narg, char* arg[])
 | 
|---|
| 118 | {
 | 
|---|
| 119 |   // Decodage des arguments et traitement 
 | 
|---|
| 120 |   string inoutpath = arg[0];    
 | 
|---|
| 121 |   int imin=0;
 | 
|---|
| 122 |   int imax=0;
 | 
|---|
| 123 |   int istep=1;
 | 
|---|
| 124 |   sscanf(arg[1],"%d,%d,%d",&imin,&imax,&istep);
 | 
|---|
| 125 |   int jf1=0;
 | 
|---|
| 126 |   int jf2=0;
 | 
|---|
| 127 |   int nfreq=0;
 | 
|---|
| 128 |   sscanf(arg[2],"%d,%d,%d",&jf1,&jf2,&nfreq);
 | 
|---|
| 129 |   int card=1;
 | 
|---|
| 130 |   vector<sa_size_t> rowlist;
 | 
|---|
| 131 |   if (narg>3) { 
 | 
|---|
| 132 |     EnumeratedSequence es;
 | 
|---|
| 133 |     int nbad;
 | 
|---|
| 134 |     es.Append(arg[3], nbad, ",");
 | 
|---|
| 135 |     for(int j=0; j<es.Size(); j++) rowlist.push_back((int)es.Value(j));
 | 
|---|
| 136 |   }
 | 
|---|
| 137 |   if (rowlist.size()<1) rowlist.push_back(0);
 | 
|---|
| 138 | 
 | 
|---|
| 139 |   cout << " ----- svv2mtx/DecodeProc - Start - InOutPath= " << inoutpath << " IMin,Max,Step=" 
 | 
|---|
| 140 |        << imin << "," << imax << "," << istep << " Card=" << card << endl;
 | 
|---|
| 141 |   cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << "  ------- " << endl;
 | 
|---|
| 142 |   cout << " RowList: " ;
 | 
|---|
| 143 |   for(int j=0; j<rowlist.size(); j++)  cout << rowlist[j] << " , " ;
 | 
|---|
| 144 |   cout << endl;
 | 
|---|
| 145 |   int rc=ProcSVFiles(inoutpath, imin, imax, istep, jf1, jf2, nfreq, rowlist);
 | 
|---|
| 146 |   return rc;
 | 
|---|
| 147 | }
 | 
|---|
| 148 | 
 | 
|---|
| 149 | // Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
 | 
|---|
| 150 | /* --Fonction-- */
 | 
|---|
| 151 | int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, 
 | 
|---|
| 152 |                 vector<sa_size_t> rowlist)
 | 
|---|
| 153 | {
 | 
|---|
| 154 |   Timer tm("ProcSVFiles");
 | 
|---|
| 155 |   
 | 
|---|
| 156 |   DataTable visdt;
 | 
|---|
| 157 |   visdt.AddDoubleColumn("mfc");
 | 
|---|
| 158 |   visdt.AddDoubleColumn("mtt");
 | 
|---|
| 159 |   visdt.AddIntegerColumn("jfreq");
 | 
|---|
| 160 |   visdt.AddIntegerColumn("numch");
 | 
|---|
| 161 |   visdt.AddFloatColumn("vre");
 | 
|---|
| 162 |   visdt.AddFloatColumn("vim");
 | 
|---|
| 163 | 
 | 
|---|
| 164 |   vector< TMatrix< complex<r_4> > > vmtf;
 | 
|---|
| 165 | 
 | 
|---|
| 166 |   if (jf1<1)  jf1=1;
 | 
|---|
| 167 |   if ((jf2<1)||(jf2<jf1))  jf2=jf1; 
 | 
|---|
| 168 |   if (nfreq<1) nfreq=1;
 | 
|---|
| 169 |   int djf=(jf2-jf1)/nfreq;
 | 
|---|
| 170 |   if (djf<1) djf=1;
 | 
|---|
| 171 | 
 | 
|---|
| 172 |   char fname[1024];
 | 
|---|
| 173 | 
 | 
|---|
| 174 |   cout << " ---- ProcSVFiles()-Begin - Reading chanum vector" << endl;
 | 
|---|
| 175 |   TVector< uint_4 > chanum;
 | 
|---|
| 176 |   {
 | 
|---|
| 177 |   sprintf(fname, "%s/chanum.ppf",inoutpath.c_str());
 | 
|---|
| 178 |   PInPersist pic(fname);
 | 
|---|
| 179 |   pic >> chanum;
 | 
|---|
| 180 |   }
 | 
|---|
| 181 | 
 | 
|---|
| 182 |   cout << " ---- ProcSVFiles()-Read chanum done " << endl;
 | 
|---|
| 183 |   sa_size_t nrows = (imax-imin+1)/istep;
 | 
|---|
| 184 |   sa_size_t kr=0;
 | 
|---|
| 185 | 
 | 
|---|
| 186 |   for(int ifile=imin; ifile<=imax; ifile+=istep) {
 | 
|---|
| 187 |     sprintf(fname, "%s/vismtx%d.ppf",inoutpath.c_str(),ifile);
 | 
|---|
| 188 |     cout << " ProcSVFiles[" << ifile << "] opening file " << fname << endl;
 | 
|---|
| 189 |     PInPersist pin(fname);
 | 
|---|
| 190 |     TMatrix< complex<r_4> > vismtx;
 | 
|---|
| 191 |     pin >> vismtx;
 | 
|---|
| 192 | 
 | 
|---|
| 193 |     if (ifile==imin) {
 | 
|---|
| 194 |       sa_size_t ncols = vismtx.NCols();
 | 
|---|
| 195 |       cout << " ProcSVFilesVJun09/Info: Output Time-Frequency matrices NRows=NFiles" 
 | 
|---|
| 196 |            << nrows << " NCols=NFreq=" << ncols << endl;
 | 
|---|
| 197 |       for(size_t j=0; j<rowlist.size(); j++)   
 | 
|---|
| 198 |         vmtf.push_back(TMatrix< complex<r_4> >(nrows, ncols) );
 | 
|---|
| 199 |     }
 | 
|---|
| 200 |     for(size_t j=0; j<rowlist.size(); j++) 
 | 
|---|
| 201 |       vmtf[j].Row(kr) = vismtx.Row(rowlist[j]);
 | 
|---|
| 202 |     kr++;
 | 
|---|
| 203 |      
 | 
|---|
| 204 | // Calcul moyenne dans des bandes en frequence 
 | 
|---|
| 205 |     FillVisDTable(visdt, vismtx, chanum, jf1, jf2, djf);
 | 
|---|
| 206 |   }
 | 
|---|
| 207 | 
 | 
|---|
| 208 |   sprintf(fname, "%s/vistfreqmtx.ppf",inoutpath.c_str());
 | 
|---|
| 209 |   cout << "ProcSVFiles: Opening file " << fname << " for writing Visi(Freq,Time) matrices" << endl;  
 | 
|---|
| 210 |   POutPersist po(fname);
 | 
|---|
| 211 |   char tagb[64];
 | 
|---|
| 212 |   for(size_t j=0; j<rowlist.size(); j++)  {
 | 
|---|
| 213 |     sprintf(tagb,"visft%d", chanum(rowlist[j]));
 | 
|---|
| 214 |     po << PPFNameTag(tagb) << vmtf[j];
 | 
|---|
| 215 |   }
 | 
|---|
| 216 |   cout << visdt;
 | 
|---|
| 217 |   sprintf(fname, "%s/svvdt.ppf",inoutpath.c_str());
 | 
|---|
| 218 |   cout << "ProcSVFile: writing visibility datatable to file " << fname << endl;
 | 
|---|
| 219 |   POutPersist pod(fname);
 | 
|---|
| 220 |   pod << visdt;
 | 
|---|
| 221 | 
 | 
|---|
| 222 |   return 0;
 | 
|---|
| 223 | }
 | 
|---|
| 224 | 
 | 
|---|
| 225 | /* --Fonction-- */
 | 
|---|
| 226 | int FillVisDTable(BaseDataTable& visdt_, TMatrix< complex<r_4> > vismtx_, TVector< uint_4> chanum_,
 | 
|---|
| 227 |                   sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_) 
 | 
|---|
| 228 | {
 | 
|---|
| 229 |   double xnt_[20];
 | 
|---|
| 230 |   xnt_[0]=(double)vismtx_.Info()["MeanFC"];  
 | 
|---|
| 231 |   xnt_[1]=(double)vismtx_.Info()["MeanTT"]/1.25e8;
 | 
|---|
| 232 | 
 | 
|---|
| 233 |   uint_4 nmean_=vismtx_.Info()["NPAQSUM"];
 | 
|---|
| 234 |   if (djf_<2) {
 | 
|---|
| 235 |     for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
 | 
|---|
| 236 |       for(sa_size_t jf=jf1_; jf<jf2_; jf++) {
 | 
|---|
| 237 |         xnt_[2]=jf;
 | 
|---|
| 238 |         xnt_[3]=chanum_(rv);
 | 
|---|
| 239 |         xnt_[4]=vismtx_(rv,jf).real()/(r_4)(nmean_);
 | 
|---|
| 240 |         xnt_[5]=vismtx_(rv,jf).imag()/(r_4)(nmean_);
 | 
|---|
| 241 |         visdt_.AddRow(xnt_);
 | 
|---|
| 242 |       }
 | 
|---|
| 243 |     }
 | 
|---|
| 244 |   }
 | 
|---|
| 245 |   else {
 | 
|---|
| 246 |     for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
 | 
|---|
| 247 |       for(sa_size_t jf=jf1_; jf<jf2_; jf+=djf_) {
 | 
|---|
| 248 |         r_4 moyreal=0.;
 | 
|---|
| 249 |         r_4 moyimag=0.;
 | 
|---|
| 250 |         sa_size_t jjfmx=jf+djf_;
 | 
|---|
| 251 |         if (jjfmx > vismtx_.NCols()) jjfmx=vismtx_.NCols();
 | 
|---|
| 252 |         for(sa_size_t jjf=jf; jjf<jjfmx; jjf++) {
 | 
|---|
| 253 |           moyreal+=vismtx_(rv,jjf).real();
 | 
|---|
| 254 |           moyimag+=vismtx_(rv,jjf).imag();
 | 
|---|
| 255 |         }
 | 
|---|
| 256 |         xnt_[2]=jf+djf_/2;
 | 
|---|
| 257 |         xnt_[3]=chanum_(rv);
 | 
|---|
| 258 |         xnt_[4]=moyreal/(r_4)(nmean_*djf_);
 | 
|---|
| 259 |         xnt_[5]=moyimag/(r_4)(nmean_*djf_);
 | 
|---|
| 260 |         visdt_.AddRow(xnt_);
 | 
|---|
| 261 |       }
 | 
|---|
| 262 |     }
 | 
|---|
| 263 |   }
 | 
|---|
| 264 |   return 0;
 | 
|---|
| 265 | }
 | 
|---|
| 266 | 
 | 
|---|
| 267 | 
 | 
|---|
| 268 | //-------------------------------------------------------------------- 
 | 
|---|
| 269 | //       Traitement fichiers produits par mcrd (V Jun09) 
 | 
|---|
| 270 | /* --Fonction-- */
 | 
|---|
| 271 | int DecodeProcVJun09(int narg, char* arg[])
 | 
|---|
| 272 | {
 | 
|---|
| 273 |   // Decodage des arguments et traitement 
 | 
|---|
| 274 |   string inoutpath = arg[0];    
 | 
|---|
| 275 |   int imin=0;
 | 
|---|
| 276 |   int imax=0;
 | 
|---|
| 277 |   int istep=1;
 | 
|---|
| 278 |   sscanf(arg[1],"%d,%d,%d",&imin,&imax,&istep);
 | 
|---|
| 279 |   int jf1=0;
 | 
|---|
| 280 |   int jf2=0;
 | 
|---|
| 281 |   int nfreq=0;
 | 
|---|
| 282 |   sscanf(arg[2],"%d,%d,%d",&jf1,&jf2,&nfreq);
 | 
|---|
| 283 |   int card=1;
 | 
|---|
| 284 |   if (narg>3) card=atoi(arg[3]); 
 | 
|---|
| 285 |   cout << " ----- svv2mtx/DecodeProcVJun09 - Start - InOutPath= " << inoutpath << " IMin,Max,Step=" 
 | 
|---|
| 286 |        << imin << "," << imax << "," << istep << " Card=" << card << endl;
 | 
|---|
| 287 |   cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << "  ------- " << endl;
 | 
|---|
| 288 |   int rc=ProcSVFilesVJun09(inoutpath, imin, imax, istep, jf1, jf2, nfreq, card);
 | 
|---|
| 289 |   return rc;
 | 
|---|
| 290 | }
 | 
|---|
| 291 | 
 | 
|---|
| 292 | 
 | 
|---|
| 293 | // Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
 | 
|---|
| 294 | /* --Fonction-- */
 | 
|---|
| 295 | int ProcSVFilesVJun09(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card)
 | 
|---|
| 296 | {
 | 
|---|
| 297 |   Timer tm("ProcSVFiles");
 | 
|---|
| 298 |   char fname[512];
 | 
|---|
| 299 | // NTuple 
 | 
|---|
| 300 |   const char* nnames[10] = {"fcsm","ttsm","jfreq","s1","s2","s12","s12re","s12im","s12phi","s12mod"};
 | 
|---|
| 301 |   NTuple nt(10, nnames);
 | 
|---|
| 302 |   double xnt[15];
 | 
|---|
| 303 |   uint_4 nmnt = 0;
 | 
|---|
| 304 |   double ms1,ms2,ms12,ms12re,ms12im,ms12phi,ms12mod;
 | 
|---|
| 305 | 
 | 
|---|
| 306 |   TMatrix<r_4> s1, s2;
 | 
|---|
| 307 |   TMatrix<r_4> v12re, v12im, v12phi,v12mod;
 | 
|---|
| 308 |   sa_size_t ncols = (imax-imin+1)/istep;
 | 
|---|
| 309 |   sa_size_t nrows = 10;
 | 
|---|
| 310 |   sa_size_t kc=0;
 | 
|---|
| 311 |   for(int ifile=imin; ifile<=imax; ifile+=istep) {
 | 
|---|
| 312 |     if (card==2) 
 | 
|---|
| 313 |       sprintf(fname, "%s/Ch34_%d.ppf",inoutpath.c_str(),ifile);
 | 
|---|
| 314 |     else 
 | 
|---|
| 315 |       sprintf(fname, "%s/Ch12_%d.ppf",inoutpath.c_str(),ifile);
 | 
|---|
| 316 |     cout << " ProcSVFilesVJun09[" << ifile << "] opening file " << fname << endl;
 | 
|---|
| 317 |     PInPersist pin(fname);
 | 
|---|
| 318 |     string tag1="specV1";
 | 
|---|
| 319 |     string tag2="specV2";
 | 
|---|
| 320 |     string tag12="visiV12";
 | 
|---|
| 321 |     if (card==2) {
 | 
|---|
| 322 |       tag1 = "specV3";
 | 
|---|
| 323 |       tag2 = "specV4";
 | 
|---|
| 324 |       tag12="visiV34";
 | 
|---|
| 325 |     }
 | 
|---|
| 326 |     TVector<r_4> sv1;
 | 
|---|
| 327 |     TVector<r_4> sv2; 
 | 
|---|
| 328 |     TVector< complex<r_4> > vv12;
 | 
|---|
| 329 |     pin >> PPFNameTag(tag1) >> sv1; 
 | 
|---|
| 330 |     pin >> PPFNameTag(tag2) >> sv2; 
 | 
|---|
| 331 |     pin >> PPFNameTag(tag12) >> vv12; 
 | 
|---|
| 332 |     if (ifile==imin) {
 | 
|---|
| 333 |       nrows = sv1.Size();
 | 
|---|
| 334 |       cout << " ProcSVFilesVJun09/Info: Output s1,s2 matrix size NRows=NFreq=" 
 | 
|---|
| 335 |            << nrows << " NCols=NFiles=" << ncols << endl;
 | 
|---|
| 336 |       s1.SetSize(nrows, ncols);
 | 
|---|
| 337 |       s2.SetSize(nrows, ncols);
 | 
|---|
| 338 |       nrows = vv12.Size();
 | 
|---|
| 339 |       cout << " ProcSVFilesVJun09/Info: Output v12 matrix size NRows=NFreq=" 
 | 
|---|
| 340 |            << nrows << " NCols=NFiles=" << ncols << endl;
 | 
|---|
| 341 |       v12re.SetSize(nrows, ncols);
 | 
|---|
| 342 |       v12im.SetSize(nrows, ncols);
 | 
|---|
| 343 |       v12phi.SetSize(nrows, ncols);
 | 
|---|
| 344 |       v12mod.SetSize(nrows, ncols);
 | 
|---|
| 345 |     }
 | 
|---|
| 346 |     s1.Column(kc) = sv1;
 | 
|---|
| 347 |     s2.Column(kc) = sv2;
 | 
|---|
| 348 |     v12re.Column(kc) = real(vv12);
 | 
|---|
| 349 |     v12im.Column(kc) = imag(vv12);
 | 
|---|
| 350 |     v12phi.Column(kc) = phase(vv12);
 | 
|---|
| 351 |     v12mod.Column(kc) = module(vv12);
 | 
|---|
| 352 | 
 | 
|---|
| 353 | // Calcul moyenne dans des bandes en frequence 
 | 
|---|
| 354 |     int deltajf=(jf2-jf1)/nfreq;
 | 
|---|
| 355 |     if (deltajf<1) deltajf=1; 
 | 
|---|
| 356 |     for(int kf=0; kf<nfreq; kf++) {
 | 
|---|
| 357 |       sa_size_t jfstart=jf1+kf*deltajf;
 | 
|---|
| 358 |       sa_size_t jfend=jfstart+deltajf;
 | 
|---|
| 359 |       if (jfend>jf2)  break;
 | 
|---|
| 360 |       nmnt=0;  ms1=ms2=ms12=ms12re=ms12im=ms12phi=ms12mod=0.;
 | 
|---|
| 361 |       for(sa_size_t jf=jfstart; jf<jfend; jf++) { 
 | 
|---|
| 362 |         ms1 += s1(jf,kc);
 | 
|---|
| 363 |         ms2 += s2(jf,kc);
 | 
|---|
| 364 |         ms12re += v12re(jf,kc);
 | 
|---|
| 365 |         ms12im += v12im(jf,kc);
 | 
|---|
| 366 |         ms12phi += v12phi(jf,kc);
 | 
|---|
| 367 |         ms12mod += v12mod(jf,kc);
 | 
|---|
| 368 |       }
 | 
|---|
| 369 |       nmnt = (jfend-jfstart);
 | 
|---|
| 370 |       if (nmnt>0)  {
 | 
|---|
| 371 |         double fnorm = (double)nmnt; 
 | 
|---|
| 372 |         xnt[0] = ((int_8)(sv1.Info()["StartFC"])+(int_8)(sv1.Info()["EndFC"]))*0.5;
 | 
|---|
| 373 |         xnt[1] = ((int_8)(sv1.Info()["StartTT"])+(int_8)(sv1.Info()["EndTT"]))*0.5;
 | 
|---|
| 374 |         xnt[2] = kf;
 | 
|---|
| 375 |         xnt[3] = ms1/fnorm;
 | 
|---|
| 376 |         xnt[4] = ms2/fnorm;
 | 
|---|
| 377 |         xnt[5] = ms12/fnorm;
 | 
|---|
| 378 |         xnt[6] = ms12re/fnorm;
 | 
|---|
| 379 |         xnt[7] = ms12im/fnorm;
 | 
|---|
| 380 |         xnt[8] = ms12phi/fnorm;
 | 
|---|
| 381 |         xnt[9] = ms12mod/fnorm;
 | 
|---|
| 382 |         nt.Fill(xnt);
 | 
|---|
| 383 |       }
 | 
|---|
| 384 |     }
 | 
|---|
| 385 |     kc++;
 | 
|---|
| 386 | 
 | 
|---|
| 387 |   }
 | 
|---|
| 388 |   if (card==2) 
 | 
|---|
| 389 |     sprintf(fname, "%s/Ch34mtx.ppf",inoutpath.c_str());
 | 
|---|
| 390 |   else 
 | 
|---|
| 391 |     sprintf(fname, "%s/Ch12mtx.ppf",inoutpath.c_str());
 | 
|---|
| 392 | 
 | 
|---|
| 393 |   cout << nt;
 | 
|---|
| 394 |   cout << "ProcSVFilesVJun09: Opening file " << fname << " for writing" << endl;  
 | 
|---|
| 395 |   POutPersist po(fname);
 | 
|---|
| 396 |   string tag1="s1";
 | 
|---|
| 397 |   string tag2="s2";
 | 
|---|
| 398 |   string tag12r="v12re";
 | 
|---|
| 399 |   string tag12i="v12im";
 | 
|---|
| 400 |   string tag12p="v12phi";
 | 
|---|
| 401 |   string tagnt="nt12";
 | 
|---|
| 402 |   if (card==2) {
 | 
|---|
| 403 |     tag1="s3";
 | 
|---|
| 404 |     tag2="s4";
 | 
|---|
| 405 |     tag12r="v34re";
 | 
|---|
| 406 |     tag12i="v34im";
 | 
|---|
| 407 |     tag12p="v34phi";
 | 
|---|
| 408 |     tagnt="nt34";
 | 
|---|
| 409 |     }
 | 
|---|
| 410 |   po <<  PPFNameTag(tag1) << s1;
 | 
|---|
| 411 |   po <<  PPFNameTag(tag2) << s2;
 | 
|---|
| 412 |   po <<  PPFNameTag(tag12r) << v12re;
 | 
|---|
| 413 |   po <<  PPFNameTag(tag12i) << v12im;
 | 
|---|
| 414 |   po <<  PPFNameTag(tag12p) << v12phi;
 | 
|---|
| 415 |   po <<  PPFNameTag(tagnt) << nt;
 | 
|---|
| 416 |   cout << "ProcSVFilesVJun09: Matrices s1, s2, v12re, v12im, v12phi, NTuple nt written to file " << fname << endl;  
 | 
|---|
| 417 |   return 0;
 | 
|---|
| 418 | }
 | 
|---|
| 419 | 
 | 
|---|
| 420 | 
 | 
|---|