Logo coherent WaveBurst  
Library Reference Guide
Logo
CWB_Plugin_waveform.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 "TComplex.h"
32 #include "TMath.h"
33 #include "mdc.hh"
34 #include <vector>
35 #include "watplot.hh"
36 #include <stdio.h>
37 
38 
39 double GetBoundaries(wavearray<double> x, double P, double& bT, double& eT);
40 
41 void
43 //!EVENT_REPORT
44 // This plugin add to the output root file the whitened recunstructed waveforms
45 
46  cout << endl;
47  cout << "-----> plugins/CWB_Plugin_waveform.C" << endl;
48  cout << "ifo " << ifo.Data() << endl;
49  cout << "type " << type << endl;
50  cout << endl;
51 
52 
53  if(type==CWB_PLUGIN_OLIKELIHOOD) {
54 
55  int i,j,k,n,m,l;
56  int nIFO = NET->ifoListSize();
57  int K = NET->nLag;
58  int M = NET->mdc__ID.size();
59  int ID;
60  char search = NET->tYPe;
61 
63 
64  bool batch = gROOT->IsBatch();
65  gROOT->SetBatch(true);
66 
67  char ifostr[20] = "";
68  char strtime[1024] = "";
69  char fname[1024];
70 
71  int minTimeDet=nIFO;
72  double minTime=1e40;
73  double eventTime[NIFO];
74  double lagTime[NIFO];
75  int ifoid[NIFO],sortid[NIFO];
76  double factor = 1;
77  detector *pD[NIFO];
78 
79 
80  netevent* EVT = new netevent(nIFO);
81 
82  // search output root file in the system list
83  TFile* froot = NULL;
84  TList *files = (TList*)gROOT->GetListOfFiles();
85  if (files) {
86  TIter next(files);
87  TSystemFile *file;
88  TString fname;
89  bool check=false;
90  while ((file=(TSystemFile*)next())) {
91  fname = file->GetName();
92  // set output root file as the current file
93  if(fname.Contains("wave_")) {froot=(TFile*)file;froot->cd();}
94  }
95  if(!froot) {
96  cout << "plugins/CWB_Plugin_waveform.C : Error - output root file not found" << endl;
97  exit(1);
98  }
99  } else {
100  cout << "plugins/CWB_Plugin_waveform.C : Error - output root file not found" << endl;
101  exit(1);
102  }
103 
104  int ndim=nIFO;
105  int run=NET->nRun;
106  int eventID;
107  float rho,netcc,frequency,likelihood;
108  float snr[NIFO];
109  double time[NIFO];
110  int lag[NIFO+1];
111  int slag[NIFO+1];
113  for(n=0; n<nIFO; n++) wave[n] = new wavearray<double>;
114  TTree* wftree = (TTree *) froot->Get("waveform");
115  if(wftree==NULL) {
116  wftree = new TTree("waveform","waveform");
117  wftree->Branch("ndim",&ndim,"ndim/I");
118  wftree->Branch("run",&run,"run/I");
119  wftree->Branch("eventID",&eventID,"eventID/I");
120  wftree->Branch("rho",&rho,"rho/F");
121  wftree->Branch("netcc",&netcc,"netcc/F");
122  wftree->Branch("likelihood",&likelihood,"likelihood/F");
123  char csnr[16];sprintf(csnr,"snr[%1d]/F",nIFO);
124  wftree->Branch("snr",snr,csnr);
125  wftree->Branch("frequency",&frequency,"frequency/F");
126  char ctime[16];sprintf(ctime,"time[%1d]/D",nIFO);
127  wftree->Branch("time",time,ctime);
128  char clag[16];sprintf(clag,"lag[%1d]/F",nIFO+1);
129  char cslag[16];sprintf(cslag,"slag[%1d]/F",nIFO+1);
130  wftree->Branch("lag",lag,clag);
131  wftree->Branch("slag",slag,cslag);
132  char label[256];
133  for(n=0; n<nIFO; n++) {
134  sprintf(label, "%s_waveform", NET->ifoName[n]);
135  wftree->Branch(label,"wavearray<double>",&wave[n],32000,0);
136  }
137  } else {
138  wftree->SetBranchAddress("ndim",&ndim);
139  wftree->SetBranchAddress("run",&run);
140  wftree->SetBranchAddress("eventID",&eventID);
141  wftree->SetBranchAddress("rho",&rho);
142  wftree->SetBranchAddress("netcc",&netcc);
143  wftree->SetBranchAddress("frequency",&frequency);
144  wftree->SetBranchAddress("likelihood",&likelihood);
145  wftree->SetBranchAddress("snr",snr);
146  wftree->SetBranchAddress("time",time);
147  wftree->SetBranchAddress("lag",lag);
148  wftree->SetBranchAddress("slag",slag);
149  char label[256];
150  for(n=0; n<nIFO; n++) {
151  wave[n] = new wavearray<double>;
152  sprintf(label, "%s_waveform", NET->ifoName[n]);
153  wftree->SetBranchAddress(label,&wave[n]);
154  }
155  }
156 
157  for(n=0; n<nIFO; n++) pD[n] = NET->getifo(n);
158 
159  int rate = int(2*NET->getifo(0)->TFmap.resolution(0)+0.5);
160 
161  for(k=0; k<K; k++) { // loop over the lags
162 
163  id = NET->getwc(k)->get(const_cast<char*>("ID"), 0, 'L', rate);
164 
165  for(j=0; j<(int)id.size(); j++) { // loop over cluster index
166 
167  ID = size_t(id.data[j]+0.5); // cluster id
168 
169  eventID = ID;
170 
171  EVT->output(NULL,NET,factor,ID,k); // get reconstructed parameters
172  if(EVT->rho[1] < cfg->cedRHO) continue; // skip events under rho threshold
173 
174  rho = EVT->rho[1];
175  netcc = EVT->netcc[0];
176  frequency = EVT->frequency[0];
177  likelihood = EVT->likelihood;
178  for(n=0; n<nIFO; n++) snr[n] = EVT->sSNR[n];
179  for(n=0; n<nIFO; n++) time[n] = EVT->time[n];
180  for(n=0; n<=nIFO; n++) lag[n] = EVT->lag[n];
181  for(n=0; n<=nIFO; n++) slag[n] = EVT->slag[n];
182 
183  int masterDet=0;
184  int lagMin=2147483647;
185  for(n=0; n<nIFO; n++) if(EVT->lag[n]<lagMin) {lagMin=int(EVT->lag[n]);masterDet=n;}
186 
187  NET->likelihood(search, NET->acor, ID, k); // exec likelihood search
188 
189  // get event network time
190  double gps_start = EVT->time[masterDet]-EVT->duration[masterDet];
191  double gps_stop = EVT->time[masterDet]+EVT->duration[masterDet];
192 
193  // create a unique label
194  //for(n=0; n<nIFO; n++) sprintf(strtime, "%s_%.3f",strtime,EVT->start[n]);
195  for(n=0; n<nIFO; n++) sprintf(strtime, "%s_%d",strtime,(int)EVT->start[n]);
196 
197  // compute event time & lags time
198  for(n=0; n<nIFO; n++) eventTime[n]=(EVT->start[n]+EVT->stop[n])/2.;
199  for(n=0; n<nIFO; n++) minTime = (eventTime[n]<minTime) ? eventTime[n] : minTime;
200  for(n=0; n<nIFO; n++) lagTime[n]=eventTime[n]-minTime;
201  for(n=0; n<nIFO; n++) minTimeDet = (eventTime[n]<=minTime) ? n : minTimeDet;
202 
203  NET->getwave(ID, k, 'W');
205  for(n=0; n<nIFO; n++) {
206  w.start(0);
207  w.rate(pD[n]->waveForm.rate());
208  w.resize(pD[n]->waveForm.size());
209  for(int i=0;i<w.size();i++) w[i]=pD[n]->waveForm.data[i];
210 
211  double bT,eT;
212  GetBoundaries(w, 0.99, bT, eT);
213  //cout << w.start() << " " << bT << " " << eT << " " << eT-bT << endl;
214 
215  // select data which contains the 0.99 of energy
216  double dt = 1./w.rate();
217  int size = (eT-bT)/dt;
218  int os = bT/dt;
219  wave[n]->resize(size);
220  wave[n]->start(w.start());
221  wave[n]->rate(w.rate());
222  for(int i=0;i<size;i++) wave[n]->data[i]=w[i+os];
223  }
224 
225  wftree->Fill();
226  } // End loop on found events
227  }
228 
229  wftree->Write("", TObject::kOverwrite);
230 
231  gROOT->SetBatch(batch); // restore batch status
232 
233  jfile->cd();
234 
235  for(n=0; n<nIFO; n++) delete wave[n];
236  delete EVT;
237  }
238 
239  if(type==CWB_PLUGIN_CLOSE_JOB) {
240  }
241 
242  return;
243 }
244 
245 
246 double
247 GetBoundaries(wavearray<double> x, double P, double& bT, double& eT) {
248 
249  if(P<0) P=0;
250  if(P>1) P=1;
251 
252  int N = x.size();
253 
254  double E = 0; // signal energy
255  double avr = 0; // average
256  for(int i=0;i<N;i++) {avr+=i*x[i]*x[i]; E+=x[i]*x[i];}
257  int M=int(avr/E); // central index
258 
259  // search range which contains percentage P of the total energy E
260  int jB=0;
261  int jE=N-1;
262  double a,b;
263  double sum = ((M>=0)&&(M<N)) ? x[M]*x[M] : 0.;
264  for(int j=1; j<N; j++) {
265  a = ((M-j>=0)&&(M-j<N)) ? x[M-j] : 0.;
266  b = ((M+j>=0)&&(M+j<N)) ? x[M+j] : 0.;
267  if(a) jB=M-j;
268  if(b) jE=M+j;
269  sum += a*a+b*b;
270  if(sum/E > P) break;
271  }
272 
273  bT = x.start()+jB/x.rate();
274  eT = x.start()+jE/x.rate();
275 
276  return eT-bT;
277 }
278 
279 
std::vector< char * > ifoName
Definition: network.hh:609
detector * getifo(size_t n)
param: detector index
Definition: network.hh:436
CWB::config * cfg
double rho
#define NIFO
Definition: wat.hh:74
size_t nLag
Definition: network.hh:573
Float_t * rho
biased null statistics
Definition: netevent.hh:112
double M
Definition: DrawEBHH.C:13
Double_t * start
cluster duration = stopW-startW
Definition: netevent.hh:99
virtual void rate(double r)
Definition: wavearray.hh:141
Float_t * duration
max cluster time relative to segment start
Definition: netevent.hh:98
double cedRHO
Definition: config.hh:298
float factor
CWB run(runID)
bool getwave(size_t, size_t, char='W')
param: cluster ID param: delay index param: time series type return: true if time series are extracte...
Definition: network.cc:3576
wavearray< double > a(hp.size())
int n
Definition: cwb_net.C:28
netevent W & wave
TString("c")
size_t nRun
Definition: network.hh:572
int ID
Definition: TestMDC.C:70
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
double frequency
wavearray< double > get(char *name, size_t index=0, char atype='R', int type=1, bool=true)
param: string with parameter name param: index in the amplitude array, which define detector param: c...
Definition: netcluster.cc:2207
string search
Definition: cWB_conf.py:110
char ifostr[64]
Long_t size
int m
Definition: cwb_net.C:28
virtual void start(double s)
Definition: wavearray.hh:137
std::vector< size_t > mdc__ID
Definition: network.hh:615
int j
Definition: cwb_net.C:28
double rate
Definition: detector.hh:341
i drho i
#define N
double GetBoundaries(wavearray< double > x, double P, double &bT, double &eT)
char ifo[NIFO_MAX][8]
size_t ifoListSize()
Definition: network.hh:431
wavearray< double > w
Definition: Test1.C:27
#define nIFO
virtual size_t size() const
Definition: wavearray.hh:145
Float_t * frequency
GPS stop time of the cluster.
Definition: netevent.hh:102
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *NET, WSeries< double > *x, TString ifo, int type)
COHERENCE.
jfile
Definition: cwb_job_obj.C:43
i() int(T_cor *100))
network NET
Definition: cwb_dump_inj.C:30
TString label
Definition: MergeTrees.C:21
Double_t * stop
GPS start time of the cluster.
Definition: netevent.hh:100
int l
TIter next(twave->GetListOfBranches())
Float_t * lag
time between consecutive events
Definition: netevent.hh:84
char fname[1024]
int k
Float_t * netcc
effective correlated SNR
Definition: netevent.hh:113
double acor
Definition: network.hh:585
Float_t * slag
time lag [sec]
Definition: netevent.hh:85
void output(TTree *=NULL, network *=NULL, double=0., size_t=0, int=-1)
Definition: netevent.cc:1193
long likelihood(char='E', double=sqrt(2.), int=0, size_t=0, int=-1, bool=false)
Definition: network.cc:4446
TFile * froot
double dt
WSeries< double > TFmap
Definition: detector.hh:354
netevent EVT(itree, nifo)
Double_t * time
beam pattern coefficients for hx
Definition: netevent.hh:94
Float_t * sSNR
data-signal correlation Xk*Sk
Definition: netevent.hh:136
double resolution(int=0)
Definition: wseries.hh:155
Float_t likelihood
Definition: netevent.hh:124
Definition: Toolbox.hh:99
netcluster * getwc(size_t n)
param: delay index
Definition: network.hh:439
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
Long_t id
double ctime
virtual void resize(unsigned int)
Definition: wavearray.cc:463
int check
char tYPe
Definition: network.hh:588
detector ** pD
exit(0)
char os[1024]
double avr