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