Logo Coherent WaveBurst  
Reference Guide
Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
readfrfile.cc
Go to the documentation of this file.
1 /*---------------------------------------------------------------------
2  * Package: Wavelet Analysis Tool
3  * File name: readfrfile.cc
4  * Authors: S.Klimenko, A.Sazonov : University of Florida, Gainesville
5  *
6  * Virgo Frame Library rev. >= 4.41 is required.
7  *---------------------------------------------------------------------
8 */
9 
10 #include "readfrfile.hh"
11 
12 // Function ReadFrFile transparently reads both WAT-compressed
13 // and standard-compressed ADC data from single frame file.
14 // wavearray<double> - output array,
15 // tlen - length of data in seconds,
16 // tskip - skip time from the begining of the file, in seconds
17 // *cname - ADC channel name
18 // *fname - frame file name
19 // seek - if true then it is allowed to seek data continuation
20 // in next file which name the function will try to guess
22  double tlen,
23  double tskip,
24  char *cname,
25  char *fname,
26  bool seek,
27  char *channel_type)
28 {
29  FrFile *iFile = FrFileINew(fname);
30  if(iFile == NULL)
31  {
32  printf(" ReadFrFile(): cannot open the input file %s\n", fname);
33  return false;
34  }
35  if(&out == NULL)
36  {
37  printf(" ReadFrFile(): no output array specified \n");
38  return false;
39  }
40 
41  bool failure = false; // true if channel not found
42  int n=0, l=0, k=0;
43  int nu, nn=0;
44  double rate=0., d;
45 
46  short ref = 0x1234; // for byte-swapping check
47  char *refC = (char *)&ref;
48 
49  int iskip = 0; // skip in index space
50  if (tskip<0) tskip=0.;
51 
53 
54  char nname[512]; // char array for computable file name
55  sprintf(nname,"%s",fname);
56 
57 // printf("%s %s\n",nname,fname);
58 
59  char gpstime[12];
60  char newgpstime[12];
61  char *ptime; // pointer to GPS time substring in file name
62 
63  FrAdcData *adc = NULL;
64  FrProcData *proc = NULL;
65  FrSerData *ser = NULL;
66  FrSimData *sim = NULL;
67  FrVect *datap = NULL;
68 
69 // frame length (should be determined from frame file )
70  double frlen = 16.;
71 
72 // this doesn't work
73 // frlen = iFile->length; // frame length
74 // cout<<frlen<<endl;
75 
76  double gt,gt0,gts,gte; // GPS time, GPS time of file end
77  gts = FrFileITStart(iFile);
78  gte = FrFileITEnd(iFile);
79 
80 // nskip is number of full frames (or multiframe files) in tskip
81  int nskip = int( tskip/frlen );
82 
83 // fskip is a fraction of tskip: 0 <= fskip < frlen
84  double fskip = tskip - nskip*frlen;
85 
86  gt = gts + nskip*frlen; // begin GPS time of the first frame to read
87  gt0 = gts + tskip; // GPS time where data begins
88 
89 //printf(" gts = %16.6f, gte = %16.6f, gt = %16.6f\n", gts, gte, gt);
90 
91  iFile->compress = 255;
92 
93  do { // read many files if names are computable
94 
95  while(gt<gte) { // read frames from a file
96 
97 //printf(" gts = %16.6f, gte = %16.6f, gt = %16.6f\n", gts, gte, gt);
98 
99 // gt+1.e-7 in function call to prevent uncertainty in frame access
100 
101  failure = false;
102 
103  if (!strcasecmp(channel_type,"adc"))
104  {
105  adc = FrAdcDataReadT(iFile,cname,gt+1.e-7);
106  if (adc == NULL) { failure = true; }
107  else { datap = adc->data; }
108  }
109  else if (!strcasecmp(channel_type,"proc"))
110  {
111  proc = FrProcDataReadT(iFile,cname,gt+1.e-7);
112  if (proc == NULL) { failure = true; }
113  else { datap = proc->data; }
114  }
115  else
116  {
117  printf(" ReadFrFile() error: unknown channel type %s\n",
118  channel_type);
119  }
120 
121  if (failure) {
122  printf(" ReadFrFile() error: channel %s is not found in file %s\n",
123  cname, nname);
124  return false;
125  }
126 
127  gt += frlen;
128 
129  if(datap == NULL ) continue;
130 
131 // printf("startx=%f%s\n",datap->startX[0],datap->unitX[0]);
132 
133  n = datap->nData;
134  nu = n;
135  if ((int)wd.size() != n) wd.resize(n);
136 
137  if (!nn) {
138 // rate = adc->sampleRate; // not available for proc
139  if (!datap->dx[0]==0.) rate = 1./datap->dx[0];
140  nn = int(tlen * rate);
141 
142  if ( nn == 0 ) //
143  {
144  printf(" ReadFrFile: time interval too short, return no data.\n");
145  return false;
146  }
147 
148  if(out.size() != size_t(nn)) out.resize(nn);
149 
150 // k counts aready done data samples
151 // l counts how much to do in current frame
152 // set initial values
153  k = 0;
154  l = nn;
155  }
156 
157  iskip = int(fskip * rate);
158  l = l < (n - iskip)? l : n - iskip;
159 
160 // do we need to swap the bytes in compressed data?
161  bool swap = (refC[0] == 0x34 && datap->compress == 255) ||
162  (refC[0] == 0x12 && datap->compress == 511);
163 
164 // Swap bytes assuming the vector is used as an int vector
165 // Note: this code perform very slowly on Alpha CPU unless
166 // compiler can take advantage of instruction sets of ev56
167 // or ev6 architecture (option -mcpu=ev56 or ev6 for gcc/g++)
168  if(swap) {
169  unsigned char local[2];
170  char *buf = datap->data;
171 
172  for(unsigned int i=0; i<datap->nBytes-3; i=i+4) {
173  local[0] = buf[3];
174  local[1] = buf[2];
175  buf[3] = buf[0];
176  buf[2] = buf[1];
177  buf[1] = local[1];
178  buf[0] = local[0];
179  buf = buf + 4;
180  }
181  }
182 
183  if( (datap->compress & 0xff) == 255 ){
184 
185  nu = unCompress(datap->dataI, wd);
186 
187  if (nu!=n) {
188  printf(" ReadFrFile: unCompress returned wrong data length\n");
189  break; // break the frame reading loop
190  }
191 
192 // round data for compatibility with uncompressed data stored in frame files
193  switch(datap->type) {
194 
195  case FR_VECT_2S:
196  for(int i=iskip; i<n; i++) {
197  d=wd.data[i];
198  wd.data[i]=float(short(d>0.? d+0.5 : d-0.5));
199  }
200  break;
201 
202  default:
203  break;
204  }
205  for(int i=iskip; i<l; i++)
206  out.data[k+i]=double(wd.data[i]);
207  }
208 
209  else {
210 
211  switch(datap->type) {
212 
213  case FR_VECT_2S:
214  for(int i=0; i<l; i++)
215  out.data[i+k] = double(datap->dataS[i+iskip]);
216  break;
217 
218  case FR_VECT_4R:
219  for(int i=0; i<l; i++)
220  out.data[i+k] = double(datap->dataF[i+iskip]);
221  break;
222 
223  default:;
224  }
225  }
226 
227  if (adc) FrAdcDataFree(adc);
228  if (proc) FrProcDataFree(proc);
229  if (sim) FrSimDataFree(sim);
230  if (ser) FrSerDataFree(ser);
231 
232  k += l;
233  l = nn - k; // how much to read?
234  fskip = 0.; // fskip applyes to first frame only
235 
236  if ( l <= 0 ) break; // break the frame reading loop
237 
238  } // end of frames reading loop
239 
240  FrFileIEnd(iFile);
241 
242  if (nn && ((nn-k)<=0)) break; // no more data to read
243 
244 // try to calculate next file name
245  sprintf(gpstime, "%9d", int(gts));
246  sprintf(newgpstime, "%9d", int(gte));
247  ptime=strstr(nname, gpstime);
248 
249  if ( ptime != NULL && atoi(ptime) == int(gts) &&
250  strlen(gpstime) == strlen(newgpstime) ){
251 
252  strncpy(ptime, newgpstime, strlen( newgpstime ));
253 // printf(" guess next file name to be %s\n",nname);
254 
255  }
256  else break; // break file search loop
257 
258  iFile = FrFileINew(nname);
259 
260  if(iFile == NULL) {
261  printf(" ReadFrFile(): cannot open next input file %s\n", nname);
262  break; // break file search loop
263  }
264 
265  gts = FrFileITStart(iFile);
266  iFile->compress = 255;
267 
268  if (gts != gte) {
269  printf(" ReadFrFile(): next input file");
270  printf(" %s doesn't provide continuous data\n", nname);
271  FrFileIEnd(iFile);
272  break; // break file search loop
273  }
274 
275  gte = FrFileITEnd(iFile);
276 
277  } while (seek); // end of file search loop
278 
279  if (out.size() != size_t(nn)) return false;
280 
281  if (( nn - k ) > 0) // fill rest of data by zeroes
282  {
283  for (int i=k; i<nn; i++) out.data[i]=0.;
284  }
285 
286 // delete nname;
287 // if (out == 0) return NULL;
288 
289 // normal return
290  out.rate(rate);
291  out.start(gt0);
292 
293  return true;
294 }
295 
296 
297 // Function ReadFrFile transparently reads both WAT-compressed
298 // and standard-compressed ADC data from single frame file.
299 // tlen - length of data in seconds,
300 // tskip - skip time from the begining of the file, in seconds
301 // *cname - ADC channel name
302 // *fname - frame file name
303 // seek - if true then it is allowed to seek data continuation
304 // in next file which name the function will try to guess
306  double tskip,
307  char *cname,
308  char *fname,
309  bool seek,
310  char *channel_type)
311 {
312  FrFile *iFile = FrFileINew(fname);
313  if(iFile == NULL)
314  {
315  printf(" ReadFrFile(): cannot open the input file %s\n", fname);
316  return NULL;
317  }
318 
319  bool failure = false; // true if channel not found
320  int n=0, l=0, k=0;
321  int nu, nn=0;
322  double rate=0., d;
323 
324  short ref = 0x1234; // for byte-swapping check
325  char *refC = (char *)&ref;
326 
327  int iskip = 0; // skip in index space
328  if (tskip<0) tskip=0.;
329 
330  wavearray<float> *out = NULL;
331  wavearray<float> wd;
332 
333  char* nname = new char[strlen(fname)]; // char array for computable file name
334  strcpy(nname,fname); // copy original file name
335 
336  char gpstime[12];
337  char newgpstime[12];
338  char *ptime; // pointer to GPS time substring in file name
339 
340  FrAdcData *adc = NULL;
341  FrProcData *proc = NULL;
342  FrSerData *ser = NULL;
343  FrSimData *sim = NULL;
344  FrVect *datap = NULL;
345 
346 // frame length (should be determined from frame file )
347  double frlen = 16.;
348 
349 // this doesn't work
350 // frlen = iFile->length; // frame length
351 // cout<<frlen<<endl;
352 
353  double gt,gt0,gts,gte; // GPS time, GPS time of file end
354  gts = FrFileITStart(iFile);
355  gte = FrFileITEnd(iFile);
356 
357 // nskip is number of full frames (or multiframe files) in tskip
358  int nskip = int( tskip/frlen );
359 
360 // fskip is a fraction of tskip: 0 <= fskip < frlen
361  double fskip = tskip - nskip*frlen;
362 
363  gt = gts + nskip*frlen; // begin GPS time of the first frame to read
364  gt0 = gts + tskip; // GPS time of data begin
365 
366 //printf(" gts = %16.6f, gte = %16.6f, gt = %16.6f\n", gts, gte, gt);
367 
368  iFile->compress = 255;
369 
370  do { // read many files if names are computable
371 
372  while(gt<gte) { // read frames from file
373 
374 //printf(" gts = %16.6f, gte = %16.6f, gt = %16.6f\n", gts, gte, gt);
375 
376 // gt+1.e-7 in function call to prevent uncertainty in frame access
377 
378  failure = false;
379 
380  if (!strcasecmp(channel_type,"adc"))
381  {
382  adc = FrAdcDataReadT(iFile,cname,gt+1.e-7);
383  if (adc == NULL) { failure = true; }
384  else { datap = adc->data; }
385  }
386  else if (!strcasecmp(channel_type,"proc"))
387  {
388  proc = FrProcDataReadT(iFile,cname,gt+1.e-7);
389  if (proc == NULL) { failure = true; }
390  else { datap = proc->data; }
391  }
392  else
393  {
394  printf(" ReadFrFile() error: unknown channel type %s\n",
395  channel_type);
396  }
397 
398  if (failure) {
399  printf(" ReadFrFile() error: channel %s is not found in file %s\n",
400  cname, nname);
401  if (out) delete out;
402  return NULL;
403  }
404 
405  gt += frlen;
406 
407  if(datap == NULL ) continue;
408 
409 // printf("startx=%f%s\n",datap->startX[0],datap->unitX[0]);
410  n = datap->nData;
411  nu = n;
412  if ((int)wd.size() != n) wd.resize(n);
413 
414  if (!out) {
415 // rate = adc->sampleRate; // not available for proc
416  if (!datap->dx[0]==0.) rate = 1./datap->dx[0];
417  nn = int(tlen * rate);
418 
419  if ( nn == 0 ) //
420  {
421  printf(" ReadFrFile: time interval too short, return no data.\n");
422  return NULL;
423  }
424 
425  out = new wavearray<float>(nn);
426 
427 // k counts aready done data samples
428 // l counts how much to do in current frame
429 // set initial values
430  k = 0;
431  l = nn;
432  }
433 
434  iskip = int(fskip * rate);
435  l = l < (n - iskip)? l : n - iskip;
436 
437 // do we need to swap the bytes in compressed data?
438  bool swap = (refC[0] == 0x34 && datap->compress == 255) ||
439  (refC[0] == 0x12 && datap->compress == 511);
440 
441 // Swap bytes assuming the vector is used as an int vector
442 // Note: this code perform very slowly on Alpha CPU unless
443 // compiler can take advantage of instruction sets of ev56
444 // or ev6 architecture (option -mcpu=ev56 or ev6 for gcc/g++)
445  if(swap) {
446  unsigned char local[2];
447  char *buf = datap->data;
448 
449  for(unsigned int i=0; i<datap->nBytes-3; i=i+4) {
450  local[0] = buf[3];
451  local[1] = buf[2];
452  buf[3] = buf[0];
453  buf[2] = buf[1];
454  buf[1] = local[1];
455  buf[0] = local[0];
456  buf = buf + 4;
457  }
458  }
459 
460  if( (datap->compress & 0xff) == 255 ){
461 
462  nu = unCompress(datap->dataI, wd);
463 
464  if (nu!=n) {
465  printf(" ReadFrFile: unCompress returned wrong data length\n");
466  break; // break the frame reading loop
467  }
468 
469 // round data for compatibility with uncompressed data stored in frame files
470  switch(datap->type) {
471 
472  case FR_VECT_2S:
473  for(int i=iskip; i<n; i++) {
474  d=wd.data[i];
475  wd.data[i]=float(short(d>0.? d+0.5 : d-0.5));
476  }
477  break;
478 
479  default:
480  break;
481  }
482  out->cpf(wd,l,iskip,k);
483  }
484 
485  else {
486 
487  switch(datap->type) {
488 
489  case FR_VECT_2S:
490  for(int i=0; i<l; i++)
491  out->data[i+k] = datap->dataS[i+iskip];
492  break;
493 
494  case FR_VECT_4R:
495  for(int i=0; i<l; i++)
496  out->data[i+k] = datap->dataF[i+iskip];
497  break;
498 
499  default:;
500  }
501  }
502 
503  if (adc) FrAdcDataFree(adc);
504  if (proc) FrProcDataFree(proc);
505  if (sim) FrSimDataFree(sim);
506  if (ser) FrSerDataFree(ser);
507 
508  k += l;
509  l = nn - k; // how much to read?
510  fskip = 0.; // fskip applyes to first frame only
511 
512  if ( l <= 0 ) break; // break the frame reading loop
513  } // end of frames reading loop
514 
515  FrFileIEnd(iFile);
516 
517  if (out && ((nn-k)<=0)) break; // no more data to read
518 
519 // try to calculate next file name
520  sprintf(gpstime, "%9d", int(gts));
521  sprintf(newgpstime, "%9d", int(gte));
522  ptime=strstr(nname, gpstime);
523 
524  if ( ptime != NULL && atoi(ptime) == int(gts) &&
525  strlen(gpstime) == strlen(newgpstime) ){
526 
527  strncpy(ptime, newgpstime, strlen( newgpstime ));
528 // printf(" guess next file name to be %s\n",nname);
529 
530  }
531  else break; // break file search loop
532 
533  iFile = FrFileINew(nname);
534 
535  if(iFile == NULL) {
536  printf(" ReadFrFile(): cannot open next input file %s\n", nname);
537  break; // break file search loop
538  }
539 
540  gts = FrFileITStart(iFile);
541  iFile->compress = 255;
542 
543  if (gts != gte) {
544  printf(" ReadFrFile(): next input file");
545  printf(" %s doesn't provide continuous data\n", nname);
546  FrFileIEnd(iFile);
547  break; // break file search loop
548  }
549 
550  gte = FrFileITEnd(iFile);
551 
552  } while (seek); // end of file search loop
553 
554  if (out==NULL) return NULL;
555 
556  if (( nn - k ) > 0) // fill rest of data by zeroes
557  {
558  for (int i=k; i<nn; i++) out->data[i]=0.;
559  }
560 
561 // delete nname;
562 
563  if (out == 0) return NULL;
564 
565 // normal return
566  out->rate(rate);
567  out->start(gt0);
568 
569  return out;
570 }
571 
virtual size_t size() const
Definition: wavearray.hh:127
printf("total live time: non-zero lags = %10.1f \n", liveTot)
virtual void rate(double r)
Definition: wavearray.hh:123
int n
Definition: cwb_net.C:10
int unCompress(int *in, wavearray< float > &out)
Definition: lossy.cc:158
virtual void start(double s)
Definition: wavearray.hh:119
i drho i
TChain sim("waveburst")
ofstream out
Definition: cwb_merge.C:196
i() int(T_cor *100))
char fname[1024]
int k
double e
bool ReadFrFile(wavearray< double > &out, double tlen, double tskip, char *cname, char *fname, bool seek, char *channel_type)
Definition: readfrfile.cc:21
strcpy(RunLabel, RUN_LABEL)
int l
Definition: cbc_plots.C:434
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
DataType_t * data
Definition: wavearray.hh:301
virtual void resize(unsigned int)
Definition: wavearray.cc:445
void cpf(const wavearray< DataType_t > &, int=0, int=0, int=0)
Definition: wavearray.cc:699