Logo coherent WaveBurst  
Library Reference Guide
Logo
netpixel.cc
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Sergey Klimenko
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 
19 //---------------------------------------------------
20 // WAT pixel class for network analysis
21 // S. Klimenko, University of Florida
22 //---------------------------------------------------
23 
24 #define NETPIXEL_CC
25 #include <time.h>
26 #include <iostream>
27 #include <stdexcept>
28 #include "netpixel.hh"
29 
30 using namespace std;
31 
32 ClassImp(netpixel)
33 
34 // constructors
35 
37  theta = phi = -1.;
38  clusterID = time = frequency = 0;
39  rate = 1.; likelihood = 0.; core = false;
40  ellipticity = 0.; polarisation = 0.;
41  data.clear();
42  neighbors.clear();
43 }
44 
46  theta = phi = -1.;
47  clusterID = time = frequency = 0;
48  rate = 1.; likelihood = 0.; core = false;
49  ellipticity = 0.; polarisation = 0.;
50  pixdata pix;
51 
52  data.clear();
53  neighbors.clear();
54  if(!n) return;
55 
56  pix.noiserms = 0; // average noise rms
57  pix.wave = 0; // vector of pixel's wavelet amplitudes
58  pix.w_90 = 0; // vector of pixel's phase shifted wavelet amplitudes
59  pix.asnr = 0; // vector of pixel's whitened amplitudes
60  pix.a_90 = 0; // vector of pixel's phase shifted amplitudes
61  pix.rank = 0; // vector of pixel's rank amplitudes
62  pix.index = 0; // index in the wavearray for other detectors
63 
64  data.push_back(pix);
65  data.reserve(n);
66  for(size_t i=1; i<n; i++) data.push_back(pix);
67 }
68 
69 // netpixel::operator =
71 {
72  this->clusterID = value.clusterID; // cluster ID
73  this->time = value.time; // time index for zero detector (accounts for time shift)
74  this->frequency = value.frequency; // frequency index (layer)
75  this->layers = value.layers; // number of frequency layers
76  this->rate = value.rate; // wavelet layer rate
77  this->likelihood = value.likelihood; // likelihood
78  this->ellipticity = value.ellipticity; // likelihood
79  this->polarisation= value.polarisation;// likelihood
80  this->theta = value.theta; // source angle theta index
81  this->phi = value.phi; // source angle phi index
82  this->core = value.core; // pixel type: true - core , false - halo
83  this->data = value.data; // copy data
84  this->neighbors = value.neighbors; // copy neighbors
85  this->tdAmp = value.tdAmp; // copy TD vectors
86 
87  return *this;
88 }
89 
90 // write pixel in file
91 bool netpixel::write(const FILE *fp)
92 {
93  size_t i,j;
94  size_t n = this->neighbors.size();
95  size_t m = this->data.size();
96  size_t k = this->tdAmp.size();
97  size_t I = n ? n : 1;
98 
99  k = (k==m) ? this->tdAmp[0].size() : 0;
100 
101  double db[14];
102  wavearray<double> wb(I+m*(7+k));
103 
104  // write metadata
105 
106  db[0] = (double)m; // number of detectors
107  db[1] = (double)I; // number of neighbors
108  db[2] = (double)this->clusterID; // cluster ID
109  db[3] = this->time; // time index for zero detector (accounts for time shift)
110  db[4] = this->frequency; // frequency index (layer)
111  db[5] = this->core ? 1. : 0.; // pixel type: true - core , false - halo
112  db[6] = this->rate; // wavelet layer rate
113  db[7] = this->likelihood; // likelihood
114  db[8] = this->theta; // source angle theta index
115  db[9] = this->phi; // source angle phi index
116  db[10]= this->ellipticity; // waveform ellipticity
117  db[11]= this->polarisation; // waveform polarisation
118  db[12]= this->layers; // number of layers
119  db[13] = (double)k; // number of TD amplitudes per detector
120 
121  if(fwrite(db, 14*sizeof(double), 1, (FILE*)fp)!=1) return false;
122 
123  // write neighbors
124  for(i=0; i<I; i++) wb.data[i] = n ? neighbors[i] : 0.;
125 
126  // write amplitudes
127  for(i=0; i<m; i++) {
128  wb.data[I++] = this->data[i].noiserms; // average noise rms
129  wb.data[I++] = this->data[i].wave; // vector of pixel's wavelet amplitudes
130  wb.data[I++] = this->data[i].w_90; // vector of pixel's wavelet amplitudes
131  wb.data[I++] = this->data[i].asnr; // vector of pixel's whitened amplitudes
132  wb.data[I++] = this->data[i].a_90; // vector of pixel's phase shifted amplitudes
133  wb.data[I++] = (double)this->data[i].rank; // vector of pixel's rank amplitudes
134  wb.data[I++] = (double)this->data[i].index; // index in wavearray
135  }
136 
137  // write TD amplitudes
138  for(i=0; i<m; i++)
139  for(j=0; j<k; j++)
140  wb.data[I++] = (double)this->tdAmp[i].data[j];
141 
142  if(fwrite(wb.data, I*sizeof(double), 1, (FILE*)fp)!=1) return false;
143  return true;
144 }
145 
146 // read pixel from file
147 bool netpixel::read(const FILE *fp)
148 {
149  size_t i;
150  int j;
151  pixdata pix;
152  double db[14];
153 
154  this->clear();
155 
156  // read metadata
157  if(fread(db, 14*sizeof(double), 1, (FILE*)fp)!=1) return false;
158 
159  this->clusterID = size_t(db[2]+0.1); // cluster ID
160  this->time = size_t(db[3]+0.1); // time index for zero detector
161  this->frequency = size_t(db[4]+0.1); // frequency index (layer)
162  this->core = db[5]>0. ? true : false; // pixel type: true - core , false - halo
163  this->rate = (float)db[6]; // wavelet layer rate
164  this->likelihood = (float)db[7]; // likelihood
165  this->theta = (float)db[8]; // source angle theta index
166  this->phi = (float)db[9]; // source angle phi index
167  this->ellipticity = (float)db[10]; // waveform ellipticity
168  this->polarisation = (float)db[11]; // waveform polarisation
169  this->layers = size_t(db[12]+0.1); // number of layers
170 
171  size_t m = size_t(db[0]+0.1); // number of detectors
172  size_t I = size_t(db[1]+0.1); // neighbor size
173  size_t k = size_t(db[13]+0.1); // number TD amplitudes per detector
174  size_t n = I + m*(7+k); // buffer size
175  wavearray<double> wb(n);
176  wavearray<float> tda(k);
177 
178  if(fread(wb.data, n*sizeof(double), 1, (FILE*)fp)!=1) return false;
179 
180  // get neighbors
181  for(i=0; i<I; i++) {
182  j = wb.data[i]<0. ? int(wb.data[i]-0.1) : int(wb.data[i]+0.1);
183  if(j) this->neighbors.push_back(j);
184  }
185 
186  // get data
187  for(i=0; i<m; i++) {
188  pix.noiserms = wb.data[I++]; // average noise rms
189  pix.wave = wb.data[I++]; // vector of pixel's wavelet amplitudes
190  pix.w_90 = wb.data[I++]; // vector of pixel's wavelet amplitudes
191  pix.asnr = wb.data[I++]; // vector of pixel's whitened amplitudes
192  pix.a_90 = wb.data[I++]; // vector of pixel's phase shifted amplitudes
193  pix.rank = float(wb.data[I++]); // vector of pixel's rank amplitudes
194  pix.index = int(wb.data[I++]+0.1); // index in wavearray
195  this->data.push_back(pix);
196  }
197 
198  // get TD amplitudes
199  for(i=0; i<m; i++) {
200  for(j=0; j<(int)k; j++) {
201  tda.data[j] = float(wb.data[I++]);
202  }
203  if(k) this->tdAmp.push_back(tda);
204  }
205 
206  return true;
207 }
208 
209 
210 
211 
212 
213 
214 
215 
216 
217 
218 
219 
220 
221 
float phi
Definition: netpixel.hh:117
par [0] value
float rank
Definition: netpixel.hh:40
size_t clusterID
Definition: netpixel.hh:109
std::vector< wavearray< float > > tdAmp
Definition: netpixel.hh:123
size_t frequency
Definition: netpixel.hh:111
float likelihood
Definition: netpixel.hh:114
int n
Definition: cwb_net.C:28
std::vector< int > neighbors
Definition: netpixel.hh:124
double frequency
std::vector< pixdata > data
Definition: netpixel.hh:122
float theta
wavearray< double > rank
Definition: Regression_H1.C:80
netpixel pix(nifo)
int layers
STL namespace.
Long_t size
size_t layers
Definition: netpixel.hh:112
int m
Definition: cwb_net.C:28
int j
Definition: cwb_net.C:28
i drho i
bool write(const FILE *)
Definition: netpixel.cc:91
wc clear()
double asnr
Definition: netpixel.hh:38
bool core
Definition: netpixel.hh:120
double w_90
Definition: netpixel.hh:37
float phi
double wave
Definition: netpixel.hh:36
i() int(T_cor *100))
float polarisation
Definition: netpixel.hh:119
int index
Definition: netpixel.hh:41
void clear()
Definition: netpixel.hh:93
int k
size_t time
Definition: netpixel.hh:110
netpixel()
Definition: netpixel.cc:36
float ellipticity
Definition: netpixel.hh:118
double noiserms
Definition: netpixel.hh:35
wavearray< int > index
float theta
Definition: netpixel.hh:116
bool read(const FILE *)
Definition: netpixel.cc:147
DataType_t * data
Definition: wavearray.hh:319
float rate
Definition: netpixel.hh:113
double a_90
Definition: netpixel.hh:39
netpixel & operator=(const netpixel &)
Definition: netpixel.cc:70