| [3646] | 1 | // Utilisation de SOPHYA pour faciliter les tests ...
|
|---|
| 2 | #include "sopnamsp.h"
|
|---|
| 3 | #include "machdefs.h"
|
|---|
| 4 |
|
|---|
| 5 | /* ----------------------------------------------------------
|
|---|
| 6 | Programme de lecture multi canaux pour BAORadio
|
|---|
| 7 | R. Ansari, C. Magneville
|
|---|
| 8 | V : Mai 2009
|
|---|
| 9 | ---------------------------------------------------------- */
|
|---|
| 10 |
|
|---|
| 11 | // include standard c/c++
|
|---|
| 12 | #include <math.h>
|
|---|
| 13 | #include <stdio.h>
|
|---|
| 14 | #include <stdlib.h>
|
|---|
| 15 | #include <string.h>
|
|---|
| 16 |
|
|---|
| 17 | #include <iostream>
|
|---|
| 18 | #include <string>
|
|---|
| 19 |
|
|---|
| 20 | #include "pexceptions.h"
|
|---|
| 21 | #include "tvector.h"
|
|---|
| 22 | #include "fioarr.h"
|
|---|
| 23 | // #include "tarrinit.h"
|
|---|
| 24 | #include "ntuple.h"
|
|---|
| 25 | #include "histinit.h"
|
|---|
| 26 | #include "matharr.h"
|
|---|
| 27 | #include "timestamp.h"
|
|---|
| 28 | /*
|
|---|
| 29 | #include "fftpserver.h"
|
|---|
| 30 | #include "fftwserver.h"
|
|---|
| 31 |
|
|---|
| 32 | #include "FFTW/fftw3.h"
|
|---|
| 33 | */
|
|---|
| 34 | // include sophya mesure ressource CPU/memoire ...
|
|---|
| 35 | #include "resusage.h"
|
|---|
| 36 | #include "ctimer.h"
|
|---|
| 37 | #include "timing.h"
|
|---|
| 38 |
|
|---|
| 39 |
|
|---|
| 40 | int Usage(bool fgshort=true);
|
|---|
| 41 | // Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
|
|---|
| [3647] | 42 | int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card=1);
|
|---|
| [3646] | 43 |
|
|---|
| 44 | //----------------------------------------------------
|
|---|
| 45 | //----------------------------------------------------
|
|---|
| 46 | int main(int narg, char* arg[])
|
|---|
| 47 | {
|
|---|
| 48 | if ((narg>1)&&(strcmp(arg[1],"-h")==0)) return Usage(false);
|
|---|
| 49 | if (narg<4) return Usage(true);
|
|---|
| 50 |
|
|---|
| 51 | HiStatsInitiator _inia;
|
|---|
| 52 | // TArrayInitiator _inia;
|
|---|
| 53 |
|
|---|
| 54 | int rc = 0;
|
|---|
| 55 | try {
|
|---|
| 56 | string inoutpath = arg[1];
|
|---|
| 57 | int imin=0;
|
|---|
| 58 | int imax=0;
|
|---|
| 59 | int istep=1;
|
|---|
| 60 | sscanf(arg[2],"%d,%d,%d",&imin,&imax,&istep);
|
|---|
| 61 | int jf1=0;
|
|---|
| 62 | int jf2=0;
|
|---|
| [3647] | 63 | int nfreq=0;
|
|---|
| 64 | sscanf(arg[3],"%d,%d,%d",&jf1,&jf2,&nfreq);
|
|---|
| [3694] | 65 | int card=1;
|
|---|
| 66 | if (narg>4) card=atoi(arg[4]);
|
|---|
| [3646] | 67 | cout << " ----- svv2mtx.cc Start - InOutPath= " << inoutpath << " IMin,Max,Step="
|
|---|
| [3694] | 68 | << imin << "," << imax << "," << istep << " Card=" << card << endl;
|
|---|
| 69 | cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << " ------- " << endl;
|
|---|
| [3646] | 70 | ResourceUsage resu;
|
|---|
| [3694] | 71 | rc=ProcSVFiles(inoutpath, imin, imax, istep, jf1, jf2, nfreq, card);
|
|---|
| [3646] | 72 | resu.Update();
|
|---|
| 73 | cout << resu;
|
|---|
| 74 | }
|
|---|
| 75 | catch (PException& exc) {
|
|---|
| 76 | cerr << " svv2mtx.cc catched MiniFITSException " << exc.Msg() << endl;
|
|---|
| 77 | rc = 77;
|
|---|
| 78 | }
|
|---|
| 79 | catch (std::exception& sex) {
|
|---|
| 80 | cerr << "\n svv2mtx.cc std::exception :"
|
|---|
| 81 | << (string)typeid(sex).name() << "\n msg= "
|
|---|
| 82 | << sex.what() << endl;
|
|---|
| 83 | rc = 78;
|
|---|
| 84 | }
|
|---|
| 85 | catch (...) {
|
|---|
| 86 | cerr << " svv2mtx.cc catched unknown (...) exception " << endl;
|
|---|
| 87 | rc = 79;
|
|---|
| 88 | }
|
|---|
| 89 |
|
|---|
| 90 | cout << ">>>> svv2mtx.cc ------- END ----------- RC=" << rc << endl;
|
|---|
| 91 | return rc;
|
|---|
| 92 |
|
|---|
| 93 | }
|
|---|
| 94 |
|
|---|
| 95 |
|
|---|
| 96 |
|
|---|
| 97 | // Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
|
|---|
| [3647] | 98 | int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card)
|
|---|
| [3646] | 99 | {
|
|---|
| 100 | Timer tm("ProcSVFiles");
|
|---|
| 101 | char fname[512];
|
|---|
| 102 | // NTuple
|
|---|
| [3647] | 103 | const char* nnames[10] = {"fcsm","ttsm","jfreq","s1","s2","s12","s12re","s12im","s12phi","s12mod"};
|
|---|
| 104 | NTuple nt(10, nnames);
|
|---|
| 105 | double xnt[15];
|
|---|
| [3646] | 106 | uint_4 nmnt = 0;
|
|---|
| [3647] | 107 | double ms1,ms2,ms12,ms12re,ms12im,ms12phi,ms12mod;
|
|---|
| [3646] | 108 |
|
|---|
| 109 | TMatrix<r_4> s1, s2;
|
|---|
| [3647] | 110 | TMatrix<r_4> v12re, v12im, v12phi,v12mod;
|
|---|
| [3646] | 111 | sa_size_t ncols = (imax-imin+1)/istep;
|
|---|
| 112 | sa_size_t nrows = 10;
|
|---|
| 113 | sa_size_t kc=0;
|
|---|
| 114 | for(int ifile=imin; ifile<=imax; ifile+=istep) {
|
|---|
| 115 | if (card==2)
|
|---|
| 116 | sprintf(fname, "%s/Ch34_%d.ppf",inoutpath.c_str(),ifile);
|
|---|
| 117 | else
|
|---|
| 118 | sprintf(fname, "%s/Ch12_%d.ppf",inoutpath.c_str(),ifile);
|
|---|
| 119 | cout << " ProcSVFiles[" << ifile << "] opening file " << fname << endl;
|
|---|
| 120 | PInPersist pin(fname);
|
|---|
| 121 | string tag1="specV1";
|
|---|
| 122 | string tag2="specV2";
|
|---|
| 123 | string tag12="visiV12";
|
|---|
| 124 | if (card==2) {
|
|---|
| 125 | tag1 = "specV3";
|
|---|
| 126 | tag2 = "specV4";
|
|---|
| 127 | tag12="visiV34";
|
|---|
| 128 | }
|
|---|
| 129 | TVector<r_4> sv1;
|
|---|
| 130 | TVector<r_4> sv2;
|
|---|
| 131 | TVector< complex<r_4> > vv12;
|
|---|
| 132 | pin >> PPFNameTag(tag1) >> sv1;
|
|---|
| 133 | pin >> PPFNameTag(tag2) >> sv2;
|
|---|
| 134 | pin >> PPFNameTag(tag12) >> vv12;
|
|---|
| 135 | if (ifile==imin) {
|
|---|
| 136 | nrows = sv1.Size();
|
|---|
| 137 | cout << " ProcSVFiles/Info: Output s1,s2 matrix size NRows=NFreq="
|
|---|
| 138 | << nrows << " NCols=NFiles=" << ncols << endl;
|
|---|
| 139 | s1.SetSize(nrows, ncols);
|
|---|
| 140 | s2.SetSize(nrows, ncols);
|
|---|
| 141 | nrows = vv12.Size();
|
|---|
| 142 | cout << " ProcSVFiles/Info: Output v12 matrix size NRows=NFreq="
|
|---|
| 143 | << nrows << " NCols=NFiles=" << ncols << endl;
|
|---|
| 144 | v12re.SetSize(nrows, ncols);
|
|---|
| 145 | v12im.SetSize(nrows, ncols);
|
|---|
| 146 | v12phi.SetSize(nrows, ncols);
|
|---|
| [3647] | 147 | v12mod.SetSize(nrows, ncols);
|
|---|
| [3646] | 148 | }
|
|---|
| 149 | s1.Column(kc) = sv1;
|
|---|
| 150 | s2.Column(kc) = sv2;
|
|---|
| 151 | v12re.Column(kc) = real(vv12);
|
|---|
| 152 | v12im.Column(kc) = imag(vv12);
|
|---|
| 153 | v12phi.Column(kc) = phase(vv12);
|
|---|
| [3647] | 154 | v12mod.Column(kc) = module(vv12);
|
|---|
| 155 |
|
|---|
| 156 | // Calcul moyenne dans des bandes en frequence
|
|---|
| 157 | int deltajf=(jf2-jf1)/nfreq;
|
|---|
| 158 | if (deltajf<1) deltajf=1;
|
|---|
| 159 | for(int kf=0; kf<nfreq; kf++) {
|
|---|
| 160 | sa_size_t jfstart=jf1+kf*deltajf;
|
|---|
| 161 | sa_size_t jfend=jfstart+deltajf;
|
|---|
| 162 | if (jfend>jf2) break;
|
|---|
| 163 | nmnt=0; ms1=ms2=ms12=ms12re=ms12im=ms12phi=ms12mod=0.;
|
|---|
| 164 | for(sa_size_t jf=jfstart; jf<jfend; jf++) {
|
|---|
| 165 | ms1 += s1(jf,kc);
|
|---|
| 166 | ms2 += s2(jf,kc);
|
|---|
| 167 | ms12re += v12re(jf,kc);
|
|---|
| 168 | ms12im += v12im(jf,kc);
|
|---|
| 169 | ms12phi += v12phi(jf,kc);
|
|---|
| 170 | ms12mod += v12mod(jf,kc);
|
|---|
| 171 | }
|
|---|
| 172 | nmnt = (jfend-jfstart);
|
|---|
| 173 | if (nmnt>0) {
|
|---|
| 174 | double fnorm = (double)nmnt;
|
|---|
| 175 | xnt[0] = ((int_8)(sv1.Info()["StartFC"])+(int_8)(sv1.Info()["EndFC"]))*0.5;
|
|---|
| 176 | xnt[1] = ((int_8)(sv1.Info()["StartTT"])+(int_8)(sv1.Info()["EndTT"]))*0.5;
|
|---|
| 177 | xnt[2] = kf;
|
|---|
| 178 | xnt[3] = ms1/fnorm;
|
|---|
| 179 | xnt[4] = ms2/fnorm;
|
|---|
| 180 | xnt[5] = ms12/fnorm;
|
|---|
| 181 | xnt[6] = ms12re/fnorm;
|
|---|
| 182 | xnt[7] = ms12im/fnorm;
|
|---|
| 183 | xnt[8] = ms12phi/fnorm;
|
|---|
| 184 | xnt[9] = ms12mod/fnorm;
|
|---|
| 185 | nt.Fill(xnt);
|
|---|
| 186 | }
|
|---|
| [3646] | 187 | }
|
|---|
| 188 | kc++;
|
|---|
| 189 |
|
|---|
| 190 | }
|
|---|
| 191 | if (card==2)
|
|---|
| 192 | sprintf(fname, "%s/Ch34mtx.ppf",inoutpath.c_str());
|
|---|
| 193 | else
|
|---|
| 194 | sprintf(fname, "%s/Ch12mtx.ppf",inoutpath.c_str());
|
|---|
| 195 |
|
|---|
| 196 | cout << nt;
|
|---|
| 197 | cout << "ProcSVFiles: Opening file " << fname << " for writing" << endl;
|
|---|
| 198 | POutPersist po(fname);
|
|---|
| 199 | string tag1="s1";
|
|---|
| 200 | string tag2="s2";
|
|---|
| 201 | string tag12r="v12re";
|
|---|
| 202 | string tag12i="v12im";
|
|---|
| 203 | string tag12p="v12phi";
|
|---|
| 204 | string tagnt="nt12";
|
|---|
| 205 | if (card==2) {
|
|---|
| 206 | tag1="s3";
|
|---|
| 207 | tag2="s4";
|
|---|
| 208 | tag12r="v34re";
|
|---|
| 209 | tag12i="v34im";
|
|---|
| 210 | tag12p="v34phi";
|
|---|
| 211 | tagnt="nt34";
|
|---|
| 212 | }
|
|---|
| 213 | po << PPFNameTag(tag1) << s1;
|
|---|
| 214 | po << PPFNameTag(tag2) << s2;
|
|---|
| 215 | po << PPFNameTag(tag12r) << v12re;
|
|---|
| 216 | po << PPFNameTag(tag12i) << v12im;
|
|---|
| 217 | po << PPFNameTag(tag12p) << v12phi;
|
|---|
| 218 | po << PPFNameTag(tagnt) << nt;
|
|---|
| 219 | cout << "ProcSVFiles: Matrices s1, s2, v12re, v12im, v12phi, NTuple nt written to file " << fname << endl;
|
|---|
| 220 | return 0;
|
|---|
| 221 | }
|
|---|
| 222 |
|
|---|
| 223 |
|
|---|
| 224 |
|
|---|
| 225 |
|
|---|
| 226 | /* --Fonction-- */
|
|---|
| 227 | int Usage(bool fgshort)
|
|---|
| 228 | {
|
|---|
| 229 | cout << " --- svv2mtx.cc : Read PPF files produced by mcrd to make matrices BAORadio" << endl;
|
|---|
| [3647] | 230 | cout << " Usage: mcrd InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq [card=1]" << endl;
|
|---|
| [3646] | 231 | if (fgshort) {
|
|---|
| 232 | cout << " mcrd -h for detailed instructions" << endl;
|
|---|
| 233 | return 1;
|
|---|
| 234 | }
|
|---|
| 235 | cout << " InOutPath : Input/Output directory name " << endl;
|
|---|
| 236 | cout << " Imin,Imax,IStep: Input PPF files sequence number \n"
|
|---|
| 237 | << " FileNames=InOutPath/Ch12_II.fits Imin<=II<=Imax II+=IStep \n"
|
|---|
| [3647] | 238 | << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for ntuple\n"
|
|---|
| [3646] | 239 | << " card=1 Ch12 , card=2 Ch34 " << endl;
|
|---|
| 240 | return 1;
|
|---|
| 241 | }
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
