Logo coherent WaveBurst  
Library Reference Guide
Logo
detector.hh
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Sergey Klimenko, Gabriele Vedovato
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 // Wavelet Analysis Tool
21 // Sergey Klimenko, University of Florida
22 // universal data container for x-correlation analysis
23 // used with DMT and ROOT
24 //**********************************************************
25 
26 #ifndef DETECTOR_HH
27 #define DETECTOR_HH
28 
29 #include <iostream>
30 #include <vector>
31 #include "complex"
32 #include "wseries.hh"
33 #include "netcluster.hh"
34 #include "skymap.hh"
35 #include "wavecomplex.hh"
36 #include "skycoord.hh"
37 #include "TString.h"
38 #include "TNamed.h"
39 #include "sseries.hh"
40 
42 typedef std::vector<int> vector_int;
43 
44 struct delayFilter {
45  std::vector<short> index; // relative wavelet array index
46  std::vector<float> value; // amplitude
47 };
48 
50  char name[32];
51  double latitude;
52  double longitude;
53  double elevation;
54  double AltX; // elevation of the x arm
55  double AzX; // azimut of the x arm (angle-deg from nord)
56  double AltY; // elevation of the y arm
57  double AzY; // azimut of the y arm (angle-deg from nord)
58 };
59 
61  TENSOR = 0,
62  SCALAR = 1
63 };
64 
66 
67 class detector : public TNamed
68 {
69  public:
70 
71  // constructors
72 
73  //: Default constructor
74  detector();
75 
76  //: detector name constructor
77  detector(char*, double=0.);
78 
79  //: user defined detector constructor
80  detector(detectorParams, double=0.);
81 
82  //: Copy constructor
83  //!param: value - object to copy from
84  detector(const detector&); // use with caution
85 
86  //: destructor
87  virtual ~detector();
88 
89  // operators
90 
91  detector& operator= (const detector&); // use with caution
92  detector& operator= (const WSeries<double> &); // use with caution
93  detector& operator<< (detector&); // copy 'from' inj/rec stuff
94  detector& operator>> (detector&); // copy 'to' inj/rec stuff
95 
96  // accessors
97 
98  //: return antenna pattern
99  //!param: source theta,phi, polarization angle psi in degrees
100  wavecomplex antenna(double, double, double=0.);
101 
102  // returns the detector parameters
103  detectorParams getDetectorParams();
104 
105  //: return detector delay
106  // time delay convention: t_detector-tau - arrival time at the center of Earth
107  //!param: source theta,phi angles in degrees
108  double getTau(double, double);
109 
110  //: set detector radius vector
111  //!param: Rx,Ry,Rz in ECEF frame
112 // void setRV(double, double, double);
113  //: set detector x-arm unit vector
114  //!param: Ex,Ey,Ez in ECEF frame
115 // void setEx(double, double, double);
116  //: set detector y-arm unit vector
117  //!param: Ex,Ey,Ez in ECEF frame
118 // void setEy(double, double, double);
119 
120  //: initialize detector tensor
121  //!param: no parameters
122  void init();
123 
124  //: initialize delay filter for non-heterodine wavelet
125  //!param: filter length
126  //!param: filter phase delay in degrees
127  //!param: number of up-sample layers
128  //!return filter size
129  size_t setFilter(size_t, double=0., size_t=0);
130 
131  //: initialize delay filter from another detector
132  //!param: detector
133  //!return filter size
134  size_t setFilter(detector&);
135 
136  //: write/read delay filter to file in the format:
137  //: first string: wavelet filter length, number of wavelet layers
138  //!param: file name
139  void writeFilter(const char*);
140  void readFilter(const char*);
141 
142  // clear filter and release memory (need for vector class)
143  inline void clearFilter() {
144  filter.clear();
145  std::vector<delayFilter>().swap(filter); // release memory
146  }
147 
148  //: apply sample delay to timeseries stored in TFmap to adjast
149  // timing for a given theta and phy sky location.
150  //!param: theta
151  //!param: phi
152  void delay(double, double);
153 
154  //: apply sample delay to input timeseries to adjast
155  // timing for a given theta and phy sky location. Do not store in the TFmap
156  //!param: input TS
157  //!param: theta
158  //!param: phi
159  void delay(wavearray<double> &, double, double);
160 
161  //: apply delay T to input timeseries.
162  //!param: input TS
163  //!param: time delay T
164  void delay(wavearray<double> &, double);
165 
166  //: apply non-heterodine delay filter to input WSeries and put result in TFmap
167  // time delay convention: + - shift TS right
168  // - - shift TS left
169  //!param: time delay in seconds
170  //!param: WSeries which should be delayed
171  void delay(double, WSeries<double> &);
172 
173  //: get pointer to time series data
174  //!param: no parameters
175  inline wavearray<double>* getHoT() { return &HoT; }
176 
177  //: get pointer to TF data
178  //!param: no parameters
179  inline WSeries<double>* getTFmap() { return &TFmap; }
180 
181  //: get sparse map index in vSS vector for input resolution (wavelet rate r)
182  inline int getSTFind(double r) {
183  for(size_t i=0; i<vSS.size(); i++) if(int(vSS[i].wrate()-r)==0) return i;
184  return vSS.size()+1;
185  }
186  //: get pointer to sparse TF data for input index n
187  inline SSeries<double>* getSTFmap(size_t n) { return n<vSS.size() ? &vSS[n] : NULL; }
188 
189  // operations with sparse maps
190  inline size_t ssize() { return vSS.size(); } // get size of sparse TF vector (number of resolutions)
191  inline void sclear() { vSS.clear(); } // clear
192  inline void addSTFmap(netcluster* pwc, double mTau=0.042);
193 
194  //: get size of sparse TF vector (number of resolutions)
195  //!param: no parameters
196 
197  //: reconstruct wavelet series for a cluster, put it in waveForm
198  //! param: input parameter is the cluster ID in the netcluster structure
199  //! param: input netcluster structure
200  //! param: amplitude type
201  //! param: amplitude index
202  //! return: cluster average noise rms
203  double getwave(int, netcluster&, char, size_t);
204 
205  //: set tau array
206  // time delay convention: t_detector-tau - arrival time at the center of Earth
207  //!param - step on phi and theta
208  //!param - theta begin
209  //!param - theta end
210  //!param - phi begin
211  //!param - phi end
212  void setTau(double,double=0.,double=180.,double=0.,double=360.);
213 
214  //: set tau array
215  // time delay convention: t_detector-tau - arrival time at the center of Earth
216  //!param - healpix order
217  void setTau(int);
218 
219  //: set antenna patterns
220  //!param - step on phi and theta
221  //!param - theta begin
222  //!param - theta end
223  //!param - phi begin
224  //!param - phi end
225  void setFpFx(double,double=0.,double=180.,double=0.,double=360.);
226 
227  //: set antenna patterns
228  //!param - healpix order
229  void setFpFx(int);
230 
231  //: return noise variance for selected wavelet layer or pixel
232  //: used in coherence() and Rank() functions
233  //!param: wavelet layer index (frequency)
234  //!param: wavelet time index
235  //!if param 2 is specified - return noise rms for specified pixel location
236  //!if param 2 is not specified - return rms averaged over the layer
237  double getNoise(size_t,int=-1);
238 
239  /* calculate and save average noise rms for each pixel in netcluster wc
240  * @memo save pixel noise rms
241  * @param input netcluster
242  * @param detector index in tne netcluster
243  */
244  bool setrms(netcluster*, size_t=0);
245 
246  //: return pointer to noise array
247 // WSeries<double>* getNoise();
248 
249  //: whiten data in TFmap and save noise RMS in nRMS
250  // param 1 - time window dT. if = 0 - dT=T, where T is wavearray duration
251  // param 2 - 0 - no whitening, 1 - single whitening, >1 - double whitening
252  // param 3 - boundary offset
253  // param 4 - noise sampling interval (window stride)
254  // the number of measurements is k=int((T-2*offset)/stride)
255  // if stride=0, then stride is set to dT
256  // output: save array of noise RMS in nRMS
257  //!what it does: see algorithm description in wseries.hh
258  void white(double dT=0, int wtype=1,double offset=0.,double stride=0.) {
259  nRMS = TFmap.white(dT, wtype, offset, stride); return;
260  }
261 
262  // set constant time shift
263  inline void shift(double s) {
264  int I = TFmap.maxLayer()+1;
265  double R = TFmap.wavearray<double>::rate();
266  if(TFmap.size()<2 || R<=0.) this->sHIFt = s;
267  else this->sHIFt = int(s*R/I+0.1)*I/R;
268  }
269  // get constant time shift
270  inline double shift() { return this->sHIFt; }
271 
272  // rotate arms in the detector plane by angle a in degrees
273  void rotate(double);
274 
275  // apply band pass filter with cut-offs specified by parameters
276  // f1>=0 && f2>=0 : band pass
277  // f1<0 && f2<0 : band cut
278  // f1<0 && f2>=0 : low pass
279  // f1>=0 && f2<0 : high pass
280  // f1==0 : f1 = TFmap.getlow()
281  // f2==0 : f2 = TFmap.gethigh()
282  void bandPass1G(double f1=0.,double f2=0.); // used by 1G
283  void bandPass(double f1,double f2, double a=0.){
284  this->TFmap.bandpass(f1,f2,a);
285  };
286  void bandCut(double f1,double f2) {bandPass(-fabs(f1),-fabs(f2));}
287  void lowPass(double f1) {bandPass(-fabs(f1),0);}
288  void highPass(double f2) {bandPass(0,-fabs(f2));}
289 
290  // calculate hrss of injected responses
291  // param 1 - injection wavelet series
292  // param 2 - array with injection time
293  // param 3 - integration window in seconds
294  // param 4 - save input waveform
295  size_t setsim(WSeries<double> &, std::vector<double>*, double=5., double=8., bool saveWF=false);
296 
297  // modify input signals (wi) at times pT according the factor pF
298  size_t setsnr(wavearray<double> &, std::vector<double>*, std::vector<double>*, double=5., double=8.);
299 
300  // print detector parameters
301  void print(); // *MENU*
302  virtual void Browse(TBrowser *b) {print();}
303 
304  bool isBuiltin() {return TString(dP.name).Sizeof()>1 ? false : true;}
305 
308 
309  inline double get_SS(){double rms=waveForm.rms(); return rms*rms*waveForm.size();}
310  inline double get_XX(){double rms=waveBand.rms(); return rms*rms*waveBand.size();}
311  inline double get_NN(){double rms=waveNull.rms(); return rms*rms*waveNull.size();}
312  inline double get_XS(){double rmx=waveBand.rms(); double rms=waveForm.rms(); return rmx*rms*waveBand.size();}
313  double getWFfreq(char atype='S');
314  double getWFtime(char atype='S');
315 
316  int wfsave() {return this->wfSAVE;}
317  void wfsave(int wfSAVE) {this->wfSAVE = (wfSAVE>=0)&&(wfSAVE<=3) ? wfSAVE : 0;}
318  // 0 : save detector definitions
319  // 1 : save injected waveforms
320  // 2 : save reconstructed waveforms
321  // 3 : save injected & reconstructed waveforms
322 
323 // data members
324 
325 //: position in ECEF frame
326 
327  char Name[16]; // detector name
328  size_t ifoID; // detector ID in the network - set up by network method
329  detectorParams dP; // user detector parameters
330  double Rv[3]; // radius vector to beam splitter
331  double Ex[3]; // vector along x-arm
332  double Ey[3]; // vector along y-arm
333  double DT[9]; // detector tenzor
334  double ED[5]; // network energy disbalance
335  double sHIFt; // time shifts for background analysis
336  double null; // unbiased null stream
337  double enrg; // total energy of PC components
338  double sSNR; // reconstructed response s-SNR
339  double xSNR; // reconstructed response x-SNR
340  double ekXk; // mean of reconstructed detector response
341  double rate; // original data rate (before downsampling)
342  size_t nDFS; // number of Delay Filter Samples
343  size_t nDFL; // number of Delay Filter Layers
344  int wfSAVE; // used in streamer method to save waveforms stuff
345 
346  skymap tau; // detector delay with respect to ECEF
347  skymap mFp; // F+ skymap
348  skymap mFx; // Fx skymap
349 
350  wavearray<double> HoT; // detector time series
351 
352  std::vector<SSeries<double> > vSS; // sparse TFmap
353 
354  WSeries<double> TFmap; // wavelet data
355  WSeries<double> waveForm; // buffer for a waveform
356  WSeries<double> waveBand; // buffer for a bandlimited waveform
357  WSeries<double> waveNull; // buffer for noise = data - signal
358  WSeries<double> nRMS; // noise RMS
359  WSeries<float> nVAR; // noise variability
360 
361  std::vector<delayFilter> filter; // delay filter
362 
363  wavearray<double> fp; // sorted F+ pattern
364  wavearray<double> fx; // sorted Fx pattern
365  wavearray<double> ffp; // sorted F+ * F+ + Fx * Fx
366  wavearray<double> ffm; // sorted F+ * F+ - Fx * Fx
367  wavearray<double> fpx; // sorted F+ * Fx * 2
368  wavearray<short> index; // index array for delayed amplitude (network)
369  wavearray<double> lagShift; // time shifts for background analysis
370 
371  wavearray<double> HRSS; // hrss of injected signals
372  wavearray<double> ISNR; // injected SNR
373  wavearray<double> FREQ; // frequency of injected signals
374  wavearray<double> BAND; // bandwith of injected signals
375  wavearray<double> TIME; // central time of injected signals
376  wavearray<double> TDUR; // duration of injected signals
377 
378  std::vector<int> IWFID; // injected waveforms ID
379  std::vector<wavearray<double>*> IWFP; // injected waveforms pointers
380 
381  std::vector<int> RWFID; // reconstructed waveforms ID
382  std::vector<wavearray<double>*> RWFP; // reconstructed waveforms pointers
383 
384  POLARIZATION polarization; // gw polarization states : TENSOR (fp,fx) SCALAR (fo)
385 
386  ClassDef(detector,4)
387 
388 }; // class detector
389 
390 
391 inline void detector::addSTFmap(netcluster* pwc, double mTau){
392  SSeries<double> SS; vSS.push_back(SS);
393  int j = vSS.size()-1;
394  vSS[j].SetMap(&TFmap);
395  vSS[j].SetHalo(mTau);
396  vSS[j].AddCore(ifoID,pwc);
397  vSS[j].UpdateSparseTable();
398 }
399 
400 inline size_t minDFindex(delayFilter &F){
401  size_t n = F.value.size();
402  if(!n) return 0;
403  size_t j = 0;
404  double x = 1.e13;
405  for(size_t i=0; i<n; i++){
406  if(x>fabs(F.value[i])) { x = fabs(F.value[i]); j=i; }
407  }
408  return j+1;
409 }
410 
411 #endif // DETECTOR_HH
412 
detectorParams dP
Definition: detector.hh:329
double sHIFt
Definition: detector.hh:335
int wfSAVE
Definition: detector.hh:344
void white(double dT=0, int wtype=1, double offset=0., double stride=0.)
what it does: see algorithm description in wseries.hh
Definition: detector.hh:258
void highPass(double f2)
Definition: detector.hh:288
wavearray< double > HoT
Definition: detector.hh:350
int offset
Definition: TestSTFT_2.C:19
double sSNR
Definition: detector.hh:338
wavearray< double > a(hp.size())
double AzX
Definition: detector.hh:55
par [0] name
int n
Definition: cwb_net.C:28
void wfsave(int wfSAVE)
Definition: detector.hh:317
void shift(double s)
Definition: detector.hh:263
TString("c")
wavearray< double > fx
Definition: detector.hh:364
WSeries< double > WSeriesD
Definition: detector.hh:65
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1) *itime+ifreq;double time=itime *dT;double freq=(ifreq >0) ? ifreq *dF :dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
std::vector< delayFilter > filter
Definition: detector.hh:361
std::vector< wavearray< double > * > RWFP
Definition: detector.hh:382
wavearray< double > HRSS
Definition: detector.hh:371
size_t nDFS
Definition: detector.hh:342
netcluster * pwc
Definition: cwb_job_obj.C:38
void lowPass(double f1)
Definition: detector.hh:287
double get_SS()
Definition: detector.hh:309
wavearray< short > index
Definition: detector.hh:368
std::vector< SSeries< double > > vSS[NIFO_MAX]
double AltX
Definition: detector.hh:54
double AzY
Definition: detector.hh:57
int polarization
void setPolarization(POLARIZATION polarization=TENSOR)
Definition: detector.hh:306
WSeries< double > waveBand
Definition: detector.hh:356
std::vector< int > vector_int
Definition: detector.hh:42
int j
Definition: cwb_net.C:28
size_t minDFindex(delayFilter &F)
Definition: detector.hh:400
double rate
Definition: detector.hh:341
i drho i
double longitude
Definition: detector.hh:52
skymap tau
Definition: detector.hh:346
double get_XX()
Definition: detector.hh:310
size_t ssize()
Definition: detector.hh:190
skymap mFx
Definition: detector.hh:348
istream & operator>>(istream &in, Time &time)
Definition: time.cc:303
wavearray< double > ffp
Definition: detector.hh:365
POLARIZATION getPolarization()
Definition: detector.hh:307
void clearFilter()
Definition: detector.hh:143
ostream & operator<<(ostream &out, Time &time)
Definition: time.cc:632
WSeries< double > nRMS
Definition: detector.hh:358
TF1 * f2
wavearray< double > TIME
Definition: detector.hh:375
r setFilter(10)
i() int(T_cor *100))
wavearray< double > * getHoT()
param: no parameters
Definition: detector.hh:175
wavearray< double > fpx
Definition: detector.hh:367
wavearray< double > fp
Definition: detector.hh:363
double get_XS()
Definition: detector.hh:312
double elevation
Definition: detector.hh:53
POLARIZATION polarization
Definition: detector.hh:384
std::vector< int > RWFID
Definition: detector.hh:381
std::vector< wavearray< double > * > IWFP
Definition: detector.hh:379
double F
Definition: skymap.hh:63
double latitude
Definition: detector.hh:51
wavearray< double > lagShift
Definition: detector.hh:369
void bandPass(double f1, double f2, double a=0.)
Definition: detector.hh:283
wavearray< double > TDUR
Definition: detector.hh:376
std::vector< int > IWFID
Definition: detector.hh:378
wavearray< double > BAND
Definition: detector.hh:374
regression r
Definition: Regression_H1.C:44
double ekXk
Definition: detector.hh:340
s s
Definition: cwb_net.C:155
char filter[1024]
wavearray< double > ffm
Definition: detector.hh:366
WSeries< double > TFmap
Definition: detector.hh:354
void sclear()
Definition: detector.hh:191
virtual void Browse(TBrowser *b)
Definition: detector.hh:302
WSeries< double > waveForm
Definition: detector.hh:355
double AltY
Definition: detector.hh:56
WSeries< double > * getTFmap()
param: no parameters
Definition: detector.hh:179
double xSNR
Definition: detector.hh:339
skymap mFp
Definition: detector.hh:347
double fabs(const Complex &x)
Definition: numpy.cc:55
wavecomplex d_complex
Definition: detector.hh:41
iD print()
double mTau
double get_NN()
Definition: detector.hh:311
std::vector< short > index
Definition: detector.hh:45
double null
Definition: detector.hh:336
wavearray< double > FREQ
Definition: detector.hh:373
double dT
Definition: testWDM_5.C:12
SSeries< double > * getSTFmap(size_t n)
Definition: detector.hh:187
POLARIZATION
Definition: detector.hh:60
std::vector< SSeries< double > > vSS
Definition: detector.hh:352
WSeries< double > waveNull
Definition: detector.hh:357
WSeries< float > nVAR
Definition: detector.hh:359
std::vector< float > value
Definition: detector.hh:46
double shift()
Definition: detector.hh:270
size_t ifoID
Definition: detector.hh:328
int getSTFind(double r)
Definition: detector.hh:182
wavearray< double > ISNR
Definition: detector.hh:372
size_t nDFL
Definition: detector.hh:343
init()
Definition: revMonster.cc:12
int wfsave()
Definition: detector.hh:316
void bandCut(double f1, double f2)
Definition: detector.hh:286
double enrg
Definition: detector.hh:337
bool isBuiltin()
Definition: detector.hh:304
void addSTFmap(netcluster *pwc, double mTau=0.042)
Definition: detector.hh:391