Logo coherent WaveBurst  
Library Reference Guide
Logo
CWB_Plugin_WRC.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 "watplot.hh"
35 #include "gwavearray.hh"
36 #include <vector>
37 
38 //#define SAVE_MRA_PLOT // enable event MRA plots
39 //#define SAVE_WHT_PLOT // enable event WHITE plots
40 //#define SAVE_STR_PLOT // enable event STRAIN plots
41 //#define SAVE_TF_PLOT // enable event WHITE TF plots
42 
43 #define MDC_OTF // enable MDC On The Fly
44 #define NOISE_OTF // enable NOISE On The Fly
45 
46 //#define SAVE_EVT_CLUSTER // save cluster to the event output tree
47 //#define SAVE_EVT_SNR // save iSNR,oSNR,ioSNR
48 
50  network* NET, netevent* EVT, int lag, int ID, int ifoID);
52  double& enINJ, double& enREC, double& xcorINJ_REC);
53 double ComputeEnergy(WSeries<double> *WS);
57  CWB::config* cfg, bool fft=false, bool strain=false);
58 void PlotMRA(netcluster* pwc, int nIFO, int ID);
59 
60 
61 void
63 //!NOISE_MDC_SIMULATION
64 // Extract whitened/strain injected & reconstructed waveforms, compute residual energy
65 // Save residual energy to the output wave root file
66 // This plugin can be use for Waveform Reconstruction Challenge (WRC) studies
67 
68  cout << endl;
69  cout << "-----> CWB_Plugin_WRC.C" << endl;
70  cout << "ifo " << ifo.Data() << endl;
71  cout << "type " << type << endl;
72  cout << endl;
73 
74  if(type==CWB_PLUGIN_CONFIG) {
75 #ifdef NOISE_OTF
76  cfg->dataPlugin=true; // disable read data from frames
77 #endif
78 #ifdef MDC_OTF
79  cfg->mdcPlugin=true; // disable read mdc from frames
80 #endif
81  cfg->outPlugin=true; // disable built-in output root file
82  }
83 
84 
85 #ifdef NOISE_OTF
86  if(type==CWB_PLUGIN_DATA) {
87 
89 
90  int seed;
91  if(ifo.CompareTo("L1")==0) seed=1000;
92  if(ifo.CompareTo("H1")==0) seed=2000;
93  if(ifo.CompareTo("V1")==0) seed=3000;
94  if(ifo.CompareTo("J1")==0) seed=4000;
95  if(ifo.CompareTo("A2")==0) seed=5000;
96  if(ifo.CompareTo("Y2")==0) seed=6000;
97  if(ifo.CompareTo("Y3")==0) seed=7000;
98 
99  TString fName;
100  if(ifo.CompareTo("L1")==0) fName="plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
101  if(ifo.CompareTo("H1")==0) fName="plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
102  if(ifo.CompareTo("V1")==0) fName="plugins/strains/advVIRGO_sensitivity_12May09_8khz_one_side.txt";
103  if(ifo.CompareTo("J1")==0) fName="plugins/strains/LCGT_sensitivity_8khz_one_side.txt";
104  if(ifo.CompareTo("A2")==0) fName="plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
105  if(ifo.CompareTo("Y2")==0) fName="plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
106  if(ifo.CompareTo("Y3")==0) fName="plugins/strains/advLIGO_NSNS_Opt_8khz_one_side.txt";
107 
108  int size=x->size();
109  double start=x->start();
110  TB.getSimNoise(*x, fName, seed, NET->nRun);
111  x->resize(size);
112  x->start(start);
113  }
114 #endif
115 
116 #ifdef MDC_OTF
117  if(type==CWB_PLUGIN_MDC) {
118 
119  char cmd[128];
120  sprintf(cmd,"network* net = (network*)%p;",NET);
121  gROOT->ProcessLine(cmd);
122 
123  CWB::mdc MDC(NET);
124 
125  // ---------------------------------
126  // read plugin config
127  // ---------------------------------
128 
129  cfg->configPlugin.Exec();
130 
131  // ---------------------------------
132  // set list of mdc waveforms
133  // ---------------------------------
134 
135  IMPORT(CWB::mdc,MDC)
136  MDC.Print();
137 
138  // ---------------------------------
139  // get mdc data
140  // ---------------------------------
141 
142  MDC.Get(*x,ifo);
143 
144  // ---------------------------------
145  // set mdc list in the network class
146  // ---------------------------------
147 
148  if(ifo.CompareTo(NET->ifoName[0])==0) {
149  NET->mdcList.clear();
150  NET->mdcType.clear();
151  NET->mdcTime.clear();
152  NET->mdcList=MDC.mdcList;
153  NET->mdcType=MDC.mdcType;
154  NET->mdcTime=MDC.mdcTime;
155  }
156 
157  cout.precision(14);
158  for(int k=0;k<(int)NET->mdcList.size();k++) cout << k << " mdcList " << MDC.mdcList[k] << endl;
159  for(int k=0;k<(int)NET->mdcTime.size();k++) cout << k << " mdcTime " << MDC.mdcTime[k] << endl;
160  for(int k=0;k<(int)NET->mdcType.size();k++) cout << k << " mdcType " << MDC.mdcType[k] << endl;
161  }
162 #endif
163 
164  if(type==CWB_PLUGIN_OLIKELIHOOD) {
165 
166  if(TString(cfg->analysis)!="2G") {
167  cout << "CWB_Plugin_WRC.C -> CWB_PLUGIN_OLIKELIHOOD implemented only for 2G" << endl;
168  gSystem->Exit(1);
169  }
170 
171  // import ifactor
172  int gIFACTOR=-1; IMPORT(int,gIFACTOR)
173  cout << "-----> CWB_Plugin_WRC.C -> " << " gIFACTOR : " << gIFACTOR << endl;
174 
175  // import slagShift
176  float* gSLAGSHIFT=NULL; IMPORT(float*,gSLAGSHIFT)
177 
178  int nIFO = NET->ifoListSize(); // number of detectors
179  int K = NET->nLag; // number of time lag
180  netevent* EVT;
182  double factor = cfg->simulation==3||cfg->simulation==4 ? -gIFACTOR : cfg->factors[gIFACTOR];
183  int rate = 0; // select all resolutions
184 
185  // search output root file in the system list
186  TFile* froot = NULL;
187  TList *files = (TList*)gROOT->GetListOfFiles();
188  TString outDump="";
189  if (files) {
190  TIter next(files);
191  TSystemFile *file;
192  TString fname;
193  bool check=false;
194  while ((file=(TSystemFile*)next())) {
195  fname = file->GetName();
196  // set output root file as the current file
197  if(fname.Contains("wave_")) {
198  froot=(TFile*)file;froot->cd();
199  outDump=fname;
200  outDump.ReplaceAll(".root.tmp",".txt");
201  cout << "output file name : " << fname << endl;
202  }
203  }
204  if(!froot) {
205  cout << "CWB_Plugin_WRC.C : Error - output root file not found" << endl;
206  gSystem->Exit(1);
207  }
208  } else {
209  cout << "CWB_Plugin_WRC.C : Error - output root file not found" << endl;
210  gSystem->Exit(1);
211  }
212 
213  netcluster* pwc = new netcluster();
214  char ciSNR[16]; sprintf(ciSNR,"_iSNR[%1d]/F",nIFO);
215  char coSNR[16]; sprintf(coSNR,"_oSNR[%1d]/F",nIFO);
216  char cioSNR[16]; sprintf(cioSNR,"_ioSNR[%1d]/F",nIFO);
217  float* iSNR = new float[nIFO];
218  float* oSNR = new float[nIFO];
219  float* ioSNR = new float[nIFO];
220 
221  TTree* net_tree = (TTree *) froot->Get("waveburst");
222  if(net_tree!=NULL) {
223  EVT = new netevent(net_tree,nIFO);
224 #ifdef SAVE_EVT_CLUSTER
225  net_tree->SetBranchAddress("netcluster",&pwc);
226 #endif
227 #ifdef SAVE_EVT_SNR
228  net_tree->SetBranchAddress("_iSNR",iSNR);
229  net_tree->SetBranchAddress("_oSNR",oSNR);
230  net_tree->SetBranchAddress("_ioSNR",ioSNR);
231 #endif
232  } else {
233  EVT = new netevent(nIFO);
234  net_tree = EVT->setTree();
235 #ifdef SAVE_EVT_CLUSTER
236  // add new netcluster leaf
237  net_tree->Branch("netcluster","netcluster",&pwc,32000,0);
238 #endif
239 #ifdef SAVE_EVT_SNR
240  net_tree->Branch("_iSNR",iSNR,ciSNR);
241  net_tree->Branch("_oSNR",oSNR,coSNR);
242  net_tree->Branch("_ioSNR",ioSNR,cioSNR);
243 #endif
244  }
245  EVT->setSLags(gSLAGSHIFT); // set slags into netevent
246 
247  for(int k=0; k<K; k++) { // loop over the lags
248 
249  id = NET->getwc(k)->get(const_cast<char*>("ID"), 0, 'L', rate);
250 
251  for(int j=0; j<(int)id.size(); j++) { // loop over cluster index
252 
253  int ID = size_t(id.data[j]+0.5);
254 
255  if(NET->getwc(k)->sCuts[ID-1]!=-1) continue; // skip rejected/processed clusters
256 
257  double ofactor=0;
258  if(cfg->simulation==4) ofactor=-factor;
259  else if(cfg->simulation==3) ofactor=-gIFACTOR;
260  else ofactor=factor;
261 
262  if(k==0) { // only for zero lag
263 
264  EVT->output2G(NULL,NET,ID,k,ofactor); // get reconstructed parameters
265 
266  double recTime = EVT->time[0]; // rec event time det=0
267  double injTh = EVT->theta[1]; // inj theta
268  double injPh = EVT->phi[1]; // inj phi
269  double recTh = EVT->theta[0]; // rec theta
270  double recPh = EVT->phi[0]; // rec phi
271 
272  injection INJ(nIFO);
273  // get inj ID
274  int M = NET->mdc__IDSize();
275  double injTime = 1.e12;
276  int injID = -1;
277  for(int m=0; m<M; m++) {
278  int mdcID = NET->getmdc__ID(m);
279  double T = fabs(recTime - NET->getmdcTime(mdcID));
280  if(T<injTime && INJ.fill_in(NET,mdcID)) {
281  injTime = T;
282  injID = mdcID;
283  }
284  }
285 
286  if(INJ.fill_in(NET,injID)) { // get inj parameters
287 
288  wavearray<double>** pwfINJ = new wavearray<double>*[nIFO];
289  wavearray<double>** pwfREC = new wavearray<double>*[nIFO];
290  detector* pd = NET->getifo(0);
291  int idSize = pd->RWFID.size();
292  int wfIndex=-1;
293  for(int mm=0; mm<idSize; mm++) if (pd->RWFID[mm]==ID) wfIndex=mm;
294 
295  // extract whitened injected & reconstructed waveforms
296  for(int n=0; n<nIFO; n++) {
297 
298  pd = NET->getifo(n);
299 
300  pwfINJ[n] = INJ.pwf[n];
301  if (pwfINJ[n]==NULL) {
302  cout << "CWB_Plugin_WRC.C : Error : Injected waveform not saved !!! : detector "
303  << NET->ifoName[n] << endl;
304  continue;
305  }
306  if (wfIndex<0) {
307  cout << "CWB_Plugin_WRC.C : Error : Reconstructed waveform not saved !!! : ID -> "
308  << ID << " : detector " << NET->ifoName[n] << endl;
309  continue;
310  }
311 
312  if (wfIndex>=0) pwfREC[n] = pd->RWFP[wfIndex];
313  double R = pd->TFmap.rate();
314  double rFactor = 1.;
315  rFactor *= factor;
316  wavearray<double>* wfINJ = pwfINJ[n]; // array of injected waveforms
317  *wfINJ*=rFactor; // injected wf is multiplied by the custom factor
318  wavearray<double>* wfREC = pwfREC[n]; // array of reconstructed waveforms
319 
320  double enINJ, enREC, xcorINJ_REC;
321  // compute residual energy in time domain
322  ComputeResidualEnergy(wfINJ, wfREC, enINJ, enREC, xcorINJ_REC);
323  // compute residual energy in TF domain at optimal resolution
324  ComputeResidualEnergyOptTF(wfINJ, wfREC, NET, EVT, k, ID, n);
325 #ifdef SAVE_WHT_PLOT
326  PlotWaveform(NET->ifoName[n], wfINJ, wfREC, cfg, false, false);
327  PlotWaveform(NET->ifoName[n], wfINJ, wfREC, cfg, true, false);
328 #endif
329 
330  iSNR[n] = enINJ; // inj whitened SNR^2 [fLow,fHigh]
331  oSNR[n] = enREC; // rec whitened SNR^2 [fLow,fHigh]
332  ioSNR[n] = xcorINJ_REC; // inj*rec whitened energy [fLow,fHigh]
333 
334  *wfINJ*=1./rFactor; // restore injected amplitude
335  }
336  delete [] pwfINJ;
337  delete [] pwfREC;
338  }
339 
340  // extract strain injected & reconstructed waveforms
341  if(INJ.fill_in(NET,-(injID+1))) { // set injections
342 
343  wavearray<double>** pwfINJ = new wavearray<double>*[nIFO];
344  wavearray<double>** pwfREC = new wavearray<double>*[nIFO];
345  detector* pd = NET->getifo(0);
346  int idSize = pd->RWFID.size();
347  int wfIndex=-1;
348  for(int mm=0; mm<idSize; mm++) if (pd->RWFID[mm]==-ID) wfIndex=mm;
349 
350  for(int n=0; n<nIFO; n++) {
351 
352  pd = NET->getifo(n);
353 
354  pwfINJ[n] = INJ.pwf[n];
355  if (pwfINJ[n]==NULL) {
356  cout << "CWB_Plugin_WRC.C : Error : Injected waveform not saved !!! : detector "
357  << NET->ifoName[n] << endl;
358  continue;
359  }
360  if (wfIndex<0) {
361  cout << "CWB_Plugin_WRC.C : Error : Reconstructed waveform not saved !!! : ID -> "
362  << ID << " : detector " << NET->ifoName[n] << endl;
363  continue;
364  }
365 
366  if (wfIndex>=0) pwfREC[n] = pd->RWFP[wfIndex];
367  double R = pd->TFmap.rate();
368  double rFactor = 1.;
369  rFactor *= factor;
370  wavearray<double>* wfINJ = pwfINJ[n]; // array of injected waveforms
371  *wfINJ*=rFactor; // injected wf is multiplied by the custom factor
372  wavearray<double>* wfREC = pwfREC[n]; // array of reconstructed waveforms
373 
374  ApplyFrequencyCuts(wfINJ, cfg);
375 #ifdef SAVE_STR_PLOT
376  PlotWaveform(NET->ifoName[n], wfINJ, wfREC, cfg, false, true);
377  PlotWaveform(NET->ifoName[n], wfINJ, wfREC, cfg, true, true);
378 #endif
379 
380  double enINJ, enREC, xcorINJ_REC;
381  // compute residual energy in time domain
382  ComputeResidualEnergy(wfINJ, wfREC, enINJ, enREC, xcorINJ_REC);
383  // compute residual energy in TF domain at optimal resolution
384  ComputeResidualEnergyOptTF(wfINJ, wfREC, NET, EVT, k, ID, n);
385 
386  iSNR[n] = enINJ; // inj whitened SNR^2 [fLow,fHigh]
387  oSNR[n] = enREC; // rec whitened SNR^2 [fLow,fHigh]
388  ioSNR[n] = xcorINJ_REC; // inj*rec whitened energy [fLow,fHigh]
389 
390  *wfINJ*=1./rFactor; // restore injected amplitude
391  }
392  delete [] pwfINJ;
393  delete [] pwfREC;
394  }
395  }
396 
397  std::vector<int> sCuts = NET->getwc(k)->sCuts; // save sCuts
398  // set sCuts=1 for the events which must be not copied with cps to pwc
399  for(int i=0; i<(int)sCuts.size(); i++) if(i!=ID-1) NET->getwc(k)->sCuts[i]=1;
400 
401  // after cpf, pwc contains only one event
402  // ID can not be used to get the event, to get event use ID=1 (ex: for watplot)
403  pwc->cpf(*(NET->getwc(k))); // note: likelihood2G delete tdAmp
404  NET->getwc(k)->sCuts = sCuts; // restore sCuts
405 
406  if(cfg->dump) EVT->dopen(outDump.Data(),const_cast<char*>("a"),false);
407  EVT->output2G(net_tree,NET,ID,k,ofactor); // get reconstructed parameters
408  if(cfg->dump) EVT->dclose();
409  if(!cfg->cedDump) NET->getwc(k)->sCuts[ID-1]=1; // mark as processed
410 
411 #ifdef SAVE_MRA_PLOT // monster event display
412  PlotMRA(pwc, nIFO, ID);
413 #endif
414  }
415  }
416 
417  jfile->cd();
418 
419  if(pwc) delete pwc;
420  if(iSNR) delete [] iSNR;
421  if(oSNR) delete [] oSNR;
422  if(ioSNR) delete [] ioSNR;
423  if(EVT) delete EVT;
424  }
425 
426  return;
427 }
428 
430  double& enINJ, double& enREC, double& xcorINJ_REC) {
431 
432  double bINJ = wfINJ->start();
433  double eINJ = wfINJ->stop();
434  double bREC = wfREC->start();
435  double eREC = wfREC->stop();
436  //cout.precision(14);
437  //cout << "bINJ : " << bINJ << " eINJ : " << eINJ << endl;
438  //cout << "bREC : " << bREC << " eREC : " << eREC << endl;
439 
440  double R=wfINJ->rate();
441 
442  int oINJ = bINJ>bREC ? 0 : int((bREC-bINJ)*R+0.5);
443  int oREC = bINJ<bREC ? 0 : int((bINJ-bREC)*R+0.5);
444  //cout << "oINJ : " << oINJ << " oREC : " << oREC << endl;
445 
446  double startXCOR = bINJ>bREC ? bINJ : bREC;
447  double endXCOR = eINJ<eREC ? eINJ : eREC;
448  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
449  //cout << "startXCOR : " << startXCOR << " endXCOR : "
450  // << endXCOR << " sizeXCOR :" << sizeXCOR << endl;
451 
452  enINJ=0;
453  enREC=0;
454  xcorINJ_REC=0;
455  if (sizeXCOR<=0) return;
456 
457  // the enINJ, enREC, xcorINJ_REC are computed in the INJ range
458  for (int i=0;i<wfINJ->size();i++) enINJ+=wfINJ->data[i]*wfINJ->data[i];
459  for (int i=0;i<sizeXCOR;i++) enREC+=wfREC->data[i+oREC]*wfREC->data[i+oREC];
460  for (int i=0;i<sizeXCOR;i++) xcorINJ_REC+=wfINJ->data[i+oINJ]*wfREC->data[i+oREC];
461 
462  double erINJ_REC = enINJ+enREC-2*xcorINJ_REC;
463  cout << "enINJ : " << enINJ << " enREC : " << enREC
464  << " xcorINJ_REC : " << xcorINJ_REC << " enINJREC : " << enINJ+enREC-2*xcorINJ_REC << endl;
465  //cout << "erINJ_REC/enINJ : " << erINJ_REC/enINJ << endl;
466 
467  return;
468 }
469 
471 
472  double f_low = cfg->fLow;
473  double f_high = cfg->fHigh;
474  cout << "f_low : " << f_low << " f_high : " << f_high << endl;
475  wf->FFTW(1);
476  double Fs=((double)wf->rate()/(double)wf->size())/2.;
477  for (int j=0;j<wf->size();j+=2) {
478  double f=j*Fs;
479  if((f<f_low)||(f>f_high)) {wf->data[j]=0.;wf->data[j+1]=0.;}
480  }
481  wf->FFTW(-1);
482  wf->resetFFTW();
483 
484  return;
485 }
486 
487 void PlotMRA(netcluster* pwc, int nIFO, int ID) {
488 
489  watplot WTS(const_cast<char*>("wts"));
490  WTS.canvas->cd();
491  char fname[1024];
492  sprintf(fname, "l_tfmap_scalogram_%d.png",ID);
493  cout << "write " << fname << endl;
494  WTS.plot(pwc, ID, nIFO, 'L', 0, const_cast<char*>("COLZ"));
495  WTS.canvas->Print(fname);
496  WTS.clear();
497 }
498 
500  CWB::config* cfg, bool fft, bool strain) {
501 
502  watplot PTS(const_cast<char*>("ptswrc"),200,20,800,500);
503 
504  //cout << "Print " << fname << endl;
505  double tmin = wfINJ->start()<wfREC->start() ? wfINJ->start() : wfREC->start();
506  wfINJ->start(wfINJ->start()-tmin);
507  wfREC->start(wfREC->start()-tmin);
508  if(fft) {
509  PTS.plot(wfINJ, const_cast<char*>("ALP"), 1, 0, 0, true, cfg->fLow, cfg->fHigh);
510  } else {
511  PTS.plot(wfINJ, const_cast<char*>("ALP"), 1, 0, 0);
512  }
513  PTS.graph[0]->SetLineWidth(1);
514  if(fft) {
515  PTS.plot(wfREC, const_cast<char*>("SAME"), 2, 0, 0, true, cfg->fLow, cfg->fHigh);
516  } else {
517  PTS.plot(wfREC, const_cast<char*>("SAME"), 2, 0, 0);
518  }
519  PTS.graph[1]->SetLineWidth(2);
520  wfINJ->start(wfINJ->start()+tmin);
521  wfREC->start(wfREC->start()+tmin);
522 
523  char label[64]="";
524  if(fft) sprintf(label,"%s","fft");
525  else sprintf(label,"%s","time");
526  if(strain) sprintf(label,"%s_%s",label,"strain");
527  else sprintf(label,"%s_%s",label,"white");
528 
529  char fname[1024];
530  sprintf(fname, "%s_wf_%s_inj_rec_gps_%d.root",ifo.Data(),label,int(tmin));
531  PTS.canvas->Print(fname);
532  cout << "write : " << fname << endl;
533  //PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
534 }
535 
537  network* NET, netevent* EVT, int lag, int ID, int ifoID) {
538 
539  double enINJ=0;
540  double enREC=0;
541  double enINJREC=0;
542 
543  // find TF optimal resolution level
544  netcluster* pwc = NET->getwc(lag);
545  double optRate = (pwc->cRate[ID-1])[0];
546  double optLayer = pwc->rate/optRate;
547  double optLevel = int(log10(optLayer)/log10(2));
548 
549  double bINJ = wfINJ->start();
550  double eINJ = wfINJ->stop();
551  double bREC = wfREC->start();
552  double eREC = wfREC->stop();
553  //cout.precision(14);
554  //cout << "bINJ : " << bINJ << " eINJ : " << eINJ << endl;
555  //cout << "bREC : " << bREC << " eREC : " << eREC << endl;
556 
557  double R=wfINJ->rate();
558 
559  int oINJ = bINJ>bREC ? 0 : int((bREC-bINJ)*R+0.5);
560  int oREC = bINJ<bREC ? 0 : int((bINJ-bREC)*R+0.5);
561  //cout << "oINJ : " << oINJ << " oREC : " << oREC << endl;
562 
563  double startXCOR = bINJ>bREC ? bINJ : bREC;
564  double endXCOR = eINJ<eREC ? eINJ : eREC;
565  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
566  //cout << "startXCOR : " << startXCOR << " endXCOR : "
567  // << endXCOR << " sizeXCOR :" << sizeXCOR << endl;
568 
569  if (sizeXCOR<=0) return;
570 
571  watplot WTS(const_cast<char*>("wtswrc"));
572  WTS.canvas->cd();
573 
574  detector* pd = NET->getifo(0);
575 
576  // length of temporary buffer for tf plots
577  int xcor_length = sizeXCOR/wfINJ->rate();
578  int wts_size = xcor_length<8 ? 16 : 16*int(1+xcor_length/8.);
579  wts_size*=wfINJ->rate();
580 
581  char fname[1024];
582  Meyer<double> S(1024,2); // set wavelet for production
583  double wts_start,wts_stop;
584 
585  // ------------------------------------------------------------------------
586  // INJECTED WAVEFORM
587  // ------------------------------------------------------------------------
588 
589  wavearray<double> xINJ(wts_size);
590  xINJ.start(wfINJ->start()-EVT->gps[0]+double(oINJ+wts_size/2)/wfINJ->rate());
591  xINJ.rate(wfINJ->rate());
592  xINJ=0.;
593  for(int m=0;m<sizeXCOR;m++)
594  if(m<(int)xINJ.size()/2) xINJ.data[m+wts_size/2]=wfINJ->data[m+oINJ];
595  WSeries<double> wINJ(xINJ,S);
596  xINJ.resize(0);
597 
598  if(NET->getwdm(optLayer+1)) wINJ.Forward(wINJ,*NET->getwdm(optLayer+1));
599 
600 #ifdef SAVE_TF_PLOT
601  //scalogram maps
602  sprintf(fname, "%s_wf_white_inj_tf.png", NET->ifoName[ifoID]);
603  cout << "Print " << fname << endl;
604  //PlotWSeries(&wINJ, fname);
605  wts_start = wINJ.start()+(double)(wts_size/2)/wINJ.rate();
606  wts_stop = sizeXCOR<wts_size/2 ? wts_start+sizeXCOR/wINJ.rate() : wts_start+(wts_size/2)/wINJ.rate();
607  WTS.plot(&wINJ, 2, wts_start, wts_stop,const_cast<char*>("COLZ"));
608  WTS.hist2D->GetYaxis()->SetRangeUser(EVT->low[ifoID], EVT->high[ifoID]);
609  WTS.canvas->Print(fname);
610  WTS.clear();
611 #endif
612 
613  // ------------------------------------------------------------------------
614  // RECONSTRUCTED WAVEFORM
615  // ------------------------------------------------------------------------
616 
617  wavearray<double> xREC(wts_size);
618  xREC.start(wfREC->start()-EVT->gps[0]+double(oREC+wts_size/2)/wfREC->rate());
619  xREC.rate(wfREC->rate());
620  xREC=0.;
621  for (int m=0;m<sizeXCOR;m++)
622  if(m<(int)xREC.size()/2) xREC.data[m+wts_size/2]=wfREC->data[m+oREC];
623  WSeries<double> wREC(xREC,S);
624  xREC.resize(0);
625 
626  if(NET->getwdm(optLayer+1)) wREC.Forward(wREC,*NET->getwdm(optLayer+1));
627 
628 #ifdef SAVE_TF_PLOT
629  //scalogram maps
630  sprintf(fname, "%s_wf_white_rec_tf.png", NET->ifoName[ifoID]);
631  cout << "Print " << fname << endl;
632  wts_start = wREC.start()+(double)(wts_size/2)/wREC.rate();
633  wts_stop = sizeXCOR<wts_size/2 ? wts_start+sizeXCOR/wREC.rate() : wts_start+(wts_size/2)/wREC.rate();
634  WTS.plot(&wREC, 2, wts_start, wts_stop,const_cast<char*>("COLZ"));
635  WTS.hist2D->GetYaxis()->SetRangeUser(EVT->low[ifoID], EVT->high[ifoID]);
636  WTS.canvas->Print(fname);
637  WTS.clear();
638 #endif
639 
640  // ------------------------------------------------------------------------
641  // INJECTED-RECONSTRUCTED WAVEFORM
642  // ------------------------------------------------------------------------
643 
644  wavearray<double> xDIF(wts_size);
645  xDIF.start(wfREC->start()-EVT->gps[0]+double(oREC+wts_size/2)/wfREC->rate());
646  xDIF.rate(wfREC->rate());
647  xDIF=0.;
648  for (int m=0;m<sizeXCOR;m++)
649  if(m<(int)xDIF.size()/2) xDIF.data[m+wts_size/2]=wfREC->data[m+oREC]-wfINJ->data[m+oINJ];
650  WSeries<double> wDIF(xDIF,S);
651  xDIF.resize(0);
652 
653 #ifdef SAVE_TF_PLOT
654  //scalogram maps
655  sprintf(fname, "%s_wf_white_dif_tf.png", NET->ifoName[ifoID]);
656  cout << "Print " << fname << endl;
657  if(NET->getwdm(optLayer+1)) wDIF.Forward(wDIF,*NET->getwdm(optLayer+1));
658  wts_start = wDIF.start()+(double)(wts_size/2)/wDIF.rate();
659  wts_stop = sizeXCOR<wts_size/2 ? wts_start+sizeXCOR/wDIF.rate() : wts_start+(wts_size/2)/wDIF.rate();
660  WTS.plot(&wDIF, 2, wts_start, wts_stop,const_cast<char*>("COLZ"));
661  WTS.hist2D->GetYaxis()->SetRangeUser(EVT->low[ifoID], EVT->high[ifoID]);
662  WTS.canvas->Print(fname);
663  WTS.clear();
664 #endif
665 
666  // compute the energy
667  enINJ = ComputeEnergy(&wINJ);
668  enREC = ComputeEnergy(&wREC);
669  enINJREC = ComputeResidualEnergy(&wINJ,&wREC);
670 
671  cout << "TF : " << "enINJ : " << enINJ << " enREC : " << enREC
672  << " enINJREC : " << enINJREC << endl;
673 
674  return;
675 }
676 
678 
679  int layers = WS->maxLayer()+1; // numbers of frequency bins (first & last bins have df/2)
680  int slices = WS->sizeZero(); // number of time bins
681 
682  float df = WS->resolution(); // frequency bin resolution (hz)
683  float dt = 1./(2*df); // time bin resolution (sec)
684 
685  int rate = int(1./dt);
686 
687 // cout << "layers : " << layers << "\t slices : " << slices << "\t rate : " << rate
688 // << "\t dt : " << dt << "\t df : " << df << endl;
689 
690  double EE=0;
691  double pixel_frequency;
692  double pixel_time;
693  for(int i=1;i<=slices;i++) {
694  pixel_time = i*dt;
695  for(int j=1;j<=layers;j++) {
696  if(j==1) pixel_frequency = df/2;
697  if(j==layers) pixel_frequency = (layers-1)*df;
698  if((j!=1)&&(j!=layers)) pixel_frequency = df/2 + (j-1)*df;
699  double ER = pow(WS->getSample(i-1,j-1),2);
700  if(j==1) EE += ER/2.;
701  if(j==layers) EE += ER/2.;
702  if((j!=1)&&(j!=layers)) EE += ER;
703  //cout << pixel_time << " : " << pixel_frequency << endl;
704  }
705  }
706 
707  return EE;
708 }
709 
711 
712  int layers = WS1->maxLayer()+1; // numbers of frequency bins (first & last bins have df/2)
713  int slices = WS1->sizeZero(); // number of time bins
714 
715  if((layers != WS2->maxLayer()+1)||(slices != WS2->sizeZero())) {
716  cout << "ComputeResidualEnergy - Error : WS1 & WS2 are not compatible !!!" << endl;
717  exit(1);
718  }
719 
720  float df = WS1->resolution(); // frequency bin resolution (hz)
721  float dt = 1./(2*df); // time bin resolution (sec)
722 
723  int rate = int(1./dt);
724 
725 // cout << "layers : " << layers << "\t slices : " << slices << "\t rate : " << rate
726 // << "\t dt : " << dt << "\t df : " << df << endl;
727 
728  double EE=0;
729  double pixel_frequency;
730  double pixel_time;
731  for(int i=1;i<=slices;i++) {
732  pixel_time = i*dt;
733  for(int j=1;j<=layers;j++) {
734  if(j==1) pixel_frequency = df/2;
735  if(j==layers) pixel_frequency = (layers-1)*df;
736  if((j!=1)&&(j!=layers)) pixel_frequency = df/2 + (j-1)*df;
737  double ER = pow(WS1->getSample(i-1,j-1)-WS2->getSample(i-1,j-1),2);
738  if(j==1) EE += ER/2.;
739  if(j==layers) EE += ER/2.;
740  if((j!=1)&&(j!=layers)) EE += ER;
741  //cout << pixel_time << " : " << pixel_frequency << endl;
742  }
743  }
744 
745  return EE;
746 }
std::vector< char * > ifoName
Definition: network.hh:609
detector * getifo(size_t n)
param: detector index
Definition: network.hh:436
CWB::config * cfg
char analysis[8]
Definition: config.hh:117
TString outDump
virtual void resize(unsigned int)
Definition: wseries.cc:901
size_t nLag
Definition: network.hh:573
std::vector< vector_int > cRate
Definition: netcluster.hh:398
Float_t * high
min frequency
Definition: netevent.hh:104
int slices
double M
Definition: DrawEBHH.C:13
std::vector< wavearray< double > > wREC[MAX_TRIALS]
void setSLags(float *slag)
Definition: netevent.cc:426
TMacro configPlugin
Definition: config.hh:362
TString Get(wavearray< double > &x, TString ifo)
Definition: mdc.cc:1529
std::vector< double > * getmdcTime()
Definition: network.hh:422
bool mdcPlugin
Definition: config.hh:365
std::vector< std::string > mdcList
Definition: mdc.hh:389
virtual void rate(double r)
Definition: wavearray.hh:141
Float_t * low
average center_of_snr frequency
Definition: netevent.hh:103
float factor
bool dataPlugin
Definition: config.hh:364
int n
Definition: cwb_net.C:28
TTree * setTree()
Definition: netevent.cc:434
TString("c")
std::vector< std::string > mdcType
Definition: mdc.hh:390
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
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
std::vector< wavearray< double > * > RWFP
Definition: detector.hh:382
double fLow
Definition: config.hh:140
netcluster * pwc
Definition: cwb_job_obj.C:38
CWB::Toolbox TB
std::vector< TGraph * > graph
Definition: watplot.hh:194
bool cedDump
Definition: config.hh:297
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *NET, WSeries< double > *x, TString ifo, int type)
COHERENCE.
int layers
std::vector< std::string > mdcType
Definition: network.hh:613
CWB::mdc * MDC
waveform wf
Long_t size
int m
Definition: cwb_net.C:28
std::vector< double > oSNR[NIFO_MAX]
virtual void start(double s)
Definition: wavearray.hh:137
int j
Definition: cwb_net.C:28
i drho i
double rate
Definition: netcluster.hh:378
std::vector< double > mdcTime
Definition: network.hh:614
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
std::vector< double > ioSNR[NIFO_MAX]
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
Definition: ComputeSNR.C:75
bool outPlugin
Definition: config.hh:369
void Print(int level=0)
Definition: mdc.cc:2736
char ifo[NIFO_MAX][8]
void dopen(const char *fname, char *mode, bool header=true)
Definition: netevent.hh:312
TH2F * hist2D
Definition: watplot.hh:193
size_t ifoListSize()
Definition: network.hh:431
WDM< double > * getwdm(size_t M)
param: number of wdm layers
Definition: network.hh:448
void clear()
Definition: watplot.hh:58
#define nIFO
virtual size_t size() const
Definition: wavearray.hh:145
TCanvas * canvas
Definition: watplot.hh:192
jfile
Definition: cwb_job_obj.C:43
void output2G(TTree *, network *, size_t, int, double)
Definition: netevent.cc:711
static void getSimNoise(wavearray< double > &u, TString fName, int seed, int run)
Definition: Toolbox.cc:5593
TTree * net_tree
Definition: mdc.hh:248
void ComputeResidualEnergyOptTF(wavearray< double > *wfINJ, wavearray< double > *wfREC, network *NET, netevent *EVT, int lag, int ID, int ifoID)
int simulation
Definition: config.hh:199
#define IMPORT(TYPE, VAR)
Definition: cwb.hh:69
watplot * WTS
Definition: ChirpMass.C:119
Double_t * gps
average center_of_gravity time
Definition: netevent.hh:95
i() int(T_cor *100))
network NET
Definition: cwb_dump_inj.C:30
TString label
Definition: MergeTrees.C:21
std::vector< std::string > mdcList
Definition: network.hh:612
TIter next(twave->GetListOfBranches())
char fname[1024]
std::vector< int > RWFID
Definition: detector.hh:381
size_t mdc__IDSize()
Definition: network.hh:414
std::vector< double > iSNR[NIFO_MAX]
void PlotWaveform(TString ifo, wavearray< double > *wfINJ, wavearray< double > *wfREC, CWB::config *cfg, bool fft=false, bool strain=false)
int k
virtual void resetFFTW()
Definition: wavearray.cc:977
size_t cpf(const netcluster &, bool=false, int=0)
Definition: netcluster.cc:117
void ComputeResidualEnergy(wavearray< double > *wfINJ, wavearray< double > *wfREC, double &enINJ, double &enREC, double &xcorINJ_REC)
TFile * froot
double dt
virtual void FFTW(int=1)
Definition: wavearray.cc:896
WSeries< double > TFmap
Definition: detector.hh:354
Float_t * phi
sqrt(h+*h+ + hx*hx)
Definition: netevent.hh:87
netevent EVT(itree, nifo)
Float_t * theta
[0]-reconstructed, [1]-injected phi angle, [2]-RA
Definition: netevent.hh:88
Double_t * time
beam pattern coefficients for hx
Definition: netevent.hh:94
double T
Definition: testWDM_4.C:11
std::vector< int > sCuts
Definition: netcluster.hh:392
void dclose()
Definition: netevent.hh:321
double fHigh
Definition: config.hh:141
Definition: Meyer.hh:36
virtual void stop(double s)
Definition: wavearray.hh:139
double fabs(const Complex &x)
Definition: numpy.cc:55
double resolution(int=0)
Definition: wseries.hh:155
void ApplyFrequencyCuts(wavearray< double > *wf, CWB::config *cfg)
int gIFACTOR
char cmd[1024]
Meyer< double > S(1024, 2)
double df
netcluster * getwc(size_t n)
param: delay index
Definition: network.hh:439
void PlotMRA(netcluster *pwc, int nIFO, int ID)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
DataType_t * data
Definition: wavearray.hh:319
Long_t id
double factors[FACTORS_MAX]
Definition: config.hh:202
DataType_t getSample(int n, double m)
Definition: wseries.hh:185
wavearray< double > wINJ[NIFO_MAX]
std::vector< double > mdcTime
Definition: mdc.hh:392
Bool_t fill_in(network *, int, bool=true)
Definition: injection.cc:349
size_t getmdc__ID(size_t n)
Definition: network.hh:426
char fName[256]
for(int i=0;i< 101;++i) Cos2[2][i]=0
pointers to detectors
wavearray< double > ** pwf
[x1,y1,z1,x2,y2,z2] components of spin vector
Definition: injection.hh:101
virtual void resize(unsigned int)
Definition: wavearray.cc:463
int check
size_t sizeZero()
Definition: wseries.hh:144
int maxLayer()
Definition: wseries.hh:139
exit(0)
double ComputeEnergy(WSeries< double > *WS)
bool dump
Definition: config.hh:295