Logo coherent WaveBurst  
Library Reference Guide
Logo
netcluster.hh
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Sergey Klimenko, Gabriele Vedovato, Valentin Necula
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 // Wavelet Analysis Tool
20 // Sergey Klimenko, University of Florida
21 // universal data container for network cluster analysis
22 // used with DMT and ROOT
23 //
24 
25 #ifndef NETCLUSTER_HH
26 #define NETCLUSTER_HH
27 
28 #include <iostream>
29 #include "wavearray.hh"
30 #include <vector>
31 #include <list>
32 #include "TH2F.h"
33 #include "TGraphErrors.h"
34 #include "TF1.h"
35 #include "WaveDWT.hh"
36 #include "wseries.hh"
37 #include "netpixel.hh"
38 #include "wat.hh"
39 #include "TNamed.h"
40 #include "TFile.h"
41 #include "TTree.h"
42 #include "WDM.hh"
43 
44 typedef std::vector<int> vector_int;
45 typedef std::vector<float> vector_float;
46 
47 class network;
48 
49 class clusterdata : public TNamed {
50 public:
53  float energy; // total cluster energy
54  float enrgsky; // cluster energy in all resolutions
55  float likenet; // signal energy
56  float netecor; // network coherent energy
57  float normcor; // normalized coherent energy
58  float netnull; // null energy in the sky loop with Gauss correction
59  float netED; // energy disbalance
60  float Gnoise; // estimated contribution of Gaussian noise
61  float likesky; // likelihood at all resolutions
62  float skycc; // network cc from the sky loop (all resolutions)
63  float netcc; // network cc for MRA or SRA analysis
64  float skyChi2; // chi2 stat from the sky loop (all resolutions)
65  float subnet; // first subNetCut statistic
66  float SUBNET; // second subNetCut statistic
67  float skyStat; // localization statistic
68  float netRHO; // coherent SNR per detector
69  float netrho; // reduced coherent SNR per detector
70  float theta; // source angle theta index
71  float phi; // source angle phi index
72  float iota; // inclination angle
73  float psi; // polarisation angle
74  float ellipticity; // waveform ellipticity
75  float cTime; // supercluster central time
76  float cFreq; // supercluster central frequency
77  float gNET; // network acceptance
78  float aNET; // network alignment
79  float iNET; // network index
80  float norm; // packet norm
81  float nDoF; // cluster degrees of freedom
82  float tmrgr; // merger time
83  float tmrgrerr; // merger time error
84  float mchirp; // chirp mass
85  float mchirperr; // chirp mass error
86  float chi2chirp; // chi2 over NDF
87  float chirpEfrac; // chirp energy fraction
88  float chirpPfrac; // chirp pixel fraction
89  float chirpEllip; // chirp ellipticity
90  int skySize; // number of sky pixels
91  int skyIndex; // index in the skymap
92  TF1 fit; //! chirp fit parameters (don't remove ! fix crash when exit from CINT)
93  TGraphErrors chirp; // chirp graph
94  wavearray<float> mchpdf; // chirp mass PDF
95  ClassDef(clusterdata,3)
96 };
97 
98 ClassImp(clusterdata)
99 
100 
101 class netcluster : public TNamed
102 {
103 public:
104 
105  // constructors
106 
107  //: Default constructor
108  netcluster();
109 
110  //: Copy constructor
111  //!param: value - object to copy from
112  netcluster(const netcluster&);
113 
114  //: destructor
115  virtual ~netcluster();
116 
117  // operators
118 
119  netcluster& operator= (const netcluster&);
120 
121  // copy function (used in = operator as x.cpf(y,false)
122  //: copy content of y into x.
123  // If no clusters are reconstructed - copy all pixels.
124  // If clusters are reconstructed - copy selected clusters.
125  // param: netcluster object y
126  // param: condition to select clusters:
127  // true - copy selected clusters with core pixels
128  // false - copy all selected clusters
129  // param: min size of BIG clusters
130  // =0 - copy all clusters regardles of their size
131  // >0 - ignore BIG clusters
132  // <0 - copy BIG clusters
133  size_t cpf(const netcluster &, bool=false, int=0);
134 
135 
136  // accessors
137 
138  //: clear lists and release memory;
139  inline void clear();
140  //: clean time delay data
141  inline void clean(int cID=0);
142  //: clean input wseries ws removing rejected clusters
143  //: output number of remaining pixels
144  int clean(WSeries<double>& ws);
145 
146  //: Get pixel list size
147  inline size_t size() { return pList.size(); }
148  //: Get pixel list capacity
149  inline size_t capacity() { return pList.capacity(); }
150  //: Get cluster list size
151  inline size_t csize() { return cList.size(); }
152  //: Get number of clusters with specified type
153  inline size_t esize(int k=2) {
154  size_t n=0;
155  for(size_t i=0; i<sCuts.size(); i++) {
156  if(k<2 && sCuts[i]==k) n++;
157  if(k>1 && sCuts[i] <1) n++;
158  }
159  return n;
160  }
161 
162  //: Get number of selected pixels
163  inline size_t psize(int k=2) {
164  size_t n=0;
165  for(size_t i=0; i<sCuts.size(); i++) {
166  if(k>1 && sCuts[i] <1) n+=cList[i].size();
167  if(k<2 && sCuts[i]==k) n+=cList[i].size();
168  }
169  return n;
170  }
171 
172  //: Get pixel pointer
173  //:parameter n - cluster ID
174  //:parameter i - pixel index in cList[n]
175  inline netpixel* getPixel(size_t n, size_t i);
176 
177  //: set black pixel probability
178  inline void setbpp(double P) { bpp = P; return; }
179  //: get black pixel probability
180  inline double getbpp() { return bpp; }
181 
182  //: set pixel core fields to be equal to the input parameter
183  // core: true/false
184  // id : set core for selected cluster id
185  size_t setcore(bool core,int id=0);
186 
187  //: set selection cuts vector used in mask(), occupancy(), getCluster()
188  //!param: cluster ID number
189  //!return void
190  inline void ignore(size_t n) {
191  if(n>0 && n<=sCuts.size()) sCuts[n-1] = 1;
192  }
193 
194  //: set sCuts vector excluding rejected clusters
195  //!param: sCuts flag
196  inline void setcuts(int n=0) {
197  for(size_t i=0; i<sCuts.size(); i++)
198  if(sCuts[i]!=1) sCuts[i] = n;
199  }
200 
201  //: remove halo pixels from pixel list
202  //!param: if false - de-cluster pixels
203  //!return size of the list
204  virtual size_t cleanhalo(bool=false);
205 
206  //: add halo pixels matching the paccket pattern to the pixel list
207  //!param: packet pattern mode
208  //!return size of the list
209  size_t addhalo(int=0);
210 
211  //: append pixel list from input cluster list
212  //: cluster metadata is lost
213  //!param: input netcluster
214  //!return size of appended pixel list
215  virtual size_t append(netcluster& wc);
216 
217  //: add netpixel object to the list
218  inline void append(netpixel& p) { pList.push_back(p); }
219 
220  // destroy supercluster neighbors links (delink superclusters)
221  // preserve links for pixels with the same wavelet resolution
222  virtual size_t delink();
223 
224  // select clusters
225  // name - parameter name
226  // thr - parameter threshold
227  virtual wavearray<double> select(char*, double);
228 
229  //:reconstruct clusters at several TF resolutions (superclusters)
230  //!param: statistic: E - excess power, L - likelihood
231  //!param: selection threshold T
232  //! for likelihood clusters, T defines a threshold on clusters
233  //! in a superclusters.
234  //!param: true - use only core pixels, false - use core & halo pixels
235  //!return size of pixel list of selected superclusters.
236  virtual size_t supercluster(char atype, double S, bool core); // used in 1G pipeline
237  virtual size_t supercluster(char atype, double S, double gap, bool core, TH1F* = NULL);
238 
239  //:merge clusters if they are close to each other
240  //! T - maximum time gap in seconds
241  //! F - maximum frequency gap in Hz
242  virtual size_t defragment(double T, double F, TH2F* = NULL);
243 
244  void PlotClusters();
245 
246  //: set clusterID field for pixels in pList vector
247  //: create cList structure - list of references to cluster's pixels
248  //!return number of clusters
249  virtual size_t cluster();
250 
251  //: recursively calculate clusterID pixels in pList vector
252  //!param: pixel pointer in pList vector
253  //!return cluster volume (total number of pixels)
254  virtual size_t cluster(netpixel* p);
255 
256  //:produce time-frequency clusters
257  //:sort pixels on index, form pixel groups,
258  //:set neighbor links for each group and cluster them
259  //!param: time gap between pixels in units of pixels
260  //!param: frequenct gap between pixels in units of pixels
261  //!return pixel number of time clusters
262  virtual size_t cluster(int kt,int kf);
263 
264  //: access function to get cluster parameters passed selection cuts
265  //!param: string with parameter name
266  //!param: index in the amplitude array, which define detector
267  //!param: character identifier for amplitude vector:
268  // 'W'-wavelet, 'S'-snr, 'R'-rank
269  //!param: rate index, if 0 ignore rate for calculation of cluster parameters
270  // if negative - extract clusters with rate = abs(rate)
271  //!return wavearray object with parameter values for clusters
272  wavearray<double> get(char* name, size_t index=0, char atype='R', int type=1, bool=true);
273 
274  //: extract WSeries for defined cluster ID
275  //!param: cluster ID
276  //!param: WSeries where to put the waveform
277  //!return: noise rms
278  double getwave(int, WSeries<double>&, char='W', size_t=0);
279 
280  // extract MRA waveforms for network net, cluster ID and detector index n
281  // works only with WDM. Create WSeries<> objects for each resolution,
282  // find principle components, fill in WSeries<>, Inverse
283  // construct waveform from MRA pixels at different resolutions
284  // atype = 'W' - get whitened detector output (Wavelet data)
285  // atype = 'w' - get detector output (Wavelet data)
286  // atype = 'S' - get whitened reconstructed response (Signal)
287  // atype = 's' - get reconstructed response (Signal)
288  // mode: -1/0/1 - return 90/mra/0 phase
289  wavearray<double> getMRAwave(network* net, int ID, size_t n, char atype='S', int mode=0);
290 
291  // write pixel structure into binary file
292  // only metadata and pixels are written
293  // first parameter is file name
294  // second parameter is the append mode 0/1 - wb/ab - new/append
295  size_t write(const char* file, int app=0);
296 
297  // write pixel structure with TD vectors attached into a file
298  // only some metadata and pixels are written, no cluster metadata is stored
299  // fp - file pointer
300  // app = 0 - store metadata
301  // app = 1 - store pixels with the TD vectors attached by setTDAmp()
302  size_t write(FILE* fp, int app=0);
303 
304  // write pixel structure with TD vectors attached into a root file
305  // froot - root file pointer
306  // tdir - internal root directory where the tree is located
307  // tname - name of tree containing the cluster
308  // app = 0 - store light netcluster
309  // app = 1 - store pixels with the TD vectors attached by setTDAmp()
310  // cycle - sim -> it is the factor id : prod -> it is the lag number
311  // irate - wavelet layer irate
312  // if irate is negative the value 'irate' is used to build the tree name
313  // this permits to create a tree for each irate
314  // cID - cluster id (cID=0 -> write all clusters)
315  size_t write(TFile *froot, TString tdir, TString tname, int app=0, int cycle=0, int irate=0, int cID=0);
316 
317  // read entire pixel structure from binary file
318  // returns number of pixels
319  size_t read(const char *);
320 
321  // read metadata and pixels stored in a file on cluster by cluster basis
322  // clusters should be contiguous in the file
323  // maxPix = 0 - read metadata
324  // maxPix > 0 - do not store TD vectors for clusters with size > maxPix
325  // maxPix < 0 - do not store TD vectors for clusters with size < maxPix
326  // clusters with no TD vectors attached are marked rejected.
327  // returns # of pixel in the cluster
328  size_t read(FILE* file, int maxPix);
329 
330  // read metadata and pixels stored in a file on cluster by cluster basis
331  // clusters should be contiguous in the file (written by write(FILE*))
332  // froot - root file pointer
333  // tdir - internal root directory where the tree is located
334  // tname - name of tree containing the cluster
335  // nmax = 0 - read metadata
336  // nmax > 0 - read no more than maxPix pixels from a cluster
337  // nmax < 0 - read all heavy pixels from a cluster
338  // nmax -2 - as for (nmax<0) & skip heavy instructions to speedup read
339  // cycle - sim -> it is factor id : prod -> it is the lag number
340  // rate - wavelet layer rate
341  // if rate is negative the value 'rate' is used to build the tree name
342  // cID - cluster ID (=0 : read all clusters)
343  std::vector<int> read(TFile* froot, TString tdir, TString tname, int nmax=0, int cycle=0, int rate=0, int cID=0);
344 
345  inline void setlow(double f) { flow=f; return; }
346  inline void sethigh(double f) { fhigh=f; return; }
347  inline double getlow() { return flow; }
348  inline double gethigh() { return fhigh; }
349 
350  // set arrays for time-delayed amplitudes in collected coherent pixels
351  // returns number of pixels to process: if count=0 - nothing to process
352  // param: input network
353  // param: 'a','A' - delayed amplitudes, 'e','E' - delayed power
354  // param: number of pixels (per resolution) to setup
355  size_t loadTDamp(network& net, char c, size_t BATCH=10000, size_t LOUD=0);
356  // fast version which uses sparse TF maps in detector class
357  size_t loadTDampSSE(network& net, char c, size_t BATCH=10000, size_t LOUD=0);
358 
359  // param 1 - cluster ID
360  // param 2 - chi2 threshold for selection of chirp pixels
361  // param 3 - tmerger cut: exclude pixels with time > tmerger+param3 - special use
362  // param 4 - threshold for pixel selection - special use
363  // returns reconstructed chirp energy
364  double mchirp(int ID, double=2.5, double=1.e20, double=0);
365 
366  // minipixels; using only P.C ; use df = sqrt(0.6* fdot)
367  double mchirpTEST(int ID);
368 
369  // draw chirp cluster
370  void chirpDraw(int id);
371 
372  // print detector parameters
373  void print(); // *MENU*
374  virtual void Browse(TBrowser *b) {print();}
375 
376  // data members
377 
378  double rate; // original Time series rate
379  double start; // interval start GPS time
380  double stop; // interval stop GPS time
381  double bpp; // black pixel probability
382  double shift; // time shift
383  double flow; // low frequency boundary
384  double fhigh; // high frequency boundary
385  size_t nPIX; // minimum number of pixels at all resolutions
386  int run; // run ID
387  bool pair; // true - 2 resolutions, false - 1 resolution
388  double nSUB; // subnetwork threshold for a single network pixel
389 
390  std::vector<netpixel> pList; // pixel list
391  std::vector<clusterdata> cData; // cluster metadata
392  std::vector<int> sCuts; /* cluster selection flags (cuts)
393  1 - rejected
394  0 - not processed / accepted
395  -1 - not complete
396  -2 - ready for processing */
397  std::vector<vector_int> cList; // cluster list defined by vector of pList references
398  std::vector<vector_int> cRate; // cluster type defined by rate
399  std::vector<float> cTime; // supercluster central time
400  std::vector<float> cFreq; // supercluster central frequency
401  std::vector<vector_float> sArea; // sky error regions
402  std::vector<vector_float> p_Map; // sky pixel map
403  std::vector<vector_int> p_Ind; // sky pixel index
404  std::vector<vector_int> nTofF; // sky time delay configuration for waveform backward correction
405 
406  ClassDef(netcluster,4)
407 
408 }; // class netcluster
409 
410 //: compare function to sort pixel objects on time
411 int compare_PIX(const void*, const void*);
412 
413 inline netpixel* netcluster::getPixel(size_t n, size_t i) {
414  if(!n) return &pList[i];
415  if(cList.size()<n) {
416  printf("getPixel(): size=%d, ID=%d\n",(int)cList.size(),(int)n);
417  return NULL;
418  }
419  if(cList[n-1].size()<=i) {
420  printf("getPixel(): size=%d, index=%d\n",(int)cList[n-1].size(),(int)i);
421  return NULL;
422  }
423  return &pList[cList[n-1][i]];
424 }
425 
426 // release vector memory
427 inline void netcluster::clear() {
428  while(!pList.empty()) pList.pop_back();
429  pList.clear(); std::vector<netpixel>().swap(pList);
430  while(!cList.empty()) cList.pop_back();
431  cList.clear(); std::vector<vector_int>().swap(cList);
432  while(!cData.empty()) cData.pop_back();
433  cData.clear(); std::vector<clusterdata>().swap(cData);
434  while(!sCuts.empty()) sCuts.pop_back();
435  sCuts.clear(); std::vector<int>().swap(sCuts);
436  while(!cRate.empty()) cRate.pop_back();
437  cRate.clear(); std::vector<vector_int>().swap(cRate);
438  while(!cTime.empty()) cTime.pop_back();
439  cTime.clear(); std::vector<float>().swap(cTime);
440  while(!cFreq.empty()) cFreq.pop_back();
441  cFreq.clear(); std::vector<float>().swap(cFreq);
442  while(!sArea.empty()) sArea.pop_back();
443  sArea.clear(); std::vector<vector_float>().swap(sArea);
444  while(!p_Map.empty()) p_Map.pop_back();
445  p_Map.clear(); std::vector<vector_float>().swap(p_Map);
446  while(!p_Ind.empty()) p_Ind.pop_back();
447  p_Ind.clear(); std::vector<vector_int>().swap(p_Ind);
448 }
449 
450 // clean time delay data
451 inline void netcluster::clean(int cID) {
452  for(int i=0; i<(int)this->cList.size(); ++i) { // loop over clusters
453  const vector_int& vint = cList[i];
454  if((cID!=0)&&((int)pList[vint[0]].clusterID!=cID)) continue;
455  for(int j=0; j<(int)vint.size(); j++) { // loop over pixels
456  if(pList[vint[j]].tdAmp.size()) pList[vint[j]].clean();
457  }
458  }
459 }
460 
461 #endif // NETCLUSTER_HH
462 
463 
464 
465 
466 
467 
468 
469 
470 
471 
472 
473 
474 
475 
476 
477 
478 
479 
void setbpp(double P)
Definition: netcluster.hh:178
float SUBNET
Definition: netcluster.hh:66
double stop
Definition: netcluster.hh:380
std::vector< int > vector_int
Definition: cluster.hh:35
std::vector< vector_int > cRate
Definition: netcluster.hh:398
int compare_PIX(const void *, const void *)
Definition: netcluster.cc:54
float likenet
Definition: netcluster.hh:55
float gNET
Definition: netcluster.hh:77
int irate
std::vector< pixel > cluster
Definition: wavepath.hh:61
void append(netpixel &p)
Definition: netcluster.hh:218
par [0] name
float energy
Definition: netcluster.hh:53
int n
Definition: cwb_net.C:28
double fhigh
Definition: netcluster.hh:384
TString("c")
float netED
Definition: netcluster.hh:59
float chirpEfrac
Definition: netcluster.hh:87
int ID
Definition: TestMDC.C:70
float mchirp
Definition: netcluster.hh:84
float tmrgrerr
Definition: netcluster.hh:83
double bpp
Definition: test_config1.C:22
float chirpPfrac
Definition: netcluster.hh:88
float normcor
Definition: netcluster.hh:57
std::vector< vector_float > sArea
Definition: netcluster.hh:401
std::vector< int > vector_int
Definition: netcluster.hh:44
double getlow()
Definition: netcluster.hh:347
double flow
Definition: netcluster.hh:383
void ignore(size_t n)
param: cluster ID number return void
Definition: netcluster.hh:190
double & gap
float tmrgr
Definition: netcluster.hh:82
Long_t size
std::vector< vector_int > cList
Definition: netcluster.hh:397
int j
Definition: cwb_net.C:28
i drho i
double rate
Definition: netcluster.hh:378
wc clear()
size_t capacity()
Definition: netcluster.hh:149
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
Definition: ComputeSNR.C:75
std::vector< vector_float > p_Map
Definition: netcluster.hh:402
void setlow(double f)
Definition: netcluster.hh:345
network ** net
NOISE_MDC_SIMULATION.
float likesky
Definition: netcluster.hh:61
void clean(int cID=0)
Definition: netcluster.hh:451
size_t mode
nc append(pix)
size_t esize(int k=2)
Definition: netcluster.hh:153
int BATCH
float ellipticity
Definition: netcluster.hh:74
std::vector< float > vector_float
Definition: netcluster.hh:45
double shift
Definition: netcluster.hh:382
void sethigh(double f)
Definition: netcluster.hh:346
double getbpp()
Definition: netcluster.hh:180
i() int(T_cor *100))
std::vector< netpixel > pList
Definition: netcluster.hh:390
double fhigh
float enrgsky
Definition: netcluster.hh:54
double start
Definition: netcluster.hh:379
printf("total live time: non-zero lags = %10.1f \, liveTot)
std::vector< vector_int > p_Ind
Definition: netcluster.hh:403
float skyStat
Definition: netcluster.hh:67
std::vector< float > cFreq
Definition: netcluster.hh:400
float netcc
Definition: netcluster.hh:63
int k
float chi2chirp
Definition: netcluster.hh:86
double F
wavearray< float > mchpdf
Definition: netcluster.hh:94
size_t size()
Definition: netcluster.hh:147
void setcuts(int n=0)
param: sCuts flag
Definition: netcluster.hh:196
virtual void Browse(TBrowser *b)
Definition: netcluster.hh:374
double gethigh()
Definition: netcluster.hh:348
TFile * froot
netpixel * getPixel(size_t n, size_t i)
Definition: netcluster.hh:413
float iota
Definition: netcluster.hh:72
float norm
Definition: netcluster.hh:80
float aNET
Definition: netcluster.hh:78
double flow
float mchirperr
Definition: netcluster.hh:85
float Gnoise
Definition: netcluster.hh:60
float iNET
Definition: netcluster.hh:79
double T
Definition: testWDM_4.C:11
TGraphErrors chirp
chirp fit parameters (don&#39;t remove ! fix crash when exit from CINT)
Definition: netcluster.hh:93
std::vector< int > sCuts
Definition: netcluster.hh:392
float netRHO
Definition: netcluster.hh:68
std::vector< float > cTime
Definition: netcluster.hh:399
int LOUD
double nSUB
Definition: netcluster.hh:388
wavearray< int > index
size_t csize()
Definition: netcluster.hh:151
iD print()
Meyer< double > S(1024, 2)
float subnet
Definition: netcluster.hh:65
float cFreq
Definition: netcluster.hh:76
float skycc
Definition: netcluster.hh:62
std::vector< clusterdata > cData
Definition: netcluster.hh:391
float cTime
Definition: netcluster.hh:75
float netnull
Definition: netcluster.hh:58
float skyChi2
Definition: netcluster.hh:64
netcluster wc
float nDoF
Definition: netcluster.hh:81
double bpp
Definition: netcluster.hh:381
size_t nPIX
Definition: netcluster.hh:385
size_t psize(int k=2)
Definition: netcluster.hh:163
float chirpEllip
Definition: netcluster.hh:89
float netecor
Definition: netcluster.hh:56
float theta
Definition: netcluster.hh:70
std::vector< vector_int > nTofF
Definition: netcluster.hh:404
float netrho
Definition: netcluster.hh:69
void clear()
Definition: netcluster.hh:427