Logo coherent WaveBurst  
Library Reference Guide
Logo
CWB_Plugin_xWRC.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 "TF2.h"
33 #include "TMath.h"
34 #include "mdc.hh"
35 #include "watplot.hh"
36 #include "gwavearray.hh"
37 #include "gskymap.hh"
38 #include <vector>
39 #include <iostream>
40 #include <algorithm>
41 
42 using namespace std;
43 
44 //#define SAVE_WHT_PLOT // enable event WHITE plots
45 //#define SAVE_WF_PLOT
46 //#define SAVE_SKYPROB
47 //#define SAVE_PLOT_WERR
48 
49 //#define SAVE_PLOT_WERR2
50 //#define SAVE_PLOT_RESIDUALS
51 //#define SAVE_MRA_PLOT // enable event MRA plots
52 //#define SAVE_WF_COMPARISON
53 
54 #define PLOT_DIR "plot"
55 
56 // ----------------------------------------------------
57 // WRC plugin functions
58 // ----------------------------------------------------
59 
61  double& enINJ, double& enREC, double& xcorINJ_REC);
62 
64 void DumpOutput(network* NET, netevent* &EVT, CWB::config* cfg, int ID, int k, int factor);
65 void GetRecWaveform(network* NET, netevent* EVT, int ID, int k, int factor);
66 void GetInjWaveform(network* NET, netevent* EVT, int ID, int k, int factor);
68 void ReadUserOptions();
69 void Clear();
70 void SetSkyMask(network* NET, CWB::config* cfg, int seed=0);
72 void Wave2Sparse(network* NET, CWB::config* cfg, char wave_type);
73 void SaveSkyProb(network* NET, CWB::config* cfg, int id);
76 double GetSparseMap(SSeries<double>* SS, bool phase, int index);
77 int _sse_mra_ps(network* NET, float* amp, float* AMP, float Eo, int K);
78 
79 void PlotFinal(int n);
80 void PlotFinal2(network* NET, int ID);
81 void PlotFinal3(network* NET, int ID, int ifoID, wavearray<double>* w1,
82  wavearray<double>* w2, TString tag, TString gtype="png");
84  CWB::config* cfg, bool fft=false, bool strain=false);
85 void PlotMRA(network* NET, int ID, int k, TString tag);
86 
87 // ----------------------------------------------------
88 // Temporary WRC structures
89 // ----------------------------------------------------
90 
92 wavearray<double> sINJ[NIFO_MAX]; // injected aligned to the same time
93 wavearray<double> rINJ[NIFO_MAX]; // injected waveform in wREC TF sub-domain
94 wavearray<double> wREC[NIFO_MAX]; // reconstructed
95 double wDINJ[NIFO_MAX]; // inj delays
96 double wDREC[NIFO_MAX]; // inj delays
97 
98 std::vector<SSeries<double> > vSS[NIFO_MAX]; // original sparse maps
99 std::vector<SSeries<double> > sSS[NIFO_MAX]; // injected aligned to the same time
100 std::vector<SSeries<double> > rSS[NIFO_MAX]; // reconstructed
101 std::vector<SSeries<double> > jSS[NIFO_MAX]; // injected
102 std::vector<WSeries<double> > vN[NIFO_MAX]; // sim noise
103 
104 std::vector<int> wTRY[NIFO_MAX];
105 std::vector<double> wAVR[NIFO_MAX];
106 std::vector<double> wRMS[NIFO_MAX];
107 
108 std::vector<wavearray<double> > vREC[NIFO_MAX]; // reconstructed
109 std::vector<double> vDREC[NIFO_MAX]; // detector rec waveform delays
110 
111 wavearray<int> sI; // shuffle index array
112 
114 TTree* net_tree = NULL;
116 bool detected=false;
117 
118 skymap rec_skyprob; // save reconstructed probability skymap
119 TH1D hist_skyprob; // used to extract random skyloc for skymask
120 
121 double inj_phi; // save injected phi
122 double inj_theta; // save injected theta
123 
124 double rec_phi; // save reconstructed phi
125 double rec_theta; // save reconstructed theta
126 
127 float rec_erA[11]; // save reconstructed error regions
128 
129 // ----------------------------------------------------
130 // Config Parameters
131 // ----------------------------------------------------
132 
133 //#define USE_wINJ // use wINJ (default wREC) (are pixels used as inj+noise in the retry)
134 //#define USE_wWHT // use wWHT (default wREC)
135 
136 #define SKYMASK_SEED 0
137 
138 #define APPLY_INJ_MRA
139 
140 
141 // user config options : default values or read from command line
148 int WRC_CED_RETRY_CONFIG=0; // dump CED at gLRETRY=WRC_CED_RETRY
151 
152 // config options used in the plugin
153 bool WRC_SKYMASK=false;
154 bool WRC_NOISE=false;
155 int WRC_RETRY=0;
156 bool WRC_AMP_CAL_ERR=false;
157 bool WRC_PHS_CAL_ERR=false;
159 int WRC_CED_RETRY=0; // dump CED at gLRETRY=WRC_CED_RETRY
160 int WRC_ID=0;
162 
163 #define AMP_CAL_ERR 0.1
164 #define PHS_CAL_ERR 10
165 
166 //#define COMPUTE_RESIDUAL_ENERGY
167 
168 // ----------------------------------------------------
169 // Output Parameters
170 // ----------------------------------------------------
171 
172 float erR[11]; // probability distribution of residuals
173 
174 #define nPAR 7
175 
176 TString parName[nPAR] = { "wrc_try",
177  "wrc_wf_noise",
178  "wrc_sm_phi",
179  "wrc_sm_theta",
180  "wrc_sm_radius",
181  "wrc_cal_amp",
182  "wrc_cal_phs"
183  };
184 float parValue[nPAR];
185 
186 float WRC_WF_NOISE=0;
189 float WRC_SM_RADIUS=0.1;
190 float WRC_CAL_AMP=0;
191 float WRC_CAL_PHS=0;
192 
193 #define WRC_PLUGIN_VERSION "v1.3"
194 
195 void
197 //!NOISE_MDC_SIMULATION
198 // This plugin is for Waveform Reconstruction Challenge (WRC) studies
199 // using the last Likelihood Stage Event Loop functionality
200 
201  cout << endl;
202  cout << "-----> CWB_Plugin_xWRC.C - " << WRC_PLUGIN_VERSION << endl;
203  cout << "ifo " << ifo.Data() << endl;
204  cout << "type " << type << endl;
205  cout << endl;
206 
207  if(type==CWB_PLUGIN_CONFIG) {
208 
209  cfg->outPlugin=true; // disable built-in output root file
210 
211  // set default config options
220 
221  // user config options : read from command line
222  ReadUserOptions();
223 
224  // print config options
225  cout << "-----------------------------------------" << endl;
226  cout << "WRC config options " << endl << endl;
227  cout << "WRC_RETRY " << WRC_RETRY << endl;
228  cout << "WRC_SKYMASK " << WRC_SKYMASK << endl;
229  cout << "WRC_NOISE " << WRC_NOISE << endl;
230  cout << "WRC_AMP_CAL_ERR " << WRC_AMP_CAL_ERR << endl;
231  cout << "WRC_PHS_CAL_ERR " << WRC_PHS_CAL_ERR << endl;
232  cout << "WRC_CED_ID " << WRC_CED_ID << endl;
233  cout << "WRC_CED_RETRY " << WRC_CED_RETRY << endl;
234  cout << "WRC_SKYMASK_REC_SKY " << WRC_SKYMASK_REC_SKY << endl;
235  cout << endl;
236  }
237 
238  if(type==CWB_PLUGIN_NETWORK) {
239  }
240 
241  if(type==CWB_PLUGIN_ILIKELIHOOD) {
242 
243  if(TString(cfg->analysis)!="2G") {
244  cout << "CWB_Plugin_xWRC.C -> CWB_PLUGIN_ILIKELIHOOD implemented only for 2G" << endl;
245  gSystem->Exit(1);
246  }
247 
248  if(NET->nLag!=1) {
249  cout << "CWB_Plugin_xWRC.C -> implemented only for zero lag" << endl;
250  gSystem->Exit(1);
251  }
252 
253  // export gLRETRY
254  EXPORT(int,gLRETRY,TString::Format("gLRETRY = %d;",WRC_RETRY).Data())
255  }
256 
257  if(type==CWB_PLUGIN_XLIKELIHOOD) {
258 
259  if(TString(cfg->analysis)!="2G") {
260  cout << "CWB_Plugin_xWRC.C -> CWB_PLUGIN_XLIKELIHOOD implemented only for 2G" << endl;
261  gSystem->Exit(1);
262  }
263 
264  // import gLRETRY
265  int gLRETRY=-1; IMPORT(int,gLRETRY)
266  cout << "-----> CWB_Plugin_xWRC.C -> " << " gLRETRY : " << gLRETRY << "/" << WRC_RETRY << " - WRC : ";
267 
268  if(WRC_NOISE) cout << "N ";
269  if(WRC_AMP_CAL_ERR) cout << "A ";
270  if(WRC_PHS_CAL_ERR) cout << "P ";
271  if(WRC_SKYMASK) cout << "S ";
272  cout << " - ";
273 #ifdef USE_wINJ // use wINJ (default wREC) (are pixels used as inj+noise in the retry)
274  cout << "wJ ";
275 #endif
276 #ifdef USE_wWHT // use wWHT (default wREC)
277  cout << "wW ";
278 #endif
279 #if !defined (USE_wINJ) && !defined (USE_wWHT)
280  cout << "wR ";
281 #endif
282  cout << endl;
283  if((WRC_NOISE||WRC_AMP_CAL_ERR||WRC_PHS_CAL_ERR) && gLRETRY<WRC_RETRY) {
284  if(rSS[0].size()>0) AddNoise2Sparse(NET,cfg, 0);
285  }
286  if(WRC_SKYMASK && gLRETRY<WRC_RETRY) SetSkyMask(NET, cfg, SKYMASK_SEED); // set new skymask
287 
288  if(gLRETRY==WRC_RETRY) { // first time
289  cfg->nSky=-999; NET->nSky=cfg->nSky;
290  } else {
291  cfg->nSky=1000; NET->nSky=cfg->nSky;
292  }
293 
294  if(WRC_CED_RETRY) {
295  if(gLRETRY==WRC_CED_RETRY) {
296  cfg->cedDump=true;
297  } else {
298  cfg->cedDump=false;
299  }
300  }
301 
302  // for each neew event the wrc config is initialized with the user input wrc values
303  if(gLRETRY==WRC_RETRY) { // first time
311 
312  // restore full skymask
313  sprintf(cfg->skyMaskFile,"--theta %f --phi %f --radius %f",0.,0.,360.);
314  //cout << cfg->skyMaskFile << endl;
315  xCWB.SetSkyMask(NET,cfg,cfg->skyMaskFile,'e');
316  }
317  }
318 
319  if(type==CWB_PLUGIN_OLIKELIHOOD) {
320 
321  if(TString(cfg->analysis)!="2G") {
322  cout << "CWB_Plugin_xWRC.C -> CWB_PLUGIN_OLIKELIHOOD implemented only for 2G" << endl;
323  gSystem->Exit(1);
324  }
325 
326  // import retry
327  int gLRETRY=-1; IMPORT(int,gLRETRY)
328 
329  // import ifactor
330  int gIFACTOR=-1; IMPORT(int,gIFACTOR)
331  double factor = cfg->simulation==3||cfg->simulation==4 ? -gIFACTOR : cfg->factors[gIFACTOR];
332  double ofactor=0;
333  if(cfg->simulation==4) ofactor=-factor;
334  else if(cfg->simulation==3) ofactor=-gIFACTOR;
335  else ofactor=factor;
336 
337  int nIFO = NET->ifoListSize(); // number of detectors
338  int rate = 0; // select all resolutions
339  int zlag = 0; // zelect only zero lag
340  netevent* EVT;
341 
342  SetOutput(NET, EVT, cfg); // set output root file
343 
344  wavearray<double>id = NET->getwc(zlag)->get(const_cast<char*>("ID"), 0, 'L', rate);
345 
346  for(int j=0; j<(int)id.size(); j++) { // loop over cluster index
347 
348  int ID = size_t(id.data[j]+0.5);
349 
350  if(NET->getwc(zlag)->sCuts[ID-1]!=-1) continue; // skip rejected/processed clusters
351 
352  if(WRC_ID>0 && ID!=WRC_ID) { // skip wrc analysis if ID!=WRC_ID
353  EXPORT(int,gLRETRY,"gLRETRY = 0;")
354  NET->getwc(zlag)->sCuts[ID-1]=1; // mark as processed
355  continue;
356  }
357 
358  if(WRC_CED_ID) { // enable/disable CED
359  if(ID==WRC_CED_ID) {
360  cfg->cedDump=true;
361  } else {
362  cfg->cedDump=false;
363  }
364  }
365 
366  if(gLRETRY==WRC_RETRY) { // mark event as detected
367  detected=true;
368  cout << "-----> Event Detected ID = " << ID << endl;
369  }
370 
371  GetRecWaveform(NET, EVT, ID, zlag, ofactor); // get & save reconstructed waveform
372 
373 #if defined (SAVE_PLOT_WERR) || defined (SAVE_PLOT_WERR2) || (SAVE_WHT_PLOT)
374  if(gLRETRY<WRC_RETRY) {
375  for(int n=0; n<nIFO; n++) {
376  wavearray<double>* wfINJ = &wINJ[n]; // injected waveform
377  wavearray<double>* wfREC = &vREC[n][vREC[n].size()-1]; // reconstructed waveform
378  ComputeErrorWF(wfINJ, wfREC, n);
379  }
380 #ifdef SAVE_WHT_PLOT
381  PlotWaveform(NET->ifoName[n], wfINJ, wfREC, cfg, false, false); // time
382  PlotWaveform(NET->ifoName[n], wfINJ, wfREC, cfg, true, false); // fft
383 #endif
384  }
385 
386  if(gLRETRY==0 && detected) { // last time
387  // final computation of the waveform errors
388  for(int n=0;n<nIFO;n++) {
389  for(int j=0;j<wAVR[n].size();j++) {
390  if(wTRY[n][j]==0) continue;
391  wAVR[n][j]/=wTRY[n][j];
392  wRMS[n][j]/=wTRY[n][j];
393  wRMS[n][j]-=wAVR[n][j]*wAVR[n][j];
394  wRMS[n][j] =sqrt(wRMS[n][j]);
395  //cout << n << " " << j << " wINJ : " << wINJ[n][j] << " wAVR : " << wAVR[n][j]
396  // << " wRMS : " << wRMS[n][j] << endl;
397  }
398 #ifdef SAVE_PLOT_WERR
399  PlotFinal(n);
400 #endif
401  }
402 #ifdef SAVE_PLOT_WERR2
403  PlotFinal2(NET, ID);
404 #endif
405  }
406 #endif
407 
408  if(gLRETRY==0 && detected) { // last time
409 #ifdef SAVE_MRA_PLOT // monster event display
410  PlotMRA(NET, ID, zlag, "last_wREC");
411 #endif
412  if(WRC_RETRY>0) ComputeErrorStatistic(NET, cfg, ID); // compute final waveform error statistic
413  DumpOutput(NET, EVT, cfg, ID, zlag, ofactor); // dump event to output root file
414  }
415 
416  // warning : must be done after DumpOutput (GetRInjWaveform destroy original pwc)
417  if(gLRETRY==WRC_RETRY) { // first time
418 #ifdef SAVE_MRA_PLOT // monster event display
419  PlotMRA(NET, ID, zlag, "wREC");
420 #endif
421 
422  GetInjWaveform(NET, EVT, ID, zlag, factor); // get & save injected waveform
423 
424  Wave2Sparse(NET,cfg,'v'); // save original sparse map to vSS
425  Wave2Sparse(NET,cfg,'r'); // rec
426  Wave2Sparse(NET,cfg,'j'); // inj
427  Wave2Sparse(NET,cfg,'s'); // wht
428  Wave2Sparse(NET,cfg,'n'); // noise
429  SaveSkyProb(NET,cfg,ID);
430 
431  GetRInjWaveform(NET, EVT, cfg, ID, zlag); // get injected waveform in wREC TF sub-domain
432  }
433 
434  // in the last loop the original event is reprocessed
435  // because must be saved on the final root file
436  if(gLRETRY==1 && detected) {
437  WRC_SKYMASK=true;
438  WRC_NOISE=false;
439  WRC_AMP_CAL_ERR=false;
440  WRC_PHS_CAL_ERR=false;
441  WRC_SKYMASK_REC_SKY=true;
442 
443  // restore original sparse maps
444  for(int n=0;n<nIFO;n++) {
445  detector* pD = NET->getifo(n);
446  pD->vSS=vSS[n];
447  }
448  }
449 
450  if(gLRETRY==0 && detected) { // last time
451 
452  // restore original sparse maps
453  for(int n=0;n<nIFO;n++) {
454  detector* pD = NET->getifo(n);
455  pD->vSS=vSS[n];
456  }
457 
458  // clean all temporary WRC structures
459  Clear();
460  detected=false;
461  }
462  }
463 
464  jfile->cd();
465 
466  if(EVT) delete EVT;
467  }
468 
469  return;
470 }
471 
473  double& enINJ, double& enREC, double& xcorINJ_REC) {
474 
475  double bINJ = wfINJ->start();
476  double eINJ = wfINJ->stop();
477  double bREC = wfREC->start();
478  double eREC = wfREC->stop();
479  //cout.precision(14);
480  //cout << "bINJ : " << bINJ << " eINJ : " << eINJ << endl;
481  //cout << "bREC : " << bREC << " eREC : " << eREC << endl;
482 
483  double R=wfINJ->rate();
484 
485  int oINJ = bINJ>bREC ? 0 : int((bREC-bINJ)*R+0.5);
486  int oREC = bINJ<bREC ? 0 : int((bINJ-bREC)*R+0.5);
487  //cout << "oINJ : " << oINJ << " oREC : " << oREC << endl;
488 
489  double startXCOR = bINJ>bREC ? bINJ : bREC;
490  double endXCOR = eINJ<eREC ? eINJ : eREC;
491  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
492  //cout << "startXCOR : " << startXCOR << " endXCOR : "
493  // << endXCOR << " sizeXCOR :" << sizeXCOR << endl;
494 
495  enINJ=0;
496  enREC=0;
497  xcorINJ_REC=0;
498  if (sizeXCOR<=0) return;
499 
500  // the enINJ, enREC, xcorINJ_REC are computed in the INJ range
501  for (int i=0;i<wfINJ->size();i++) enINJ+=wfINJ->data[i]*wfINJ->data[i];
502  for (int i=0;i<sizeXCOR;i++) enREC+=wfREC->data[i+oREC]*wfREC->data[i+oREC];
503  for (int i=0;i<sizeXCOR;i++) xcorINJ_REC+=wfINJ->data[i+oINJ]*wfREC->data[i+oREC];
504 
505  double erINJ_REC = enINJ+enREC-2*xcorINJ_REC;
506  //cout << "enINJ : " << enINJ << " enREC : " << enREC
507  // << " xcorINJ_REC : " << xcorINJ_REC << " enINJREC : " << enINJ+enREC-2*xcorINJ_REC << endl;
508  //cout << "erINJ_REC/enINJ : " << erINJ_REC/enINJ << endl;
509 
510  return;
511 }
512 
514 
515  double bINJ = wfINJ->start();
516  double eINJ = wfINJ->stop();
517  double bREC = wfREC->start();
518  double eREC = wfREC->stop();
519  //cout.precision(14);
520  //cout << "bINJ : " << bINJ << " eINJ : " << eINJ << endl;
521  //cout << "bREC : " << bREC << " eREC : " << eREC << endl;
522 
523  double R=wfINJ->rate();
524 
525  int oINJ = bINJ>bREC ? 0 : int((bREC-bINJ)*R+0.5);
526  int oREC = bINJ<bREC ? 0 : int((bINJ-bREC)*R+0.5);
527  //cout << "oINJ : " << oINJ << " oREC : " << oREC << endl;
528 
529  double startXCOR = bINJ>bREC ? bINJ : bREC;
530  double endXCOR = eINJ<eREC ? eINJ : eREC;
531  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
532  //cout << "startXCOR : " << startXCOR << " endXCOR : "
533  // << endXCOR << " sizeXCOR :" << sizeXCOR << endl;
534 
535  for (int i=0;i<sizeXCOR;i++) {
536  wTRY[ifoID][i+oINJ]+=1;
537  wAVR[ifoID][i+oINJ]+=wfREC->data[i+oREC];
538  wRMS[ifoID][i+oINJ]+=wfREC->data[i+oREC]*wfREC->data[i+oREC];
539  }
540 
541  return;
542 }
543 
544 void PlotMRA(network* NET, int ID, int k, TString tag) {
545 
546  if(NET->getwc(k)->sCuts[ID-1]!=-1) return; // no events
547 
548  // set sCuts=1 for the events which must be not copied with cps to pwc
549  std::vector<int> sCuts = NET->getwc(k)->sCuts; // save sCuts
550  for(int i=0; i<(int)sCuts.size(); i++) if(i!=ID-1) NET->getwc(k)->sCuts[i]=1;
551  // after cpf, pwc contains only one event
552  // ID can not be used to get the event, to get event use ID=1 (ex: for watplot)
553  netcluster* pwc = new netcluster();
554  pwc->cpf(*(NET->getwc(k))); // note: likelihood2G delete tdAmp
555 
556  int nIFO = NET->ifoListSize(); // number of detectors
557  watplot WTS(const_cast<char*>("wts"));
558  WTS.canvas->cd();
559  char fname[1024];
560  sprintf(fname, "%s/%s_l_tfmap_scalogram_%d.png",PLOT_DIR,tag.Data(),ID);
561  cout << "write " << fname << endl;
562  WTS.plot(pwc, 1, nIFO, 'L', 0, const_cast<char*>("COLZ"));
563  WTS.canvas->Print(fname);
564  WTS.clear();
565  sprintf(fname, "%s/%s_n_tfmap_scalogram_%d.png",PLOT_DIR,tag.Data(),ID);
566  cout << "write " << fname << endl;
567  WTS.plot(pwc, 1, nIFO, 'N', 0, const_cast<char*>("COLZ"));
568  WTS.canvas->Print(fname);
569  WTS.clear();
570 
571  if(pwc) delete pwc;
572  NET->getwc(k)->sCuts = sCuts; // restore sCuts
573 }
574 
576  CWB::config* cfg, bool fft, bool strain) {
577 
578  watplot PTS(const_cast<char*>("ptswrc"),200,20,800,500);
579 
580  //cout << "Print " << fname << endl;
581  double tmin = wfINJ->start()<wfREC->start() ? wfINJ->start() : wfREC->start();
582  wfINJ->start(wfINJ->start()-tmin);
583  wfREC->start(wfREC->start()-tmin);
584  if(fft) {
585  PTS.plot(wfINJ, const_cast<char*>("ALP"), 1, 0, 0, true, cfg->fLow, cfg->fHigh);
586  } else {
587  PTS.plot(wfINJ, const_cast<char*>("ALP"), 1, 0, 0);
588  }
589  PTS.graph[0]->SetLineWidth(1);
590  if(fft) {
591  PTS.plot(wfREC, const_cast<char*>("SAME"), 2, 0, 0, true, cfg->fLow, cfg->fHigh);
592  } else {
593  PTS.plot(wfREC, const_cast<char*>("SAME"), 2, 0, 0);
594  }
595  PTS.graph[1]->SetLineWidth(2);
596  wfINJ->start(wfINJ->start()+tmin);
597  wfREC->start(wfREC->start()+tmin);
598 
599  char label[64]="";
600  if(fft) sprintf(label,"%s","fft");
601  else sprintf(label,"%s","time");
602  if(strain) sprintf(label,"%s_%s",label,"strain");
603  else sprintf(label,"%s_%s",label,"white");
604 
605  char fname[1024];
606  sprintf(fname, "%s/%s_wf_%s_inj_rec_gps_%d.root",PLOT_DIR,ifo.Data(),label,int(tmin));
607  PTS.canvas->Print(fname);
608  cout << "write : " << fname << endl;
609  //PTS.canvas->Write(REPLACE(fname,dirCED,gtype));
610 }
611 
613 
614  // NEW SKYMAP
615 
616  gRandom->SetSeed(seed);
617  int isky = (int)hist_skyprob.GetRandom();
618  WRC_SM_THETA = 90-rec_skyprob.getTheta(isky);
619  WRC_SM_PHI = rec_skyprob.getPhi(isky);
620 
621  if(WRC_SKYMASK_REC_SKY) {
622  WRC_SM_THETA = 90-rec_theta;
624  }
625 
626 /*
627  // test code : use injected sky location for skymask
628  WRC_SM_THETA = 90-inj_theta;
629  WRC_SM_PHI = inj_phi;
630 */
631 
632  // --------------------------------------------------------
633  // define SetSkyMask
634  // --------------------------------------------------------
635  sprintf(cfg->skyMaskFile,"--theta %f --phi %f --radius %f",WRC_SM_THETA,WRC_SM_PHI,WRC_SM_RADIUS);
636  cout << cfg->skyMaskFile << endl;
637  xCWB.SetSkyMask(NET,cfg,cfg->skyMaskFile,'e');
638 
639 }
640 
641 void Wave2Sparse(network* NET, CWB::config* cfg, char wave_type) {
642 
643  if(wave_type!='j' && wave_type!='r' && wave_type!='s' && wave_type!='n' && wave_type!='v') {
644  cout << "Wave2Sparse - Error : wrong wave type" << endl; exit(1);
645  }
646 
647  gRandom->SetSeed(0);
648 
649  //cout << "-----------------> Wave2Sparse" << endl;
650 
651  // init waveform sparse maps
652  int nIFO = NET->ifoListSize(); // number of detectors
653  for(int n=0;n<nIFO;n++) {
654  detector* pD = NET->getifo(n);
655  if(wave_type=='v') vSS[n]=pD->vSS;
656  if(wave_type=='r') rSS[n]=pD->vSS;
657  if(wave_type=='j') jSS[n]=pD->vSS;
658  if(wave_type=='s') sSS[n]=pD->vSS;
659  }
660  if(wave_type=='v') return;
661 
662  int nRES = cfg->l_high-cfg->l_low+1; // number of frequency resolution levels
663 
664  // build waveform array
666  for(int n=0; n<nIFO; n++) {
667 
668  if(wave_type=='s') {
669  WF[n].resize(sSS[n][0].wdm_nSTS);
670  WF[n].rate(sSS[n][0].wdm_rate);
671  WF[n].start(sSS[n][0].wdm_start);
672  WF[n]=0;
673  int wos = double(sINJ[n].start()-WF[n].start())*WF[n].rate();
674  for (int i=0;i<sINJ[n].size();i++) WF[n][wos+i] = sINJ[n][i];
675  }
676  if(wave_type=='j') {
677  WF[n].resize(jSS[n][0].wdm_nSTS);
678  WF[n].rate(jSS[n][0].wdm_rate);
679  WF[n].start(jSS[n][0].wdm_start);
680  WF[n]=0;
681  int wos = double(wINJ[n].start()-WF[n].start())*WF[n].rate();
682  for (int i=0;i<wINJ[n].size();i++) WF[n][wos+i] = wINJ[n][i];
683  }
684  if(wave_type=='r') {
685  WF[n].resize(rSS[n][0].wdm_nSTS);
686  WF[n].rate(rSS[n][0].wdm_rate);
687  WF[n].start(rSS[n][0].wdm_start);
688  WF[n]=0;
689  int wos = double(wREC[n].start()-WF[n].start())*WF[n].rate();
690  for (int i=0;i<wREC[n].size();i++) WF[n][wos+i] = wREC[n][i];
691  }
692  if(wave_type=='n') { // noise
693  detector* pD = NET->getifo(n);
694  WF[n].resize(pD->vSS[0].wdm_nSTS);
695  WF[n].rate(pD->vSS[0].wdm_rate);
696  WF[n].start(pD->vSS[0].wdm_start);
697  WF[n]=0;
698  for (int i=0;i<WF[n].size();i++) WF[n][i] = gRandom->Gaus(0,1);
699  }
700 
701 #ifdef SAVE_WF_PLOT
702  gwavearray<double> gw(&WF[n]);
703  gw.Draw(GWAT_TIME);
704  watplot* plot = gw.GetWATPLOT();
705  TString gfile;
706  if(wave_type=='r') gfile=TString::Format("%s/REC_%s.png",PLOT_DIR,NET->ifoName[n]);
707  if(wave_type=='j') gfile=TString::Format("%s/INJ_%s.png",PLOT_DIR,NET->ifoName[n]);
708  if(wave_type=='s') gfile=TString::Format("%s/sINJ_%s.png",PLOT_DIR,NET->ifoName[n]);
709  if(wave_type!='n') (*plot) >> gfile;
710 #endif
711  }
712 
713  // fill waveform sparse maps
715  WDM<double>* pwdm = NULL;
716  double r00,r90;
717  for(int n=0; n<nIFO; n++) {
718  for(int i=0; i<nRES; i++) {
719  int level = i+cfg->l_low;
720  int layers = level>0 ? 1<<level : 0;
721  pwdm = NET->getwdm(layers+1);
722  w.Forward(WF[n],*pwdm);
723 
724  if((wave_type=='n')&&(i==nRES-1)) { // time shuffled index for noise array vN
725  sI.resize(w.sizeZero());
726  for(int j=0;j<sI.size();j++) sI[j]=j;
727  }
728 
729  int size;
730  if(wave_type=='r') size = rSS[n][i].sparseMap00.size();
731  if(wave_type=='j') size = jSS[n][i].sparseMap00.size();
732  if(wave_type=='s') size = sSS[n][i].sparseMap00.size();
733  if(wave_type=='n') vN[n].push_back(w);
734 
735  for(int j=0; j<size; j++) {
736  int index;
737  if(wave_type=='r') index = rSS[n][i].sparseIndex[j];
738  if(wave_type=='j') index = jSS[n][i].sparseIndex[j];
739  if(wave_type=='s') index = sSS[n][i].sparseIndex[j];
740  double v00 = w.pWavelet->pWWS[index];
741  double v90 = w.pWavelet->pWWS[index+w.maxIndex()+1];
742  if(wave_type=='r') {
743  rSS[n][i].sparseMap00[j]=v00;
744  rSS[n][i].sparseMap90[j]=v90;
745  }
746  if(wave_type=='j') {
747  jSS[n][i].sparseMap00[j]=v00;
748  jSS[n][i].sparseMap90[j]=v90;
749  }
750  if(wave_type=='s') {
751  sSS[n][i].sparseMap00[j]=v00;
752  sSS[n][i].sparseMap90[j]=v90;
753  }
754  }
755  }
756  }
757 }
758 
760 
761  //cout << "-----------------> AddNoise2Sparse" << endl;
762 
763  gRandom->SetSeed(seed);
764 
765  TF2 *gaus2d = new TF2("gaus2d","exp(-x*x/2)*exp(-y*y/2)", -5,5,-5,5);
766 
767  double d2r = TMath::Pi()/180.;
768  double C,S;
769 
770  // copy waveform sparse to detector sparse map
771  int nIFO = NET->ifoListSize(); // number of detectors
772  for(int n=0;n<nIFO;n++) {
773  detector* pD = NET->getifo(n);
774 #ifdef USE_wWHT // use wWHT (default wREC)
775  if(vSS[n].size()>0) pD->vSS=vSS[n];
776 #else
777  if(rSS[n].size()>0) pD->vSS=rSS[n];
778 #endif
779  }
780 
781  int nRES = cfg->l_high-cfg->l_low+1; // number of frequency resolution levels
782 
783  // random shuffle of vN[i] time indexes
784  std::srand(time(0));
785  random_shuffle(&(sI.data[0]),&(sI.data[sI.size()]));
786 
787  // add noise/mis-calibration to sparse map
788  double r00,r90;
789  for(int n=0; n<nIFO; n++) {
790  detector* pD = NET->getifo(n);
791  if(WRC_AMP_CAL_ERR) {
792  WRC_CAL_AMP = gRandom->Uniform(1-AMP_CAL_ERR,1+AMP_CAL_ERR);
793  cout << "WRC_CAL_AMP : " << WRC_CAL_AMP << endl;
794  }
795  if(WRC_PHS_CAL_ERR) {
796  WRC_CAL_PHS = gRandom->Uniform(1-PHS_CAL_ERR,1+PHS_CAL_ERR);
797  cout << "WRC_CAL_PHS : " << WRC_CAL_PHS << endl;
798  C = cos(-WRC_CAL_PHS*d2r);
799  S = sin(-WRC_CAL_PHS*d2r);
800  }
801  for(int i=0; i<nRES; i++) {
802  int size = pD->vSS[i].sparseMap00.size();
803  for(int j=0; j<size; j++) {
804  if(WRC_AMP_CAL_ERR) {
805  pD->vSS[i].sparseMap00[j]*=WRC_CAL_AMP;
806  pD->vSS[i].sparseMap90[j]*=WRC_CAL_AMP;
807  }
808  if(WRC_PHS_CAL_ERR) {
809  double aa = pD->vSS[i].sparseMap00[j];
810  double AA = pD->vSS[i].sparseMap90[j];
811 
812  pD->vSS[i].sparseMap00[j] = C*aa + S*AA;
813  pD->vSS[i].sparseMap90[j] = -S*aa + C*AA ;
814  }
815  if(WRC_NOISE) {
816  int index = pD->vSS[i].sparseIndex[j];
817 
818  int layers = vN[n][i].maxLayer()+1; // numbers of frequency bins (first & last bins have df/2)
819  int slices = vN[n][i].sizeZero(); // number of time bins
820  //cout << "layers " << layers << " slices " << slices << endl;
821 
822  int levels = (1<<(nRES-1-i));
823 
824  // index = itime*layers + ifreq;
825  int ifreq = index%layers;
826  int itime = (index-ifreq)/layers;
827  itime = levels*sI[itime/levels];
828  index = itime*layers + ifreq; // shuffled index
829 
830  r00 = vN[n][i].pWavelet->pWWS[index];
831  r90 = vN[n][i].pWavelet->pWWS[index+vN[n][i].maxIndex()+1];
832  //gaus2d->GetRandom2(r00,r90);
833 
834  pD->vSS[i].sparseMap00[j]+=r00;
835  pD->vSS[i].sparseMap90[j]+=r90;
836  }
837  }
838  }
839  }
840 }
841 
843 
844  std::vector<float>* vP;
845  std::vector<int>* vI;
846 
847  // save the probability skymap
848  if(NET->wc_List[0].p_Map.size()) {
849 
850  vP = &(NET->wc_List[0].p_Map[id-1]);
851  vI = &(NET->wc_List[0].p_Ind[id-1]);
852 
853  rec_skyprob = NET->getifo(0)->tau;
854  rec_skyprob = 0.;
855 
856  for(int j=0; j<int(vP->size()); j++) {
857  int i = (*vI)[j];
858  double th = rec_skyprob.getTheta(i);
859  double ph = rec_skyprob.getPhi(i);
860  int k=rec_skyprob.getSkyIndex(th, ph);
861  rec_skyprob.set(k,(*vP)[j]);
862  }
863  }
864 
865  hist_skyprob.SetBins(rec_skyprob.size(),0,rec_skyprob.size()-1);
866 
867  double prob=0;
868  for(int l=0;l<rec_skyprob.size();l++) prob+=rec_skyprob.get(l);
869  cout << "PROB : " << prob << endl;
870  for(int l=0;l<rec_skyprob.size();l++) hist_skyprob.SetBinContent(l,rec_skyprob.get(l));
871 
872 #ifdef SAVE_SKYPROB
873  // plot rec_skyprob
875  //gSM.SetOptions("hammer","Celestial",2);
876  gSM.SetZaxisTitle("probability");
877  gSM.Draw(1);
878  //TH2D* hsm = gSM.GetHistogram();
879  //hsm->GetZaxis()->SetTitleOffset(0.85);
880  //hsm->GetZaxis()->SetTitleSize(0.03);
881  gSM.Print(TString::Format("%s/rec_skyprob.png",PLOT_DIR));
882  exit(0);
883 #endif
884 }
885 
886 void PlotFinal(int n) {
887 
888  // draw envelope
889  gwavearray<double> geINJ=wINJ[n];
890  CWB::mdc::PhaseShift(geINJ,90);
891  for(int j=0;j<geINJ.size();j++) geINJ[j]=sqrt(wINJ[n][j]*wINJ[n][j]+geINJ[j]*geINJ[j]);
892  geINJ.Draw(GWAT_TIME);
893 
894  gwavearray<double> geAVR=wINJ[n];
895  for(int j=0;j<wAVR[n].size();j++) geAVR[j]=wAVR[n][j];
896  CWB::mdc::PhaseShift(geAVR,90);
897  for(int j=0;j<geAVR.size();j++) geAVR[j]=sqrt(wAVR[n][j]*wAVR[n][j]+geAVR[j]*geAVR[j]);
898  geINJ.Draw(&geAVR,GWAT_TIME,"SAME",kRed);
899 
900  gwavearray<double> geRMS=wINJ[n];
901  for(int j=0;j<wRMS[n].size();j++) geRMS[j]=wRMS[n][j];
902  geINJ.Draw(&geRMS,GWAT_TIME,"SAME",kBlue);
903 
904  watplot* plot = geINJ.GetWATPLOT();
905  TString gfile=TString::Format("%s/eAVR_%d.root",PLOT_DIR,n);
906  (*plot) >> gfile;
907 
908 
909  // draw frequency
910  double dt=1./wINJ[n].rate();
911  double Pi = TMath::Pi();
912  wavearray<double> xx=wINJ[n];
913  wavearray<double> XX=wINJ[n];
914  CWB::mdc::PhaseShift(XX,90);
915  gwavearray<double> gF=wINJ[n];
916  for(int j=1;j<gF.size();j++) {
917 //double angle = atan2(XX[j],xx[j])*180/Pi;
918 //angle += ceil( -angle / 360 ) * 360;
919 //gF[j]=angle;
920  gF[j]=atan2(XX[j],xx[j])-atan2(XX[j-1],xx[j-1]);
921  if(gF[j]<-Pi) gF[j]+=2*Pi;
922  if(gF[j]> Pi) gF[j]-=2*Pi;
923  gF[j]/=2*Pi*dt;
924  }
925  gF.Draw(GWAT_TIME);
926 
927  wavearray<double> yy=wINJ[n];
928  for(int j=0;j<yy.size();j++) yy[j]=wAVR[n][j];
930  CWB::mdc::PhaseShift(YY,90);
931  gwavearray<double> gFF=wINJ[n];
932  for(int j=1;j<gFF.size();j++) {
933  gFF[j]=atan2(YY[j],yy[j])-atan2(YY[j-1],yy[j-1]);
934  if(gFF[j]<-Pi) gFF[j]+=2*Pi;
935  if(gFF[j]> Pi) gFF[j]-=2*Pi;
936  gFF[j]/=2*Pi*dt;
937  }
938  gF.Draw(&gFF,GWAT_TIME,"SAME",kRed);
939 
940  plot = gF.GetWATPLOT();
941  gfile=TString::Format("%s/FREQ_%d.root",PLOT_DIR,n);
942  (*plot) >> gfile;
943 
944 
945  // draw envelope difference
946  gwavearray<double> geDIF=wINJ[n];
947  for(int j=0;j<geDIF.size();j++) geDIF[j]=geAVR[j]-geINJ[j];
948  geDIF.Draw(GWAT_TIME);
949  geDIF.Draw(&geRMS,GWAT_TIME,"SAME",kBlue);
950 
951  plot = geDIF.GetWATPLOT();
952  gfile=TString::Format("%s/eDIF_%d.root",PLOT_DIR,n);
953  (*plot) >> gfile;
954 
955 
956  // draw waveforms INJ vs AVR
957  gwavearray<double> gwAVR=wINJ[n];
958  gwAVR.Draw(GWAT_TIME);
959  for(int j=0;j<wAVR[n].size();j++) gwAVR[j]=wAVR[n][j];
960  gwAVR.Draw(GWAT_TIME,"SAME",kRed);
961  for(int j=0;j<wAVR[n].size();j++) gwAVR[j]=wRMS[n][j];
962  gwAVR.Draw(GWAT_TIME,"SAME",kBlue);
963 
964  plot = gwAVR.GetWATPLOT();
965  gfile=TString::Format("%s/wAVR_%d.root",PLOT_DIR,n);
966  (*plot) >> gfile;
967 
968  // draw waveforms REC vs AVR/INJ
969 
970  double bINJ = wINJ[n].start();
971  double eINJ = wINJ[n].stop();
972  double bREC = wREC[n].start();
973  double eREC = wREC[n].stop();
974  cout.precision(14);
975  cout << "bINJ : " << bINJ << " eINJ : " << eINJ << endl;
976  cout << "bREC : " << bREC << " eREC : " << eREC << endl;
977 
978  double R=wINJ[n].rate();
979 
980  int oINJ = bINJ>bREC ? 0 : int((bREC-bINJ)*R+0.5);
981  int oREC = bINJ<bREC ? 0 : int((bINJ-bREC)*R+0.5);
982  cout << "oINJ : " << oINJ << " oREC : " << oREC << endl;
983 
984  double startXCOR = bINJ>bREC ? bINJ : bREC;
985  double endXCOR = eINJ<eREC ? eINJ : eREC;
986  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
987  cout << "sizeXCOR : " << sizeXCOR << endl;
988  cout << "wINJ[n].size() : " << wINJ[n].size() << endl;
989  cout << "wREC[n].size() : " << wREC[n].size() << endl;
990 
991  gwavearray<double> gwREC = wINJ[n];
992  for(int j=0;j<sizeXCOR;j++) gwREC[j+oINJ]=wREC[n][j+oREC];
993  gwREC.Draw(GWAT_TIME);
994  //for(int j=0;j<wAVR[n].size();j++) gwREC[j]=wAVR[n][j];
995  for(int j=0;j<wINJ[n].size();j++) gwREC[j]=wINJ[n][j];
996  gwREC.Draw(GWAT_TIME,"SAME",kRed);
997  for(int j=0;j<wAVR[n].size();j++) gwREC[j]=wRMS[n][j];
998  gwREC.Draw(GWAT_TIME,"SAME",kBlue);
999 
1000  plot = gwREC.GetWATPLOT();
1001  gfile=TString::Format("%s/wREC_%d.root",PLOT_DIR,n);
1002  (*plot) >> gfile;
1003 }
1004 
1005 double GetSparseMap(SSeries<double>* SS, bool phase, int index) {
1006 
1007  int layer = SS->GetLayer(index);
1008  int start = SS->sparseLookup[layer]; // sparse table layer offset
1009  int end = SS->sparseLookup[layer+1]-1; // sparse table layer+1 offset
1010 
1011  int i = SS->binarySearch(SS->sparseIndex.data,start,end,index);
1012 
1013  if(i<0) return 0;
1014 
1015  return phase ? SS->sparseMap00[i] : SS->sparseMap90[i];
1016 }
1017 
1018 void PlotFinal2(network* NET, int ID) {
1019 
1020  double tMin,tMax;
1021  double start; // use ifoID=0 as offset
1022 
1023  int nIFO = NET->ifoListSize(); // number of detectors
1024 
1025  for(int n=0;n<nIFO;n++) {
1026 
1027  // draw waveforms INJ vs AVR
1028  gwavearray<double> gw=wINJ[n];
1029  //gwavearray<double> gw=rINJ[n];
1030 
1031  if(n==0) {
1032  start = gw.start();
1033  gw.GetTimeRange(tMin, tMax, 0.998);
1034  cout << "tMin " << tMin << " tMax " << tMax << endl;
1035  }
1036  gw.SetTimeRange(start+tMin, start+tMax);
1037  gw.Draw(GWAT_TIME, "CUSTOM");
1038  for(int j=0;j<wAVR[n].size();j++) gw[j]=wAVR[n][j];
1039  gw.Draw(GWAT_TIME,"SAME",kRed);
1040  for(int j=0;j<wAVR[n].size();j++) gw[j]=wRMS[n][j];
1041  gw.Draw(GWAT_TIME,"SAME",kBlue);
1042 
1043  watplot* plot = gw.GetWATPLOT();
1044 
1045  char gtitle[256];
1046  sprintf(gtitle,"%s : %10.3f (gps) - injected(Black) - average(Red) - rms(Blue)",gw.start(),NET->ifoName[n]);
1047  plot->gtitle(gtitle,"time(sec)","amplitude");
1048 
1049  TString gfile=TString::Format("%s/inj_vs_avr_%s_%d.root",PLOT_DIR,NET->ifoName[n],ID);
1050  //TString gfile=TString::Format("%s/rinj_vs_avr_%s_%d.root",PLOT_DIR,NET->ifoName[n],ID);
1051  (*plot) >> gfile;
1052 
1053  }
1054 }
1055 
1057 
1058  double tMin,tMax;
1059  double start; // use ifoID=0 as offset
1060 
1061  // draw waveforms INJ vs AVR
1062  gwavearray<double> gw=*w1;
1063 
1064  start = gw.start();
1065  gw.GetTimeRange(tMin, tMax, 0.998);
1066  cout << "tMin " << tMin << " tMax " << tMax << endl;
1067 
1068  gw.SetTimeRange(start+tMin, start+tMax);
1069  gw.Draw(GWAT_TIME, "CUSTOM");
1070  gw.Draw(w2,GWAT_TIME,"SAME",kRed);
1071 
1072  watplot* plot = gw.GetWATPLOT();
1073 
1074  char gtitle[256];
1075  sprintf(gtitle,"%s : %10.3f (gps) - %s",gw.start(),NET->ifoName[ifoID],tag.Data());
1076  plot->gtitle(gtitle,"time(sec)","amplitude");
1077 
1078  TString gfile=TString::Format("%s/%s_%s_%d.%s",PLOT_DIR,tag.Data(),NET->ifoName[ifoID],ID,gtype.Data());
1079  (*plot) >> gfile;
1080 
1081 }
1082 
1084 
1085  int nIFO = NET->ifoListSize(); // number of detectors
1086 
1087  double enINJ[NIFO_MAX], enREC[NIFO_MAX], xcorINJ_REC[NIFO_MAX];
1088 
1089  // shift all data to the first ifo arrival time
1090  for(int n=0;n<nIFO;n++) {
1091  CWB::mdc::TimeShift(wREC[n],-wDREC[n]);
1092  if(cfg->simulation) CWB::mdc::TimeShift(wINJ[n],-wDINJ[n]);
1093  if(cfg->simulation) CWB::mdc::TimeShift(rINJ[n],-wDREC[n]);
1094  for(int i=0;i<vREC[n].size();i++) {
1095  CWB::mdc::TimeShift(vREC[n][i],-vDREC[n][i]);
1096  }
1097  }
1098 
1099  double rnum=0;
1100  double rden=0;
1101 /*
1102  for(int n=0;n<nIFO;n++) {
1103  ComputeResidualEnergy(&wINJ[n], &wREC[n], enINJ[n], enREC[n], xcorINJ_REC[n]);
1104  //cout << n << " RESIDUAL : " << (enINJ[n]+enREC[n]-2*xcorINJ_REC[n])/enINJ[n] << endl;
1105  rnum+=(enINJ[n]+enREC[n]-2*xcorINJ_REC[n]);
1106  rden+=enINJ[n];
1107  }
1108  cout << "NET RESIDUAL : " << rnum/rden << endl;
1109 */
1110  // copy wINJ (only the wREC time range) into winj
1112  if(cfg->simulation) {
1113  for(int n=0;n<nIFO;n++) {
1114  winj[n]=wREC[n];
1115  winj[n]=0;
1116 
1117  double bINJ = wINJ[n].start();
1118  double eINJ = wINJ[n].stop();
1119  double bREC = wREC[n].start();
1120  double eREC = wREC[n].stop();
1121 
1122  double R=wINJ[0].rate();
1123 
1124  int oINJ = bINJ>bREC ? 0 : int((bREC-bINJ)*R+0.5);
1125  int oREC = bINJ<bREC ? 0 : int((bINJ-bREC)*R+0.5);
1126 
1127  double startXCOR = bINJ>bREC ? bINJ : bREC;
1128  double endXCOR = eINJ<eREC ? eINJ : eREC;
1129  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
1130 
1131  if (sizeXCOR<=0) continue;
1132 
1133  for (int j=0;j<sizeXCOR;j++) winj[n][j+oREC]+=wINJ[n].data[j+oINJ];
1134  }
1135  }
1136 
1137  // copy vREC (only the wREC time range) into wavr
1139  for(int n=0;n<nIFO;n++) {
1140  wavr[n]=wREC[n];
1141  wavr[n]=0;
1142  double bINJ = wREC[n].start();
1143  double eINJ = wREC[n].stop();
1144  for(int i=0;i<vREC[n].size();i++) {
1145 
1146  double bREC = vREC[n][i].start();
1147  double eREC = vREC[n][i].stop();
1148  //cout.precision(14);
1149  //cout << "bINJ : " << bINJ << " eINJ : " << eINJ << endl;
1150  //cout << "bREC : " << bREC << " eREC : " << eREC << endl;
1151 
1152  double R=wREC[n].rate();
1153 
1154  int oINJ = bINJ>bREC ? 0 : int((bREC-bINJ)*R+0.5);
1155  int oREC = bINJ<bREC ? 0 : int((bINJ-bREC)*R+0.5);
1156  //cout << "oINJ : " << oINJ << " oREC : " << oREC << endl;
1157 
1158  double startXCOR = bINJ>bREC ? bINJ : bREC;
1159  double endXCOR = eINJ<eREC ? eINJ : eREC;
1160  int sizeXCOR = int((endXCOR-startXCOR)*R+0.5);
1161  //cout << "startXCOR : " << startXCOR << " endXCOR : "
1162  // << endXCOR << " sizeXCOR :" << sizeXCOR << endl;
1163 
1164  if (sizeXCOR<=0) continue;
1165 
1166  for (int j=0;j<sizeXCOR;j++) wavr[n][j+oINJ]+=vREC[n][i].data[j+oREC];
1167  }
1168  wavr[n]*=1./vREC[n].size();
1169  }
1170 
1171  wavearray<float> vRES(vREC[0].size());
1172 
1173  for(int i=0;i<vREC[0].size();i++) {
1174  rnum=0;
1175  rden=0;
1176  for(int n=0;n<nIFO;n++) {
1177  ComputeResidualEnergy(&wavr[n], &vREC[n][i], enINJ[n], enREC[n], xcorINJ_REC[n]);
1178  //cout << i << " " << n << " RESIDUAL : " << (enINJ[n]+enREC[n]-2*xcorINJ_REC[n])/enINJ[n] << endl;
1179  rnum+=(enINJ[n]+enREC[n]-2*xcorINJ_REC[n]);
1180  rden+=enINJ[n];
1181  }
1182  //cout << "vRECvsAVR NET RESIDUAL : " << rnum/rden << endl;
1183  vRES[i]=rnum/rden;
1184  }
1185 
1186  double jres=0;
1187  if(cfg->simulation) {
1188  rnum=0;
1189  rden=0;
1190  for(int n=0;n<nIFO;n++) {
1191  ComputeResidualEnergy(&wavr[n], &winj[n], enINJ[n], enREC[n], xcorINJ_REC[n]);
1192  //cout << n << " RESIDUAL : " << (enINJ[n]+enREC[n]-2*xcorINJ_REC[n])/enINJ[n] << endl;
1193  rnum+=(enINJ[n]+enREC[n]-2*xcorINJ_REC[n]);
1194  rden+=enINJ[n];
1195  }
1196  //cout << "INJvsAVR NET RESIDUAL : " << rnum/rden << endl;
1197 
1198  rnum=0;
1199  rden=0;
1200  for(int n=0;n<nIFO;n++) {
1201  ComputeResidualEnergy(&wavr[n], &rINJ[n], enINJ[n], enREC[n], xcorINJ_REC[n]);
1202  //cout << n << " RESIDUAL : " << (enINJ[n]+enREC[n]-2*xcorINJ_REC[n])/enINJ[n] << endl;
1203  rnum+=(enINJ[n]+enREC[n]-2*xcorINJ_REC[n]);
1204  rden+=enINJ[n];
1205  }
1206  jres=rnum/rden;
1207  //cout << "rINJvsAVR NET RESIDUAL : " << rnum/rden << endl;
1208  }
1209 
1210  rnum=0;
1211  rden=0;
1212  for(int n=0;n<nIFO;n++) {
1213  ComputeResidualEnergy(&wavr[n], &wREC[n], enINJ[n], enREC[n], xcorINJ_REC[n]);
1214  //cout << n << " RESIDUAL : " << (enINJ[n]+enREC[n]-2*xcorINJ_REC[n])/enINJ[n] << endl;
1215  rnum+=(enINJ[n]+enREC[n]-2*xcorINJ_REC[n]);
1216  rden+=enINJ[n];
1217  }
1218  //cout << "RECvsAVR NET RESIDUAL : " << rnum/rden << endl;
1219 
1220 #ifdef SAVE_WF_COMPARISON
1221  for(int n=0;n<nIFO;n++) {
1222 // PlotFinal3(NET, ID, n, &wavr[n], &wavr[n], "AVRvsAVR","root");
1223 // PlotFinal3(NET, ID, n, &winj[n], &winj[n], "INJvsINJ","root");
1224 // PlotFinal3(NET, ID, n, &wREC[n], &wREC[n], "RECvsREC","root");
1225 
1226  PlotFinal3(NET, ID, n, &rINJ[n], &winj[n], "rINJvsINJ","root");
1227  PlotFinal3(NET, ID, n, &rINJ[n], &wREC[n], "rINJvsREC","root");
1228  PlotFinal3(NET, ID, n, &rINJ[n], &wavr[n], "rINJvsAVR","root");
1229 
1230  PlotFinal3(NET, ID, n, &winj[n], &wavr[n], "INJvsAVR","root");
1231  PlotFinal3(NET, ID, n, &wREC[n], &wavr[n], "RECvsAVR","root");
1232  }
1233 #endif
1234 
1235  wavearray<int> index(vRES.size());
1236  TMath::Sort(int(vRES.size()),const_cast<float*>(vRES.data),index.data,false);
1237 
1238  erR[0] = jres;
1239 
1240  for(int i=1;i<10;i++) erR[i]=vRES[index[i*int(vRES.size()/10.)]];
1241 
1242  erR[10]=0;
1243  for(int i=0;i<vRES.size();i++) {
1244  if(vRES[i]<=jres) erR[10]+=1.;
1245  }
1246  erR[10]/=vRES.size();
1247 
1248  cout.precision(3);
1249  cout << "erR : ";
1250  for(int i=0;i<10;i++) cout << erR[i] << "/";
1251  cout << erR[10] << endl;
1252 
1253 #ifdef SAVE_PLOT_RESIDUALS
1254 
1255  // plot residuals
1256  TCanvas* canvas2= new TCanvas("canvas2", "canvas2", 200, 20, 600, 600);
1257  TH1F* hres = new TH1F("residuals","residuals",10,0,1);
1258  for(int i=0;i<vRES.size();i++) hres->Fill(vRES[i]);
1259 
1260  hres->Draw();
1261  canvas2->SetLogy();
1262 
1263  float ymax = hres->GetMaximum();
1264  TLine *line = new TLine(jres,0,jres,ymax);
1265  line->SetLineColor(kRed);
1266  line->Draw();
1267 
1268  gStyle->SetLineColor(kWhite);
1269  hres->SetTitle("Distribution of residuals");
1270  hres->SetStats(kTRUE);
1271  canvas2->Print(TString::Format("%s/residuals_%d.png",PLOT_DIR,ID));
1272 
1273  delete hres;
1274  delete canvas2;
1275 
1276 #endif
1277 }
1278 
1279 void Clear() {
1280 
1281  sI.resize(0);
1282  for(int n=0;n<NIFO_MAX;n++) {
1283 
1284  wINJ[n].resize(0);
1285  sINJ[n].resize(0);
1286  wREC[n].resize(0);
1287 
1288  while(!vSS[n].empty()) vSS[n].pop_back();
1289  vSS[n].clear(); std::vector<SSeries<double> >().swap(vSS[n]);
1290  while(!sSS[n].empty()) sSS[n].pop_back();
1291  sSS[n].clear(); std::vector<SSeries<double> >().swap(sSS[n]);
1292  while(!rSS[n].empty()) rSS[n].pop_back();
1293  rSS[n].clear(); std::vector<SSeries<double> >().swap(rSS[n]);
1294  while(!jSS[n].empty()) jSS[n].pop_back();
1295  jSS[n].clear(); std::vector<SSeries<double> >().swap(jSS[n]);
1296 
1297  while(!vN[n].empty()) vN[n].pop_back();
1298  vN[n].clear(); std::vector<WSeries<double> >().swap(vN[n]);
1299 
1300  while(!wTRY[n].empty()) wTRY[n].pop_back();
1301  wTRY[n].clear(); std::vector<int >().swap(wTRY[n]);
1302 
1303  while(!wAVR[n].empty()) wAVR[n].pop_back();
1304  wAVR[n].clear(); std::vector<double >().swap(wAVR[n]);
1305 
1306  while(!wRMS[n].empty()) wRMS[n].pop_back();
1307  wRMS[n].clear(); std::vector<double >().swap(wRMS[n]);
1308 
1309  while(!vDREC[n].empty()) vDREC[n].pop_back();
1310  vDREC[n].clear(); std::vector<double >().swap(vDREC[n]);
1311 
1312  while(!vREC[n].empty()) vREC[n].pop_back();
1313  vREC[n].clear(); std::vector<wavearray<double> >().swap(vREC[n]);
1314 
1315  }
1316 }
1317 
1319 
1320  TString cwb_inet_options=TString(gSystem->Getenv("CWB_INET_OPTIONS"));
1321  if(cwb_inet_options.CompareTo("")!=0) {
1322  TString inet_options = cwb_inet_options;
1323  cout << inet_options << endl;
1324  if(!inet_options.Contains("--")) { // parameters are used only by cwb_inet
1325 
1326  TObjArray* token = TString(inet_options).Tokenize(TString(' '));
1327  for(int j=0;j<token->GetEntries();j++){
1328 
1329  TObjString* tok = (TObjString*)token->At(j);
1330  TString stok = tok->GetString();
1331 
1332  if(stok.Contains("wrc_id=")) {
1333  TString wrc_id=stok;
1334  wrc_id.Remove(0,wrc_id.Last('=')+1);
1335  if(wrc_id.IsDigit()) WRC_ID=wrc_id.Atoi();
1337  cout << "WRC_ID : " << WRC_ID << endl;
1338  }
1339 
1340  if(stok.Contains("wrc_retry=")) {
1341  TString wrc_retry=stok;
1342  wrc_retry.Remove(0,wrc_retry.Last('=')+1);
1343  if(wrc_retry.IsDigit()) WRC_RETRY=wrc_retry.Atoi();
1345  cout << "WRC_RETRY : " << WRC_RETRY << endl;
1346  }
1347 
1348  if(stok.Contains("wrc_ced_id=")) {
1349  TString wrc_ced_id=stok;
1350  wrc_ced_id.Remove(0,wrc_ced_id.Last('=')+1);
1351  if(wrc_ced_id.IsDigit()) WRC_CED_ID=wrc_ced_id.Atoi();
1353  cout << "WRC_CED_ID : " << WRC_CED_ID << endl;
1354  }
1355 
1356  if(stok.Contains("wrc_ced_retry=")) {
1357  TString wrc_ced_retry=stok;
1358  wrc_ced_retry.Remove(0,wrc_ced_retry.Last('=')+1);
1359  if(wrc_ced_retry.IsDigit()) WRC_CED_RETRY=wrc_ced_retry.Atoi();
1361  cout << "WRC_CED_RETRY : " << WRC_CED_RETRY << endl;
1362  }
1363 
1364  if(stok.Contains("wrc_skymask=")) {
1365  TString wrc_skymask=stok;
1366  wrc_skymask.Remove(0,wrc_skymask.Last('=')+1);
1367  if(wrc_skymask=="true") WRC_SKYMASK=true;
1368  if(wrc_skymask=="false") WRC_SKYMASK=false;
1370  }
1371 
1372  if(stok.Contains("wrc_noise=")) {
1373  TString wrc_noise=stok;
1374  wrc_noise.Remove(0,wrc_noise.Last('=')+1);
1375  if(wrc_noise=="true") WRC_NOISE=true;
1376  if(wrc_noise=="false") WRC_NOISE=false;
1378  }
1379 
1380  }
1381  }
1382  }
1383 }
1384 
1385 int _sse_mra_ps(network* NET, float* amp, float* AMP, float Eo, int K) {
1386 // fast multi-resolution analysis inside sky loop
1387 // select max E pixel and either scale or skip it based on the value of residual
1388 // pointer to 00 phase amplitude of monster pixels
1389 // pointer to 90 phase amplitude of monster pixels
1390 // Eo - energy threshold
1391 // K - max number of principle components to extract
1392 // returns number of MRA pixels
1393 
1394  int j,n,mm,J;
1395  int k = 0;
1396  int m = 0;
1397  int f = NIFO/4;
1398  int V = (int)NET->rNRG.size();
1399  float* ee = NET->rNRG.data; // residual energy
1400  float* c = NULL;
1401  float E2 = Eo/2; // threshold
1402  float E;
1403  float mam[NIFO];
1404  float mAM[NIFO];
1405 
1406  __m128* _m00 = (__m128*) mam;
1407  __m128* _m90 = (__m128*) mAM;
1408  __m128* _amp = (__m128*) amp;
1409  __m128* _AMP = (__m128*) AMP;
1410  __m128* _a00 = (__m128*) NET->a_00.data;
1411  __m128* _a90 = (__m128*) NET->a_90.data;
1412 
1413  while(k<K){
1414 
1415  for(j=0; j<V; ++j) if(ee[j]>ee[m]) m=j; // find max pixel
1416  if(ee[m]<=Eo) break; mm = m*f;
1417 
1418  //cout<<k<<" "<<" V= "<<V<<" m="<<m<<" ee[m]="<<ee[m];
1419 
1420  float cc = 0.;
1421 
1422  _sse_zero_ps(_m00);
1423  _sse_zero_ps(_m90);
1424 
1425  J = NET->wdmMRA.size()/7;
1426  c = NET->wdmMRA.getXTalk(m); // c1*c2+c3*c4=c1*c3+c2*c4=0
1427  for(j=0; j<J; j++) {
1428  n = int(c[0]+0.1);
1429  if(ee[n]>Eo) {
1430  _sse_rotadd_ps(_a00+n*f,c[1],_a90+n*f,c[3],_m00); // construct 00 vector
1431  _sse_rotadd_ps(_a00+n*f,c[2],_a90+n*f,c[4],_m90); // construct 90 vector
1432  }
1433  if(ee[n]>0) cc += c[5];
1434  c += 7;
1435  }
1436  _sse_mul_ps(_m00,cc>1?1./cc:0.7);
1437  _sse_mul_ps(_m90,cc>1?1./cc:0.7);
1438  E = _sse_abs_ps(_m00,_m90); // get PC energy
1439 
1440  if(E > ee[m]) { // correct overtuning PC
1441  _sse_cpf_ps(mam,_a00+mm); // store a00 for max pixel
1442  _sse_cpf_ps(mAM,_a90+mm); // store a90 for max pixel
1443  }
1444  _sse_add_ps(_amp+mm,_m00); // update 00 PC
1445  _sse_add_ps(_AMP+mm,_m90); // update 90 PC
1446 
1447  J = NET->wdmMRA.size()/7;
1448  c = NET->wdmMRA.getXTalk(m); // c1*c2+c3*c4=c1*c3+c2*c4=0
1449  for(j=0; j<J; j++) {
1450  n = int(c[0]+0.1);
1451  if(E<E2 && n!=m) {c+=7; continue;}
1452  if(ee[n]>Eo) {
1453  ee[n] = _sse_rotsub_ps(_m00,c[1],_m90,c[2],_a00+n*f); // subtract PC from a00
1454  ee[n]+= _sse_rotsub_ps(_m00,c[3],_m90,c[4],_a90+n*f); // subtract PC from a90
1455  ee[n]+= 1.e-6;
1456  }
1457  c += 7;
1458  }
1459  //cout<<" "<<ee[m]<<" "<<k<<" "<<E<<" "<<EE<<" "<<endl;
1460  //cout<<" "<<m<<" "<<cc<<" "<<ee[m]<<" "<<E<<" "<<pp[m]<<endl;
1461  k++;
1462  }
1463  return k;
1464 }
1465 
1467 
1468  // import slagShift
1469  float* gSLAGSHIFT=NULL; IMPORT(float*,gSLAGSHIFT)
1470 
1471  int nIFO = NET->ifoListSize(); // number of detectors
1472 
1473  // search output root file in the system list
1474  TFile* froot = NULL;
1475  TList *files = (TList*)gROOT->GetListOfFiles();
1476  outDump="";
1477  if (files) {
1478  TIter next(files);
1479  TSystemFile *file;
1480  TString fname;
1481  bool check=false;
1482  while ((file=(TSystemFile*)next())) {
1483  fname = file->GetName();
1484  // set output root file as the current file
1485  if(fname.Contains("wave_")) {
1486  froot=(TFile*)file;froot->cd();
1487  outDump=fname;
1488  outDump.ReplaceAll(".root.tmp",".txt");
1489  //cout << "output file name : " << fname << endl;
1490  }
1491  }
1492  if(!froot) {
1493  cout << "CWB_Plugin_xWRC.C : Error - output root file not found" << endl;
1494  gSystem->Exit(1);
1495  }
1496  } else {
1497  cout << "CWB_Plugin_xWRC.C : Error - output root file not found" << endl;
1498  gSystem->Exit(1);
1499  }
1500 
1501  net_tree = (TTree *) froot->Get("waveburst");
1502  if(net_tree!=NULL) {
1503  EVT = new netevent(net_tree,nIFO);
1504  for(int j=0;j<nPAR;j++) net_tree->SetBranchAddress(parName[j],&parValue[j]);
1505  net_tree->SetBranchAddress("erR",erR);
1506  } else {
1507  EVT = new netevent(nIFO);
1508  net_tree = EVT->setTree();
1509  for(int j=0;j<nPAR;j++)
1510  net_tree->Branch(parName[j],&parValue[j],TString::Format("%s/F",parName[j].Data()));
1511  net_tree->Branch("erR",erR,TString::Format("erR[%i]/F",11));
1512  }
1513  EVT->setSLags(gSLAGSHIFT); // set slags into netevent
1514  EVT->Psave=cfg->Psave;
1515 }
1516 
1517 void DumpOutput(network* NET, netevent* &EVT, CWB::config* cfg, int ID, int k, int factor) {
1518 
1519  for(int j=0;j<11;j++) NET->getwc(k)->sArea[ID-1][j]=rec_erA[j]; // restore erA
1520  if(cfg->dump) EVT->dopen(outDump.Data(),const_cast<char*>("a"),false);
1521  EVT->output2G(net_tree,NET,ID,k,factor); // get reconstructed parameters
1522  if(cfg->dump) EVT->dclose();
1523  if(!cfg->cedDump) NET->getwc(k)->sCuts[ID-1]=1; // mark as processed
1524 }
1525 
1526 void GetRecWaveform(network* NET, netevent* EVT, int ID, int k, int factor) {
1527 
1528  // import retry
1529  int gLRETRY=-1; IMPORT(int,gLRETRY)
1530 
1531  int nIFO = NET->ifoListSize(); // number of detectors
1532 
1533  EVT->output2G(NULL,NET,ID,k,factor); // get reconstructed parameters
1534 
1535  cout << "rec_theta : " << EVT->theta[0] << " rec_phi : " << EVT->phi[0] << endl;
1536 
1537  wavearray<double>** pwfREC = new wavearray<double>*[nIFO];
1538  for(int n=0; n<nIFO; n++) {
1539 
1540  detector* pd = NET->getifo(n);
1541  int idSize = pd->RWFID.size();
1542  int wfIndex=-1;
1543  for(int mm=0; mm<idSize; mm++) if (pd->RWFID[mm]==ID) wfIndex=mm;
1544 
1545  if(wfIndex<0) {
1546  cout << "CWB_Plugin_xWRC.C : Error : Reconstructed waveform not saved !!! : ID -> "
1547  << ID << " : detector " << NET->ifoName[n] << endl;
1548  continue;
1549  }
1550  if(wfIndex>=0) pwfREC[n] = pd->RWFP[wfIndex];
1551  wavearray<double>* wfREC = pwfREC[n]; // array of reconstructed waveforms
1552 
1553  double delay = EVT->time[n]-EVT->time[0];
1554  //cout << "DEB DELAY : " << n << " " << delay << endl;
1555 
1556  // save tst rec waveforms
1557  if(gLRETRY<WRC_RETRY) {
1558  vREC[n].push_back(*wfREC);
1559  vDREC[n].push_back(delay);
1560  }
1561  // save inj/rec waveforms
1562  if(gLRETRY==WRC_RETRY) { // first time
1563  wREC[n]=*wfREC;
1564  wDREC[n]=delay;
1565 
1566  // save injected/reconstruced sky location
1567  rec_theta = EVT->theta[0]; // rec theta
1568  rec_phi = EVT->phi[0]; // rec phi
1569 
1570  for(int j=0;j<11;j++) rec_erA[j] = EVT->erA[j];
1571  }
1572 
1573  // save wrc parameters
1574  parValue[0]=WRC_RETRY-gLRETRY;
1575  if(WRC_NOISE) parValue[1]=1; else parValue[1]=0;
1576  parValue[2]=WRC_SM_PHI;
1579  parValue[5]=WRC_CAL_AMP;
1580  parValue[6]=WRC_CAL_PHS;
1581  }
1582  delete [] pwfREC;
1583 }
1584 
1585 void GetInjWaveform(network* NET, netevent* EVT, int ID, int k, int factor) {
1586 
1587  // import retry
1588  int gLRETRY=-1; IMPORT(int,gLRETRY)
1589 
1590  int nIFO = NET->ifoListSize(); // number of detectors
1591 
1592  double recTime = EVT->time[0]; // rec event time det=0
1593 
1594  injection INJ(nIFO);
1595  // get inj ID
1596  int M = NET->mdc__IDSize();
1597  double injTime = 1.e12;
1598  int injID = -1;
1599  for(int m=0; m<M; m++) {
1600  int mdcID = NET->getmdc__ID(m);
1601  double T = fabs(recTime - NET->getmdcTime(mdcID));
1602  if(T<injTime && INJ.fill_in(NET,mdcID)) {
1603  injTime = T;
1604  injID = mdcID;
1605  }
1606  }
1607 
1608  if(INJ.fill_in(NET,injID)) { // get inj parameters
1609 
1610  wavearray<double>** pwfINJ = new wavearray<double>*[nIFO];
1611 
1612  // extract whitened injected & reconstructed waveforms
1613  for(int n=0; n<nIFO; n++) {
1614 
1615  detector* pd = NET->getifo(n);
1616 
1617  pwfINJ[n] = INJ.pwf[n];
1618  if (pwfINJ[n]==NULL) {
1619  cout << "CWB_Plugin_xWRC.C : Error : Injected waveform not saved !!! : detector "
1620  << NET->ifoName[n] << endl;
1621  continue;
1622  }
1623 
1624  double rFactor = 1.;
1625  rFactor *= factor>0 ? factor : 1;
1626  wavearray<double>* wfINJ = pwfINJ[n]; // array of injected waveforms
1627  *wfINJ*=rFactor; // injected wf is multiplied by the custom factor
1628 
1629  // save inj waveforms
1630  wINJ[n]=*wfINJ;
1631 
1632  //double delay = INJ.time[n]-INJ.time[0];
1633  double delay = EVT->time[n]-EVT->time[0];
1634  sINJ[n]=*wfINJ; CWB::mdc::TimeShift(sINJ[n],-delay);
1635  wDINJ[n]=delay;
1636 
1637  wTRY[n].resize(wINJ[n].size());
1638  wAVR[n].resize(wINJ[n].size());
1639  wRMS[n].resize(wINJ[n].size());
1640  for(int j=0;j<wINJ[n].size();j++) {wTRY[n][j]=0;wAVR[n][j]=0;wRMS[n][j]=0;};
1641 
1642  // save injected/reconstruced sky location
1643  inj_theta = EVT->theta[1]; // inj theta
1644  inj_phi = EVT->phi[1]; // inj phi
1645 
1646  *wfINJ*=1./rFactor; // restore injected amplitude
1647  }
1648  delete [] pwfINJ;
1649  }
1650 }
1651 
1653 
1654  int nIFO = NET->ifoListSize(); // number of detectors
1655 
1656  netcluster* pwc = NET->getwc(k);
1657 
1658  std::vector<int>* vint = &(pwc->cList[ID-1]); // pixel list
1659  int V = vint->size(); // cluster size
1660  if(!V) return;
1661 
1662  int irate=0;
1663  if(cfg->optim) { // extract optimal rate
1664  vector_int* pv = pwc->cRate.size() ? &(pwc->cRate[ID-1]) : NULL;
1665  irate = pv!=NULL ? (*pv)[0] : 0;
1666  }
1667  bool isPCs = !(cfg->optim&&std::isupper(cfg->search)); // are Principal Components ?
1668 
1669  detector* pD = NET->getifo(0);
1670  double R = pD->TFmap.rate();
1671 
1672  // ------------------------------------------------------------------
1673  // Find max pixel parameters
1674  // ------------------------------------------------------------------
1675 
1676  int K=0;
1677  for(int n=0; n<V; n++) {
1678  netpixel* pix = pwc->getPixel(ID,n);
1679  if(pix->layers%2==0) {
1680  cout << "CWB_Plugin_eWRC.C - Error : is enabled only for WDM (2G)" << endl;
1681  exit(1);
1682  }
1683  for(int m=0; m<nIFO; m++) {
1684 #ifdef APPLY_INJ_MRA
1685  NET->a_00[n*NIFO+m]=0;
1686  NET->a_90[n*NIFO+m]=0;
1687 #endif
1688  pix->setdata(0,'S',m);
1689  pix->setdata(0,'P',m);
1690  }
1691 #ifdef APPLY_INJ_MRA
1692  NET->rNRG[n]=0;
1693 #endif
1694  if(!pix->core) continue; // select only the principal components pixels
1695  if((irate)&&(irate != int(pix->rate+0.5))) continue; // if irate!=0 skip rate!=irate
1696  for(int m=0; m<nIFO; m++) {
1697  int index = pix->data[m].index;
1698  int ires = int(TMath::Log2(R/pix->rate))-cfg->l_low;
1699  double aa = GetSparseMap(&sSS[m][ires], true, index);
1700  double AA = GetSparseMap(&sSS[m][ires], false , index);
1701  pix->setdata(aa,'S',m);
1702  pix->setdata(AA,'P',m);
1703 #ifdef APPLY_INJ_MRA
1704  NET->a_00[n*NIFO+m]=aa;
1705  NET->a_90[n*NIFO+m]=AA;
1706  NET->rNRG[n]+=aa*aa+AA*AA;
1707 #endif
1708  }
1709 K++;
1710  }
1711 
1712 #ifdef SAVE_MRA_PLOT // monster event display
1713  PlotMRA(NET, ID, k, "rINJ");
1714 #endif
1715 
1716 #ifdef APPLY_INJ_MRA
1717 
1718  int V4 = V + (V%4 ? 4 - V%4 : 0);
1719  int V44 = V4 + 4;
1720 
1721  wavearray<float> amp(NIFO*V44); amp=0;
1722  wavearray<float> AMP(NIFO*V44); AMP=0;
1723  float Eo = 1e-3; // must be > 1e-6
1724  //float Eo = 2*cfg->Acore*cfg->Acore*nIFO;
1725 
1726  _sse_mra_ps(NET, amp.data, AMP.data, Eo, K);
1727  for(int n=0; n<V; n++) {
1728  netpixel* pix = pwc->getPixel(ID,n);
1729  if(pix->layers%2==0) {
1730  cout << "CWB_Plugin_eWRC.C - Error : is enabled only for WDM (2G)" << endl;
1731  exit(1);
1732  }
1733  for(int m=0; m<nIFO; m++) {
1734  pix->setdata(0,'S',m);
1735  pix->setdata(0,'P',m);
1736  }
1737  if(!pix->core) continue; // select only the principal components pixels
1738  if((irate)&&(irate != int(pix->rate+0.5))) continue; // if irate!=0 skip rate!=irate
1739  for(int m=0; m<nIFO; m++) {
1740  double aa = amp[n*NIFO+m];
1741  double AA = AMP[n*NIFO+m];
1742  pix->setdata(aa,'S',m);
1743  pix->setdata(AA,'P',m);
1744  }
1745  }
1746 
1747 #ifdef SAVE_MRA_PLOT // monster event display
1748  PlotMRA(NET, ID, k, "mra_rINJ");
1749 #endif
1750 
1751 #endif
1752 
1753  //NET->getMRAwave(ID,0,'S',0,true);
1754  NET->getMRAwave(ID,0,'S',0,false);
1755 
1756  for(int n=0; n<nIFO; n++) {
1757  detector* pD = NET->getifo(n);
1758  rINJ[n]=pD->waveForm;
1759  rINJ[n].start(pwc->start+pD->waveForm.start());
1760  CWB::mdc::TimeShift(rINJ[n],wDINJ[n]);
1761  }
1762  return;
1763 }
int WRC_CED_RETRY
monster wdmMRA
list of pixel pointers for MRA
Definition: network.hh:651
wavearray< double > gwREC[NIFO_MAX]
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
gskymap * gSM
static float _sse_abs_ps(__m128 *_a)
Definition: watsse.hh:137
TString outDump
std::vector< int > vector_int
Definition: cluster.hh:35
double wDINJ[NIFO_MAX]
static const double C
Definition: GNGen.cc:28
void AddNoise2Sparse(network *NET, CWB::config *cfg, int seed=0)
#define NIFO
Definition: wat.hh:74
void PlotFinal(int n)
size_t nLag
Definition: network.hh:573
std::vector< int > wTRY[NIFO_MAX]
std::vector< vector_int > cRate
Definition: netcluster.hh:398
bool WRC_NOISE
#define PLOT_DIR
int WRC_ID
void gtitle(TString title="", TString xtitle="", TString ytitle="")
Definition: watplot.cc:1256
int slices
void ComputeErrorWF(wavearray< double > *wfINJ, wavearray< double > *wfREC, int ifoID)
bool optim
Definition: config.hh:122
double M
Definition: DrawEBHH.C:13
void PlotWaveform(TString ifo, wavearray< double > *wfINJ, wavearray< double > *wfREC, CWB::config *cfg, bool fft=false, bool strain=false)
std::vector< netcluster > wc_List
Definition: network.hh:610
void SaveSkyProb(network *NET, CWB::config *cfg, int id)
void setSLags(float *slag)
Definition: netevent.cc:426
double GetSparseMap(SSeries< double > *SS, bool phase, int index)
std::vector< WSeries< double > > vN[NIFO_MAX]
void PlotFinal2(network *NET, int ID)
std::vector< double > * getmdcTime()
Definition: network.hh:422
virtual void rate(double r)
Definition: wavearray.hh:141
int irate
float factor
void Clear()
#define SKYMASK_SEED
char skyMaskFile[1024]
Definition: config.hh:308
int n
Definition: cwb_net.C:28
static void _sse_zero_ps(__m128 *_p)
Definition: watsse.hh:44
float WRC_CAL_AMP
wavearray< float > a_90
buffer for cluster sky 00 amplitude
Definition: network.hh:653
float gFF[3]
TTree * setTree()
Definition: netevent.cc:434
TString("c")
int WRC_CED_ID_CONFIG
int ID
Definition: TestMDC.C:70
size_t maxIndex()
Definition: wseries.hh:149
static void PhaseShift(wavearray< double > &x, double pShift=0.)
Definition: mdc.cc:2955
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
std::vector< pixdata > data
Definition: netpixel.hh:122
double phase
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
bool WRC_SKYMASK_REC_SKY
double fLow
Definition: config.hh:140
std::vector< vector_float > sArea
Definition: netcluster.hh:401
netpixel pix(nifo)
netcluster * pwc
Definition: cwb_job_obj.C:38
wavearray< double > rINJ[NIFO_MAX]
TH2F * ph
void Wave2Sparse(network *NET, CWB::config *cfg, char wave_type)
int GetLayer(int index)
Definition: sseries.hh:122
std::vector< TGraph * > graph
Definition: watplot.hh:194
plot gtitle(gtitle,"frequency (Hz)","strain/#sqrt{Hz}")
bool cedDump
Definition: config.hh:297
wavearray< int > sI
static float _sse_rotsub_ps(__m128 *_u, float c, __m128 *_v, float s, __m128 *_a)
Definition: watsse.hh:399
int layers
STL namespace.
double getTheta(size_t i)
Definition: skymap.hh:224
void DumpOutput(network *NET, netevent *&EVT, CWB::config *cfg, int ID, int k, int factor)
void Draw(int dpaletteId=1, Option_t *option="colfz")
Definition: gskymap.cc:460
static void _sse_cpf_ps(float *a, __m128 *_p)
Definition: watsse.hh:307
Long_t size
size_t layers
Definition: netpixel.hh:112
int m
Definition: cwb_net.C:28
DataType_t * pWWS
Definition: WaveDWT.hh:141
std::vector< vector_int > cList
Definition: netcluster.hh:397
virtual void start(double s)
Definition: wavearray.hh:137
int j
Definition: cwb_net.C:28
bool getMRAwave(size_t ID, size_t lag, char atype='S', int mode=0, bool tof=false)
Definition: network.cc:3666
i drho i
skymap tau
Definition: detector.hh:346
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
wavearray< int > sparseIndex
Definition: sseries.hh:180
static void _sse_rotadd_ps(__m128 *_u, float c, __m128 *_v, float s, __m128 *_a)
Definition: watsse.hh:388
int WRC_ID_CONFIG
TString parName[nPAR]
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
Definition: ComputeSNR.C:75
bool outPlugin
Definition: config.hh:369
bool core
Definition: netpixel.hh:120
char ifo[NIFO_MAX][8]
void dopen(const char *fname, char *mode, bool header=true)
Definition: netevent.hh:312
double rec_theta
size_t ifoListSize()
Definition: network.hh:431
bool WRC_PHS_CAL_ERR
cwb xCWB
WDM< double > * getwdm(size_t M)
param: number of wdm layers
Definition: network.hh:448
int Psave
Definition: config.hh:193
void clear()
Definition: watplot.hh:58
wavearray< double > w
Definition: Test1.C:27
#define nIFO
size_t getSkyIndex(double th, double ph)
param: theta param: phi
Definition: skymap.cc:720
virtual size_t size() const
Definition: wavearray.hh:145
TCanvas * canvas
Definition: watplot.hh:192
gWSeries< double > gw(w)
x plot
bool WRC_SKYMASK
jfile
Definition: cwb_job_obj.C:43
void output2G(TTree *, network *, size_t, int, double)
Definition: netevent.cc:711
wavearray< double > xx
Definition: TestFrame1.C:11
double rec_phi
TTree * net_tree
void GetInjWaveform(network *NET, netevent *EVT, int ID, int k, int factor)
Int_t Psave
number of detectors
Definition: netevent.hh:66
int simulation
Definition: config.hh:199
int l_high
Definition: config.hh:156
int WRC_CED_RETRY_CONFIG
std::vector< double > wRMS[NIFO_MAX]
#define IMPORT(TYPE, VAR)
Definition: cwb.hh:69
char search
Definition: config.hh:121
watplot * WTS
Definition: ChirpMass.C:119
bool WRC_PHS_CAL_ERR_CONFIG
xtalk getXTalk(int nLay1, size_t indx1, int nLay2, size_t indx2)
param: numbers of layers identifying the resolution of the first pixel param: TF map index of the fir...
Definition: monster.cc:351
wavearray< int > sparseLookup
Definition: sseries.hh:178
i() int(T_cor *100))
network NET
Definition: cwb_dump_inj.C:30
TString label
Definition: MergeTrees.C:21
double Pi
std::vector< SSeries< double > > rSS[NIFO_MAX]
const int NIFO_MAX
Definition: wat.hh:22
std::vector< SSeries< double > > sSS[NIFO_MAX]
wavearray< double > wREC[NIFO_MAX]
void ComputeErrorStatistic(network *NET, CWB::config *cfg, int ID)
double wDREC[NIFO_MAX]
wavearray< float > a_00
wdm multi-resolution analysis
Definition: network.hh:652
bool WRC_SKYMASK_CONFIG
void PlotFinal3(network *NET, int ID, int ifoID, wavearray< double > *w1, wavearray< double > *w2, TString tag, TString gtype="png")
double start
Definition: netcluster.hh:379
int l
bool detected
TIter next(twave->GetListOfBranches())
char fname[1024]
void SetOutput(network *NET, netevent *&EVT, CWB::config *cfg)
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *NET, WSeries< double > *x, TString ifo, int type)
COHERENCE.
std::vector< int > RWFID
Definition: detector.hh:381
size_t mdc__IDSize()
Definition: network.hh:414
bool WRC_SKYMASK_REC_SKY_CONFIG
double inj_phi
int k
double GetTimeRange(double &tMin, double &tMax, double efraction=0.9999999)
Definition: gwavearray.cc:359
int _sse_mra_ps(network *NET, float *amp, float *AMP, float Eo, int K)
size_t cpf(const netcluster &, bool=false, int=0)
Definition: netcluster.cc:117
std::vector< double > wAVR[NIFO_MAX]
TObjArray * token
void SetTimeRange(double tMin, double tMax)
Definition: gwavearray.hh:75
#define nPAR
int size()
Definition: monster.hh:91
Definition: skymap.hh:63
long nSky
Definition: config.hh:299
skymap rec_skyprob
float WRC_SM_PHI
double e
float WRC_SM_RADIUS
wavearray< double > yy
Definition: TestFrame5.C:12
wavearray< double > sINJ[NIFO_MAX]
std::vector< SSeries< double > > vSS[NIFO_MAX]
int ifreq
char tag[256]
Definition: cwb_merge.C:92
TFile * froot
netpixel * getPixel(size_t n, size_t i)
Definition: netcluster.hh:413
double dt
float rec_erA[11]
int WRC_CED_ID
static void TimeShift(wavearray< double > &x, double tShift=0.)
Definition: mdc.cc:2903
int l_low
Definition: config.hh:155
#define EXPORT(TYPE, VAR, CMD)
Definition: config.cc:92
#define AMP_CAL_ERR
double getPhi(size_t i)
Definition: skymap.hh:164
WSeries< double > TFmap
Definition: detector.hh:354
void GetRecWaveform(network *NET, netevent *EVT, int ID, int k, int factor)
void ComputeResidualEnergy(wavearray< double > *wfINJ, wavearray< double > *wfREC, double &enINJ, double &enREC, double &xcorINJ_REC)
Float_t * phi
sqrt(h+*h+ + hx*hx)
Definition: netevent.hh:87
void PlotMRA(network *NET, int ID, int k, TString tag)
bool WRC_AMP_CAL_ERR
netevent EVT(itree, nifo)
void ReadUserOptions()
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
static void _sse_add_ps(__m128 *_a, __m128 *_b)
Definition: watsse.hh:248
std::vector< int > sCuts
Definition: netcluster.hh:392
WSeries< double > waveForm
Definition: detector.hh:355
void dclose()
Definition: netevent.hh:321
double fHigh
Definition: config.hh:141
double inj_theta
bool WRC_AMP_CAL_ERR_CONFIG
wavearray< int > index
virtual void stop(double s)
Definition: wavearray.hh:139
double fabs(const Complex &x)
Definition: numpy.cc:55
void SetSkyMask(network *NET, CWB::config *cfg, int seed=0)
wavearray< float > sparseMap00
Definition: sseries.hh:181
int SetSkyMask(network *net, CWB::config *cfg, char *options, char skycoord, double skyres=-1)
Definition: cwb.cc:2484
wavearray< float > sparseMap90
Definition: sseries.hh:182
int gIFACTOR
TH1D hist_skyprob
Float_t * erA
noise rms
Definition: netevent.hh:108
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Definition: wseries.cc:246
Meyer< double > S(1024, 2)
netcluster * getwc(size_t n)
param: delay index
Definition: network.hh:439
std::vector< double > vDREC[NIFO_MAX]
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
int WRC_RETRY_CONFIG
float erR[11]
TString gfile
void set(size_t i, double a)
param: sky index param: value to set
Definition: skymap.hh:122
#define PHS_CAL_ERR
static void _sse_mul_ps(__m128 *_a, float b)
Definition: watsse.hh:56
watplot * GetWATPLOT()
Definition: gwavearray.hh:88
Definition: cwb.hh:136
DataType_t * data
Definition: wavearray.hh:319
wavearray< float > rNRG
buffer for cluster sky 90 amplitudes
Definition: network.hh:654
long nSky
Definition: network.hh:574
const int nRES
Definition: revMonster.cc:7
double * itime
int WRC_RETRY
double factors[FACTORS_MAX]
Definition: config.hh:202
double get(size_t i)
param: sky index
Definition: skymap.cc:699
bool WRC_NOISE_CONFIG
WaveDWT< DataType_t > * pWavelet
Definition: wseries.hh:456
wavearray< double > wINJ[NIFO_MAX]
char line[1024]
Bool_t fill_in(network *, int, bool=true)
Definition: injection.cc:349
float rate
Definition: netpixel.hh:113
std::vector< SSeries< double > > vSS
Definition: detector.hh:352
float parValue[nPAR]
size_t size()
Definition: skymap.hh:136
size_t getmdc__ID(size_t n)
Definition: network.hh:426
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
void Print(TString pname)
Definition: gskymap.cc:1122
int check
size_t sizeZero()
Definition: wseries.hh:144
#define WRC_PLUGIN_VERSION
float WRC_SM_THETA
std::vector< wavearray< double > > vREC[NIFO_MAX]
void Draw(GWAT_DRAW type=GWAT_TIME, TString options="ALP", Color_t color=kBlack)
Definition: gwavearray.cc:89
detector ** pD
void GetRInjWaveform(network *NET, netevent *EVT, CWB::config *cfg, int ID, int k)
std::vector< double > vRES
float WRC_CAL_PHS
void SetZaxisTitle(TString zAxisTitle)
Definition: gskymap.hh:150
int binarySearch(int array[], int start, int end, int key)
Definition: sseries.cc:461
exit(0)
std::vector< SSeries< double > > jSS[NIFO_MAX]
bool setdata(double a, char type='R', size_t n=0)
Definition: netpixel.hh:58
float WRC_WF_NOISE
bool dump
Definition: config.hh:295