Logo coherent WaveBurst  
Library Reference Guide
Logo
CWB_Plugin_HandleWhiteDataWDM.C
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 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 #define XIFO 4
20 
21 #pragma GCC system_header
22 
23 #include "cwb.hh"
24 #include "config.hh"
25 #include "network.hh"
26 #include "wavearray.hh"
27 #include "TString.h"
28 #include "TObjArray.h"
29 #include "TObjString.h"
30 #include "TRandom.h"
31 #include "Toolbox.hh"
32 #include "watplot.hh"
33 
34 
35 void
37 //!DATA_HANDLING
38 // Plugin to handle whitened input data after its WDM transform (only for 2G analysis !!!)
39 
40 // #define PLOT_WDM_TF // Plot WDM Scalogram
41 // #define DATA_COMPARISON // Compare data before/after WDM data handling
42  #define DATA_HANDLING // example of data handling
43  #define COPY_HANDLED_DATA // copy handled data to into the original x array
44 
45  cout << endl;
46  cout << "-----> CWB_Plugin_HandleWhiteDataWDM.C" << endl;
47  cout << "ifo " << ifo.Data() << endl;
48  cout << "type " << type << endl;
49  cout << endl;
50 
51  if(type==CWB_PLUGIN_WHITE) {
52 
53  if(TString(cfg->analysis)!="2G") {
54  cout << "CWB_Plugin_HandleWhiteDataWDM.C : Error - analysis type must be 2G !!!" << endl;
55  gSystem->Exit(1);
56  }
57 
58  wavearray<double> y = *x; // copy original data to y (contains the final handled data)
59  WSeries<double> w; // temporary array for data manipulation
60  WDM<double>* pwdm = NULL;
61 
62  for(int level=cfg->l_high; level>=cfg->l_low; level--) { // loop over levels
63 
64  int layers = level>0 ? 1<<level : 0; // get layers
65  cout << "level : " << level << " layers : " << layers << endl;
66 
67  pwdm = net->getwdm(layers+1); // get pointer to wdm transform
68  if(pwdm==NULL) {
69  cout << "CWB_Plugin_HandleWhiteDataWDM.C : Error - WDM not defined !!!" << endl;
70  gSystem->Exit(1);
71  }
72 
73  w.Forward(y,*pwdm); // apply wdm transformation
74 
75 #ifdef DATA_HANDLING
76  // example of data handling
77  double R = x->rate();
78 
79  double dF = R/(2*layers); // frequency resolution
80  double dT = layers/R; // time resolution
81 
82  cout << "WDM Decomposition Level : " << w.getLevel()
83  << " dF : " << dF << " Hz " << " dT : " << 1000*dT << " ms " <<endl;
84 
85  // Compute Energy
86  int M = w.getLevel();
87  int L = pwdm->maxLayer();
88  cout << "WDM Decomposition Level : " << M << " Layers " << L << endl;
89 
90  int mF = int(w.size()/pwdm->nSTS);
91  int nTC = w.size()/(M+1)/mF; // # of Time Coefficients
92  // Extract map00 & map90 : size = (M+1)*nTC
93  double* map00 = pwdm->pWWS;
94  double* map90 = map00 + (mF-1)*(M+1)*nTC;
95 
98  double E00=0;
99  double E90=0;
100  // this example show the corrispondence between map00,map90 and xx,XX
101  for(int j=0;j<M+1;j++) { // loop over the layers
102  w.getLayer(xx,j); // extract phase 00 @ layer j
103  w.getLayer(XX,-j); // extract phase 90 @ layer j
104  // for phase 00 -> (xx[i] = map00[i*(M+1)+j])
105  // for phase 90 -> (XX[i] = map90[i*(M+1)+j])
106  for(int i=0;i<xx.size();i++) {
107  //cout << j << " P00 " << xx[i] << " " << map00[i*(M+1)+j]
108  //<< " P90 " << XX[i] << " " << map90[i*(M+1)+j] << endl;
109  }
110  for(int i=0;i<xx.size();i++) E00+=xx[i]*xx[i];
111  for(int i=0;i<XX.size();i++) E90+=XX[i]*XX[i];
112  }
113  cout << "E00 = " << E00 << endl;
114  cout << "E90 = " << E90 << endl;
115 #endif
116 
117 #ifdef PLOT_WDM_TF
118  // Plot WDM Scalogram
119  watplot WTS(const_cast<char*>("wdm"));
120  //scalogram maps
121  double start = w.start();
122  double stop = w.start()+w.size()/w.rate();
123  double flow = 64;
124  double fhigh = 2048;
125  WTS.plot(&w, 2, start, stop,const_cast<char*>("COLZ"));
126  WTS.hist2D->GetYaxis()->SetRangeUser(flow, fhigh);
127  // dump spectrum
128  char fname[1024];
129  sprintf(fname,"%s_wdm_scalogram_%d_lev%d.root", ifo.Data(), int(w.start()),w.getLevel());
130  cout << endl << "Dump WDM Scalogram : " << fname << endl << endl;
131  WTS.canvas->Print(fname);
132 #endif
133 
134  w.Inverse(); // inverse transform
135  y = w; // copy manipulated data to y
136  }
137 
138 #ifdef DATA_COMPARISON
139  // Compare data before/after WDM data handling
141 
142  watplot plot(const_cast<char*>("plot"),200,20,800,500);
143  char gtitle[256];
144  sprintf(gtitle,"WDM : Original(black) - After WDM Forward/Inverse Tranform(Red)");
145  plot.gtitle(gtitle,"time(sec)","amplitude");
146  double start = z.start()+cfg->segEdge;
147  double stop = z.start()+z.size()/z.rate()-cfg->segEdge;
148  plot.goptions("alp", 1, start, stop);
149  // draw signal
150  z-=*x; // difference
151  //*x >> plot;
152  z >> plot;
153 
154  int nEdge = cfg->segEdge*y.rate();
155  double ez=0;
156  for(int i=nEdge;i<(int)y.size()-nEdge;i++) ez+=y[i]*y[i];
157  double ex=0;
158  for(int i=nEdge;i<(int)x->size()-nEdge;i++) ex+=x->data[i]*x->data[i];
159  cout << "Energy of Difference / Energy of Input Data " << ez/ex << endl;
160 
161  //TString gfile=ifo+"_WDM_data_comparison.png"; // save plot the difference to png file
162  TString gfile=ifo+"_WDM_data_comparison.root"; // save plot the difference to root file
163  plot >> gfile;
164 #endif
165 
166 #ifdef COPY_HANDLED_DATA
167  // copy handled data to into the original x array
168  for(int i=0;i<(int)x->size();i++) x->data[i]=y[i];
169 #endif
170  }
171 
172  return;
173 }
CWB::config * cfg
char analysis[8]
Definition: config.hh:117
void gtitle(TString title="", TString xtitle="", TString ytitle="")
Definition: watplot.cc:1256
double M
Definition: DrawEBHH.C:13
virtual void rate(double r)
Definition: wavearray.hh:141
wavearray< double > z
Definition: Test10.C:32
TString("c")
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
plot gtitle(gtitle,"frequency (Hz)","strain/#sqrt{Hz}")
int layers
DataType_t * pWWS
Definition: WaveDWT.hh:141
virtual void start(double s)
Definition: wavearray.hh:137
int j
Definition: cwb_net.C:28
#define COPY_HANDLED_DATA
i drho i
void plot(wavearray< double > &, char *=NULL, int=1, double=0., double=0., bool=false, float=0., float=0., bool=false, float=0., bool=false)
Definition: watplot.cc:150
char ifo[NIFO_MAX][8]
network ** net
NOISE_MDC_SIMULATION.
TH2F * hist2D
Definition: watplot.hh:193
virtual int maxLayer()
Definition: Wavelet.hh:89
WDM< double > * getwdm(size_t M)
param: number of wdm layers
Definition: network.hh:448
wavearray< double > w
Definition: Test1.C:27
double segEdge
Definition: config.hh:164
virtual size_t size() const
Definition: wavearray.hh:145
TCanvas * canvas
Definition: watplot.hh:192
int getLevel()
Definition: wseries.hh:109
x plot
jfile
Definition: cwb_job_obj.C:43
wavearray< double > xx
Definition: TestFrame1.C:11
unsigned long nSTS
Definition: WaveDWT.hh:143
int l_high
Definition: config.hh:156
watplot * WTS
Definition: ChirpMass.C:119
i() int(T_cor *100))
double fhigh
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
Definition: wseries.cc:193
char fname[1024]
void goptions(char *opt=NULL, int col=1, double t1=0., double t2=0., bool fft=false, float f1=0., float f2=0., bool psd=false, float t3=0., bool oneside=false)
Definition: watplot.cc:1221
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *net, WSeries< double > *x, TString ifo, int type)
COHERENCE.
double flow
int l_low
Definition: config.hh:155
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Definition: wseries.cc:246
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
TString gfile
DataType_t * data
Definition: wavearray.hh:319
double dT
Definition: testWDM_5.C:12
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
Definition: wseries.cc:291
wavearray< double > y
Definition: Test10.C:31