Logo coherent WaveBurst  
Library Reference Guide
Logo
cwb2G.cc
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Gabriele Vedovato, Sergey Klimenko
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 #include "cwb2G.hh"
20 #include "TKey.h"
21 #include "regression.hh"
22 #include "watplot.hh"
23 #include "sseries.hh"
24 #include "gwavearray.hh"
25 #include <iomanip>
26 
27 // minimun skymap resolution used for subNetCut
28 #define MIN_SKYRES_HEALPIX 4
29 #define MIN_SKYRES_ANGLE 3
30 
31 // regression parameters
32 #define REGRESSION_FILTER_LENGTH 8
33 #define REGRESSION_MATRIX_FRACTION 0.95
34 #define REGRESSION_SOLVE_EIGEN_THR 0.
35 #define REGRESSION_SOLVE_EIGEN_NUM 10
36 #define REGRESSION_SOLVE_REGULATOR 'h'
37 #define REGRESSION_APPLY_THR 0.8
38 
39 // select function parameter
40 #define SELECT_SUBRHO 5.0
41 #define SELECT_SUBNET 0.1
42 
43 // WDM default parameters
44 #define WDM_BETAORDER 6 // beta function order for Meyer
45 #define WDM_PRECISION 10 // wavelet precision
46 
47 #define EXIT(ERR) gSystem->Exit(ERR) // better exit handling for ROOT stuff
48 
49 ClassImp(cwb2G)
50 
51 cwb2G::~cwb2G() {
52 //
53 // Destructor
54 //
55 
56  for(int i=0; i<NRES_MAX; i++) {
57  if(pwdm[i]!=NULL) delete pwdm[i];
58  }
59 }
60 
61 void cwb2G::Init() {
62 //
63 // - Initialize the WDM<double> tranforms for each resolution level used in the analysis
64 // - Load the xTalk Catalog (network::setMRAcatalog) used for the Multi Resolution Analysis
65 //
66 
67  // export this (cwb2G) object (can be used by plugins)
68  char cmd[256];
69  sprintf(cmd,"gCWB2G = (void*)%p;",this); EXPORT(void*,gCWB2G,cmd);
70 
71  // check & set analysis stage
72  if(TString(cfg.analysis)!="2G")
73  {cout << "cwb2G::Init - Error : analysis must be 2G" << endl;EXIT(1);}
74  if(istage!=CWB_STAGE_FULL && iname=="")
75  {cout << "cwb2G::Init - Error : if initial stage!=CWB_STAGE_FULL "
76  << "then input file name must be a job file " << endl;EXIT(1);}
77 
78  nRES = cfg.l_high-cfg.l_low+1; // number of frequency resolution levels
79 
80  for(int i=0; i<NIFO_MAX; i++) hot[i]=NULL;
81  for(int i=0; i<NRES_MAX; i++) pwdm[i]=NULL;
82 
83  // loading MRA catalog
84  char MRAcatalog[1024];
85  sprintf(MRAcatalog,"%s/%s",cfg.filter_dir,cfg.wdmXTalk);
86  cout << "cwb2G::Init - Loading catalog of WDM cross-talk coefficients ... " << endl;
87  cout << MRAcatalog << endl;
88  TB.checkFile(MRAcatalog);
89  NET.setMRAcatalog(MRAcatalog);
90 
91  int BetaOrder = WDM_BETAORDER; // beta function order for Meyer
92  int precision = WDM_PRECISION; // wavelet precision
93  if(NET.wdmMRA.tag!=0) { // new catalog format : read BetaOrder,precision from catalog
94  BetaOrder=NET.wdmMRA.BetaOrder;
95  precision=NET.wdmMRA.precision;
96  }
97  // print catalog infos
98  char info[256];
99  sprintf(info,"-Tag:%d-BetaOrder:%d-precision:%d",NET.wdmMRA.tag,BetaOrder,precision);
100  PrintAnalysisInfo(CWB_STAGE_INIT,"cwb2G::Init",info,true,false);
101 
102  //cout << "cwb2G::Init - Create WDM ... " << endl;
103  for(int level=cfg.l_high; level>=cfg.l_low; level--) {
104  int layers = level>0 ? 1<<level : 0;
105  pwdm[cfg.l_high-level] = new WDM<double>(layers,layers,BetaOrder,precision);
106  // check if filter lenght is less than cwb scratch length
107  double wdmFLen = double(pwdm[cfg.l_high-level]->m_H)/rateANA; // sec
108  if(wdmFLen > cfg.segEdge+0.001) {
109  cout << endl;
110  cout << "cwb2G::Init : Error - filter length must be <= segEdge !!!" << endl;
111  cout << "filter length : " << wdmFLen << " sec" << endl;
112  cout << "cwb scratch : " << cfg.segEdge << " sec" << endl;
113  EXIT(1);
114  } else {
115  cout << "Filter length = " << wdmFLen << " (sec)" << endl;
116  }
117  // check if the length for time delay amplitudes is less than cwb scratch length
118  // the factor 1.5 is used to avoid to use pixels on the border which could be distorted
119  double rate = rateANA>>level;
120  if(cfg.segEdge<int(1.5*(cfg.TDSize/rate)+0.5)) {
121  cout << endl;
122  cout << "cwb2G::Init : Error - segEdge must be > "
123  << "1.5x the length for time delay amplitudes!!!" << endl;
124  cout << "TD length : " << cfg.TDSize/rate << " sec" << endl;
125  cout << "segEdge : " << cfg.segEdge << " sec" << endl;
126  cout << "Select segEdge > " << int(1.5*(cfg.TDSize/rate)+0.5) << endl << endl;
127  EXIT(1);
128  }
129  // add WDM to network
130  NET.add(pwdm[cfg.l_high-level]); // network vector must be filled starting from max resolution level
131  }
132 
133  // check if analysis layers are contained in the MRAcatalog
134  // level : is the decomposition level
135  // layes : are the number of layers along the frequency axis rateANA/(rateANA>>level)
136  int check_layers=0;
137  for(int level=cfg.l_high; level>=cfg.l_low; level--) {
138  int layers = level>0 ? 1<<level : 0;
139  for(int j=0;j<NET.wdmMRA.nRes;j++) if(layers==NET.wdmMRA.layers[j]) check_layers++;
140  }
141 
142  if(check_layers!=nRES) {
143  cout << "cwb2G::Init - Error : analysis layers do not match the MRA catalog" << endl;
144  cout << endl << "analysis layers : " << endl;
145  for(int level=cfg.l_high; level>=cfg.l_low; level--) {
146  int layers = level>0 ? 1<<level : 0;
147  cout << "level : " << level << " layers : " << layers << endl;
148  }
149  cout << endl << "MRA catalog layers : " << endl;
150  for(int i=0;i<NET.wdmMRA.nRes;i++)
151  cout << "layers : " << NET.wdmMRA.layers[i] << endl;
152  EXIT(1);
153  }
154  else {
155  cout << endl;
156  for(int level=cfg.l_high; level>=cfg.l_low; level--) {
157  int layers = level>0 ? 1<<level : 0;
158  int rate = rateANA>>level;
159  cout << "level : " << level << "\t rate(hz) : " << rate << "\t layers : " << layers
160  << "\t df(hz) : " << rateANA/2./double(1<<level)
161  << "\t dt(ms) : " << 1000./rate << endl;
162  }
163  cout << endl;
164  }
165 
166  // Check if lagStep compatible with WDM parity
167  // This condition is necessary to avoid mixing between odd
168  // and even pixels when circular buffer is used for lag shift
169  // The MRAcatalog distinguish odd and even pixels
170  int rate_min = rateANA>>cfg.l_high;
171  double dt_max = 1./rate_min;
172  if(fmod(rate_min,1.)) {
173  cout << "cwb2G::Init - Error : rate min=" << rate_min << "(Hz) is not integer" << endl << endl;
174  EXIT(1);
175  }
176  if(int(cfg.lagStep*rate_min+0.001)&1) {
177  cout << "cwb2G::Init - Error : lagStep=" << cfg.lagStep << "(sec)"
178  << " is not a multple of 2*max_time_resolution=" << 2*dt_max << "(sec)" << endl << endl;
179  EXIT(1);
180  }
181  if(int(cfg.segEdge*rate_min+0.001)&1) {
182  cout << "cwb2G::Init - Error : segEdge=" << cfg.segEdge << "(sec)"
183  << " is not a multple of 2*max_time_resolution=" << 2*dt_max << "(sec)" << endl << endl;
184  EXIT(1);
185  }
186  if(int(cfg.segMLS*rate_min+0.001)&1) {
187  cout << "cwb2G::Init - Error : segMLS=" << cfg.segMLS << "(sec)"
188  << " is not a multple of 2*max_time_resolution=" << 2*dt_max << "(sec)" << endl << endl;
189  EXIT(1);
190  }
191 
192  // time-delay filter rate
193  if(cfg.fResample>0) { // RESAMPLING
195  } else {
197  }
198 
199  return;
200 }
201 
202 double cwb2G::ReadData(double mdcShift, int ifactor) {
203 //
204 // Read Noise & MDC data from frame file or "On The Fly" from plugin
205 //
206 // Loop over detectors
207 // - Read noise from frames or "On The Fly" from plugin
208 // - Resampling data
209 // - Read injections from frames or "On The Fly" from plugin (config::simulation>0)
210 // - Resampling data
211 // - if(simulation==2) MDC are rescaled (detector::setsnr) with a fixed
212 // network SNR according to the config::factors
213 // - Store noise & MDC to job file
214 //
215 
216  // if 2G analysis istage>=CWB_STAGE_STRAIN then skip read data
217  if(istage>=CWB_STAGE_STRAIN) return 0.;
218 
219  PrintStageInfo(CWB_STAGE_STRAIN,"cwb2G::ReadData");
220 
221  // data are loaded from root file
222  if((istage!=CWB_STAGE_INIT)&&(iname!="")) return cwb::ReadData(iname);
223 
224  Meyer<double> B(1024); // set wavelet for resampling
225  wavearray<double> x,z,w; // temporary time series
226  std::vector<wavearray<double> > y; // temporary time series for snr mode
227  y.resize(nIFO);
228  WSeries<double> wM; // mdc WSeries
230 
231  int layers_high = 1<<cfg.l_high;
232  WDM<double> WDMwhite(layers_high,layers_high,6,10); // set whitening WDM
233 
234  // check if wdm filter lenght is less than cwb scratch
235  double wdmFLen = double(WDMwhite.m_H)/rateANA; // sec
236  if(wdmFLen > cfg.segEdge+0.001) {
237  cout << endl;
238  cout << "cwb2G::ReadData : Error - filter scratch must be <= cwb scratch!!!" << endl;
239  cout << "filter length : " << wdmFLen << " sec" << endl;
240  cout << "cwb scratch : " << cfg.segEdge << " sec" << endl;
241  EXIT(1);
242  } else {
243  cout << "WDM filter length for regression = " << wdmFLen << " (sec)" << endl;
244  }
245 
246  int xsize=0.;
247  double xrate=0.;
248  double xstart=0.;
249 
250  // data are loaded from frame files
251  jfile = new TFile(jname,"UPDATE");
252  if(jfile==NULL||!jfile->IsOpen())
253  {cout << "cwb2G::ReadData - Error opening root file : " << jname << endl;EXIT(1);}
254  TDirectory* cdstrain = (TDirectory*)jfile->Get("strain");
255  if(cdstrain==NULL) cdstrain = jfile->mkdir("strain");
256 
257  for(int i=0; i<nIFO; i++) {
258 
259  if(cfg.simulation==4) { // sim4 -> read ifo strain from job file
260  px = (wavearray<double>*)jfile->Get(TString("strain/")+ifo[i]);
261  }
262  if((cfg.simulation==4)&&(px!=NULL)) { // sim4 -> use strain from job file
263  x = *px; delete px;
264  } else {
265  if(cfg.dataPlugin) { // data are provided by the user plugin
266  x.rate(cfg.inRate); x.start(FRF[i].start); x.resize(int(x.rate()*(FRF[i].stop-FRF[i].start)));
267  } else { // data are read from frame
269  }
271  if(TMath::IsNaN(x.mean()))
272  {cout << "cwb2G::ReadData - Error : found NaN in strain data !!!" << endl;EXIT(1);}
273 
274  if(x.rate()!=cfg.inRate)
275  {cout << "cwb2G::ReadData - input rate from frame " << x.rate()
276  << " do not match the one defined in config : " << cfg.inRate << endl;EXIT(1);}
277 
278  x.start(x.start()+cfg.dataShift[i]); // dataShift
279  x.start(x.start()-cfg.segLen*(segID[i]-segID[0])); // SLAG
281  if(cfg.dcCal[i]>0.) x*=cfg.dcCal[i]; // DC correction
282  if(cfg.fResample>0) x.Resample(cfg.fResample); // RESAMPLING
283  x.Resample(x.rate()/(1<<cfg.levelR)); // resampling
284  x*=sqrt(1<<cfg.levelR); // rescaling
285 
286  if(cfg.simulation==2) w = x; // snr mode - save strain
287 
288  // save ifo data to temporary job file
289  cdstrain->cd();gwavearray<double> gx(x);gx.Write(ifo[i],TObject::kOverwrite);
290  }
291 
292  if(i==0) {xrate=x.rate();xstart=x.start();xsize=x.size();}
293 
294  fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(),x.size()/x.rate(),x.rate());
295  if(i>0 && xstart != x.start()) {
296  cout << "cwb2G::ReadData - Error : ifo noise data not synchronized" << endl;
297  cout << ifo[i] << " " << x.start() << " != " << ifo[0] << " " << xstart << endl;
298  EXIT(1);
299  }
300  if(i>0 && xrate != x.rate()) {
301  cout << "cwb2G::ReadData - Error : ifo noise data have different rates" << endl;
302  cout << ifo[i] << " " << x.rate() << " != " << ifo[0] << " " << xrate << endl;
303  EXIT(1);
304  }
305 
306  if(cfg.simulation) {
307  TDirectory* cdmdc = (TDirectory*)jfile->Get("mdc");
308  if(cdmdc==NULL) cdmdc = jfile->mkdir("mdc");
309 
310  if(cfg.mdcPlugin) { // mdc are provided by the user plugin
311  x.rate(cfg.inRate); x.start(FRF[i+nIFO].start);
312  x.resize(int(x.rate()*(FRF[i+nIFO].stop-FRF[i+nIFO].start)));
313  } else { // mdc are read from frame
314  fr[nIFO].readFrames(FRF[i+nIFO],cfg.channelNamesMDC[i],x);
315  if(x.rate()!=cfg.inRate)
316  {cout << "cwb2G::ReadData - input rate from frame " << x.rate()
317  << " do not match the one defined in config : " << cfg.inRate << endl;EXIT(1);}
318  }
320  if(TMath::IsNaN(x.mean()))
321  {cout << "cwb2G::ReadData - Error : found NaN in MDC data !!!" << endl;EXIT(1);}
322 
323  x.start(x.start()+mdcShift);
324  if(cfg.fResample>0) x.Resample(cfg.fResample); // RESAMPLING
325  y[i] = x; // save inj for snr mode
326  x.Resample(x.rate()/(1<<cfg.levelR)); // resampling
327  x*=sqrt(1<<cfg.levelR); // rescaling
328 
329  if(cfg.simulation==2) { // snr mode
330 
331  // calculate noise rms
332  pTF[i] = pD[i]->getTFmap();
333  pTF[i]->Forward(w,WDMwhite);
334  pTF[i]->setlow(cfg.fLow);
335  pTF[i]->sethigh(cfg.fHigh);
337 
338  // compute snr
339  wM.Forward(x,WDMwhite);
340  wM.setlow(cfg.fLow);
341  wM.sethigh(cfg.fHigh);
342  // zero f<fLow to avoid whitening issues when psd noise is not well defined for f<fLow
343  int layers = wM.maxLayer();
344  for(int j=0;j<layers;j++) if(wM.frequency(j)<cfg.fLow) {
345  double layer = j+0.01; // -epsilon select 0 layer for 90 phase
346  wM.getLayer(z, layer);z=0;wM.putLayer(z, layer); // 0 phase
347  wM.getLayer(z,-layer);z=0;wM.putLayer(z,-layer); // 90 phase
348  }
349 
350  // compute mdc snr
351  pD[i]->setsim(wM,NET.getmdcTime(),cfg.iwindow/2.,cfg.segEdge,false);
352  }
353  else { // save mdc data to temporary job file
354  cdmdc->cd();gwavearray<double> gx(x);gx.Write(ifo[i],TObject::kOverwrite);
355  }
356 
357  fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(),x.size()/x.rate(),x.rate());
358  if(xstart != x.start()) {
359  cout << "cwb2G::ReadData - Error : mdc/noise data with different start time" << endl;
360  printf("start time : noise = %10.6f - mdc = %10.6f\n",xstart,x.start());
361  EXIT(1);
362  }
363  if(xrate != x.rate()) {
364  cout << "cwb2G::ReadData - Error : mdc/noise data with different rate" << endl;
365  printf("rate : noise = %10.6f - mdc = %10.6f\n",xrate,x.rate());
366  EXIT(1);
367  }
368  if(xsize != x.size()) {
369  cout << "cwb2G::ReadData - Error : mdc/noise data with different buffer size" << endl;
370  printf("buffer size : noise = %d - mdc = %lu\n",xsize,x.size());
371  EXIT(1);
372  }
373  }
374  if(singleDetector) break;
375  }
376 
377  // for snr mode the factors parameters set to be the mdc snr
378  if(cfg.simulation==2) {
379  TDirectory* cdmdc = (TDirectory*)jfile->Get("mdc");
380  if(cdmdc==NULL) cdmdc = jfile->mkdir("mdc");
381 
382  // compute rescale factor -> snr=1
383  std::vector<double> mdcFactor;
384  for (int k=0;k<(int)pD[0]->ISNR.size();k++) {
385  double snr=0;
386  if(singleDetector) {
387  for(int i=0; i<nIFO; i++) snr+=pD[0]->ISNR.data[k];
388  } else {
389  for(int i=0; i<nIFO; i++) snr+=pD[i]->ISNR.data[k];
390  }
391  snr=sqrt(snr);
392  if(snr>0) mdcFactor.push_back(1./snr); else mdcFactor.push_back(0.);
393  }
394 
395  size_t K = mdcFactor.size();
396  for(int k=0; k<(int)K; k++) {
397  if(mdcFactor[k]) cout << k << " mdcFactor : " << mdcFactor[k] << endl;
398  }
399 
400  // rescale mdc snr to 1
401  for(int i=0; i<nIFO; i++) {
402  pD[i]->setsnr(y[i],NET.getmdcTime(),&mdcFactor,cfg.iwindow/2.,cfg.segEdge);
404 
405  wM.Forward(y[i],B,cfg.levelR);
406  wM.getLayer(y[i],0);
407  cdmdc->cd();gwavearray<double> gyi(y[i]);gyi.Write(ifo[i],TObject::kOverwrite);
408 
409  // rescaled saved injected waveforms
410  for(int k=0; k<(int)pD[i]->IWFP.size(); k++) {
411  wavearray<double>* pwf = pD[i]->IWFP[k];
412  int iwfid = pD[i]->IWFID[k];
413  for(int j=0;j<(int)pwf->size();j++) pwf->data[j]*=mdcFactor[iwfid];
414  }
415  for(int k=0; k<(int)K; k++) {
416  pD[i]->HRSS[k]*=mdcFactor[k];
417  pD[i]->ISNR[k]*=mdcFactor[k];
418  }
419 
420  if(singleDetector) break;
421  }
422  // rescale amplitudes stored in the mdcList
423  for(int k=0; k<(int)K; k++) {
424  int ilog[5] = {1,3,12,13,14};
425  for(int l=0;l<5;l++) {
426  double mfactor = l<2 ? mdcFactor[k] : mdcFactor[k]*mdcFactor[k];
427  TString slog = TB.GetMDCLog(NET.mdcList[k], ilog[l]);
428  NET.mdcList[k]=TB.SetMDCLog(NET.mdcList[k], ilog[l], mfactor*slog.Atof());
429  }
430  }
431  }
432 
433  // mdc data are copied to detectors HoT to be accessible by the user plugin
434  for(int i=0; i<nIFO; i++) pD[i]->HoT = y[i];
436  // clean mdc from detectors HoT
437  for(int i=0; i<nIFO; i++) {pD[i]->HoT.resize(0);y[i].resize(0);}
438 
439  jfile->Close();
440 
441  x.resize(0);
442  z.resize(0);
443  w.resize(0);
444  wM.resize(0);
445 
446  return x.rate();
447 }
448 
449 void
451 //
452 // Apply regression to remove lines & whiten data
453 //
454 // Loop over detectors
455 // - read ifo strain from job file
456 // - read MDC data from temporary job file (config::simulation>0)
457 // - if(config::simulation==1) MDC are rescaled according to the config::factors
458 // - Add MDC to noise
459 // - Apply regression to remove lines
460 // - Use detector::white to estimate noise (detector::nRMS)
461 // - Use the estimated noise to whiten data (WSeries<double>::white)
462 // - Store injected waveforms (SaveWaveforms)
463 // - Store whitened data (detector::HoT) to job file (jfile)
464 // - Store estimated noise to job file (detector::nRMS)
465 //
466 
467  // data are loaded from root file
469  return DataConditioning(iname, ifactor);
470  // if 2G analysis istage==CWB_STAGE_SUPERCLUSTER then skip data conditioning
471  if(istage==CWB_STAGE_SUPERCLUSTER) return;
472 
473  PrintStageInfo(CWB_STAGE_CSTRAIN,"cwb2G::DataConditioning");
474 
475  char info[256];
480  TDirectory* cdrms=NULL;
481  TDirectory* cdcstrain=NULL;
482  double factor=cfg.factors[ifactor];
483 
484  int layers_high = 1<<cfg.l_high;
485  WDM<double> WDMwhite(layers_high,layers_high,6,10); // set whitening WDM
486  int layers = rateANA/8;
487  WDM<double> WDMlpr(layers,layers,6,10); // set LPE filter WDM
488 
489  // check if whitening WDM filter lenght is less than cwb scratch
490  double wdmFLen = double(WDMwhite.m_H)/rateANA; // sec
491  if(wdmFLen > cfg.segEdge+0.001) {
492  cout << endl;
493  cout << "cwb2G::DataConditioning : Error - filter scratch must be <= cwb scratch!!!" << endl;
494  cout << "filter length : " << wdmFLen << " sec" << endl;
495  cout << "cwb scratch : " << cfg.segEdge << " sec" << endl;
496  EXIT(1);
497  } else {
498  cout << "WDM filter max length = " << wdmFLen << " (sec)" << endl;
499  }
500 
501  // open input job file
502  TFile* ifile = ((istage==CWB_STAGE_STRAIN)&&(iname!="")) ? new TFile(iname) : NULL;
503  // open temporary job file
504  jfile = new TFile(jname, "UPDATE");
505  if(jfile==NULL||!jfile->IsOpen())
506  {cout << "cwb2G::DataConditioning - Error : file " << jname << " not found" << endl;EXIT(1);}
507 
508  // create cstrain,rms root dirs
510  TDirectory* dcstrain = (TDirectory*)jfile->Get("cstrain");
511  if(dcstrain==NULL) dcstrain=jfile->mkdir("cstrain");
512  char cdcstrain_name[32];sprintf(cdcstrain_name,"cstrain-f%d",ifactor);
513  cdcstrain=dcstrain->mkdir(cdcstrain_name);
514 
515  TDirectory* drms = (TDirectory*)jfile->Get("rms");
516  if(drms==NULL) drms=jfile->mkdir("rms");
517  char cdrms_name[32];sprintf(cdrms_name,"rms-f%d",ifactor);
518  cdrms=drms->mkdir(cdrms_name);
519  }
520 
521  for(int i=0; i<nIFO; i++) {
522  // read ifo strain from temporary job file
523  if(ifile!=NULL) px = (wavearray<double>*)ifile->Get(TString("strain/")+ifo[i]);
524  else px = (wavearray<double>*)jfile->Get(TString("strain/")+ifo[i]);
525  if(TMath::IsNaN(px->mean()))
526  {cout << "cwb2G::DataConditioning - Error : found NaN in strain data !!!" << endl;EXIT(1);}
527  hot[i] = pD[i]->getHoT();
528  *hot[i] = *px; delete px;
529 
530  if(cfg.simulation) {
531  // read mdc data from temporary job file
532  if(ifile!=NULL) px = (wavearray<double>*)ifile->Get(TString("mdc/")+ifo[i]);
533  else px = (wavearray<double>*)jfile->Get(TString("mdc/")+ifo[i]);
534  if(TMath::IsNaN(px->mean()))
535  {cout << "cwb2G::DataConditioning - Error : found NaN in MDC data !!!" << endl;EXIT(1);}
536 
537  if(cfg.simulation==3) { // time shift : factor is the shift time
538  int nshift = int(factor*px->rate()); // number of shifted samples
539  wavearray<double> pX = *px;(*px)=0; // temporary array
540  int jstart = nshift<0 ? -nshift : 0;
541  int jstop = nshift<0 ? px->size() : px->size()-nshift;
542  for(int j=jstart;j<jstop;j++) px->data[j+nshift] = pX.data[j];
543 
544  double tshift=nshift/px->rate(); // time shift (sec)
545  // take into account of the previous applied time shift
546  tshift = ifactor==0 ? tshift : tshift-int(cfg.factors[ifactor-1]*px->rate())/px->rate();
547 
548  // tshift saved injected waveforms
549  for(int k=0; k<(int)pD[i]->IWFP.size(); k++) {
550  wavearray<double>* pwf = pD[i]->IWFP[k];
551  pwf->start(pwf->start()+tshift);
552  }
553  // tshift saved central times
554  for(int k=0; k<(int)pD[i]->TIME.size(); k++) pD[i]->TIME[k]+=tshift;
555 
556  // shift times stored in the NET.mdcList & NET.mdcTime
557  if(i==0) {
558  vector<TString> ifos(nIFO);
559  for(int n=0;n<nIFO;n++) ifos[n]=ifo[n];
560  TB.shiftBurstMDCLog(NET.mdcList, ifos, tshift);
561  for(int k=0;k<(int)NET.mdcTime.size();k++) NET.mdcTime[k]+=tshift;
562  }
563  xM = *px; // copy MDC to temporary wavearray
564  } else if(cfg.simulation==4) {
565  xM = *px; // copy MDC to temporary wavearray
566  } else {
567  xM = *px; // copy MDC to temporary wavearray
568  (*px)*=factor;
569  }
570  hot[i]->add(*px);
571  delete px;
572  }
573 
574  if(bplugin)
576 
577  // 2G Data Conditioning & Whitening
578 
579  if(!cfg.dcPlugin) { // built in data conditioning
580  pTF[i] = pD[i]->getTFmap();
581 
582  // regression
583  pTF[i]->Forward(*hot[i],WDMlpr);
584  regression rr(*pTF[i],const_cast<char*>("target"),1.,cfg.fHigh);
585  rr.add(*hot[i],const_cast<char*>("target"));
590  *hot[i] = rr.getClean();
591 
592  // whitening
593  pTF[i]->Forward(*hot[i],WDMwhite);
594  pTF[i]->setlow(cfg.fLow);
595  pTF[i]->sethigh(cfg.fHigh);
596  pD[i]->white(cfg.whiteWindow,0,cfg.segEdge,cfg.whiteStride); // calculate noise rms
597  pD[i]->nRMS.bandpass(16.,0.,1); // high pass filtering at 16Hz
598  pTF[i]->white(pD[i]->nRMS,1); // whiten 0 phase WSeries
599  pTF[i]->white(pD[i]->nRMS,-1); // whiten 90 phase WSeries
600  if(cfg.simulation) { // estimated whitened MDC parms
601  wM.Forward(xM,WDMwhite);
602  wM.setlow(cfg.fLow);
603  wM.sethigh(cfg.fHigh);
604  // zero f<fLow to avoid whitening issues when psd noise is not well defined for f<fLow
605  int layers = wM.maxLayer();
606  for(int j=0;j<layers;j++) if(wM.frequency(j)<cfg.fLow) {
607  double layer = j+0.01; // -epsilon select 0 layer for 90 phase
608  wM.getLayer(x, layer);x=0;wM.putLayer(x, layer); // 0 phase
609  wM.getLayer(x,-layer);x=0;wM.putLayer(x,-layer); // 90 phase
610  }
611  // compute mdc snr and save whiten waveforms
612  pD[i]->setsim(wM,NET.getmdcTime(),cfg.iwindow/2.,cfg.segEdge,true);
613  }
614  WSeries<double> wtmp = *pTF[i];
615  pTF[i]->Inverse();
616  wtmp.Inverse(-2);
617  *hot[i] = *pTF[i];
618  *hot[i] += wtmp;
619  *hot[i] *= 0.5;
620  // add infos to history
621  sprintf(info,"-IFO:%d-RMS:%g",i,hot[i]->rms());
622  PrintAnalysisInfo(CWB_STAGE_CSTRAIN,"cwb2G::DataConditioning",info,false);
623  } else { // data conditioning is provided by the user plugin
624  char cmd[128];
625  // export to CINT variables
626  sprintf(cmd,"gRATEANA = %lu;",rateANA); EXPORT(size_t,gRATEANA,cmd);
627  sprintf(cmd,"gMDC = %p;",&xM); EXPORT(void*,gMDC,cmd);
628  if(bplugin)
630  }
631 
632  if(bplugin)
634 
636  SaveWaveforms(jfile, pD[i], ifactor); // write inj waveforms to output job file
637  }
638 
639  if(jobfOptions&CWB_JOBF_SAVE_CSTRAIN) {
640  cdcstrain->cd();gwavearray<double> ghi(*hot[i]);ghi.Write(ifo[i]);
641  cdrms->cd();pD[i]->nRMS.Write(ifo[i]);
642  }
643 
644  cout<<"After "<<ifo[i]<<" data conditioning"<<endl;
645  gSystem->Exec("/bin/date"); GetProcInfo();
646 
647  if(singleDetector) {
648  *pD[1]=*pD[0];
649  // copy detector data not implemented in the copy operator
650  pD[1]->HRSS = pD[0]->HRSS;
651  pD[1]->ISNR = pD[0]->ISNR;
652  pD[1]->FREQ = pD[0]->FREQ;
653  pD[1]->BAND = pD[0]->BAND;
654  pD[1]->TIME = pD[0]->TIME;
655  pD[1]->TDUR = pD[0]->TDUR;
656  pD[1]->IWFID = pD[0]->IWFID;
657  pD[1]->IWFP = pD[0]->IWFP;
658  pD[1]->RWFID = pD[0]->RWFID;
659  pD[1]->RWFP = pD[0]->RWFP;
660  break;
661  }
662  }
663  if(ifile!=NULL) ifile->Close();
664  jfile->Close();
665 
666  xM.resize(0); wM.resize(0); x.resize(0);
667 
668  // strains and mdc data are removed if not set in the jobfOptions (only for the last factor)
669  int ioffset = (cfg.simulation==4) ? int(cfg.factors[0]) : 0; // ifactor offset
670  if(ifactor==ioffset+cfg.nfactor-1) { // the last factor
671  vector<TString> delObjList;
672  if(!(jobfOptions&CWB_JOBF_SAVE_STRAIN)) delObjList.push_back("strain");
673  if(cfg.simulation && !(jobfOptions&CWB_JOBF_SAVE_MDC)) delObjList.push_back("mdc");
674  FileGarbageCollector(jname,"",delObjList);
675  }
676 
677  return;
678 }
679 
680 void
682 //
683 // Read Conditioned Data from Job File (from the previous processed stage)
684 //
685 // Loop over detectors
686 // - Read conditioned data
687 //
688 
689  PrintStageInfo(CWB_STAGE_CSTRAIN,"cwb2G::DataConditioning from file");
690 
691  char cdrms_name[32];sprintf(cdrms_name,"rms-f%d",ifactor);
692  char cdcstrain_name[32];sprintf(cdcstrain_name,"cstrain-f%d",ifactor);
693 
694  // data are loaded from input root file
695  TFile* ifile = new TFile(fName);
696  if(ifile==NULL)
697  {cout << "cwb2G::DataConditioning - Error opening root file : " << fName << endl;EXIT(1);}
698  // open temporary job file
699  jfile = new TFile(jname,"UPDATE");
700  if(jfile==NULL||!jfile->IsOpen())
701  {cout << "cwb::DataConditioning - Error opening root file : " << jname << endl;EXIT(1);}
702 
703  // open job name
704  TDirectory* jcdrms = NULL;
705  TDirectory* jcdcstrain = NULL;
707  TDirectory* dcstrain = (TDirectory*)jfile->Get("cstrain");
708  if(dcstrain==NULL) dcstrain=jfile->mkdir("cstrain");
709  jcdcstrain=dcstrain->mkdir(cdcstrain_name);
710 
711  TDirectory* drms = (TDirectory*)jfile->Get("rms");
712  if(drms==NULL) drms=jfile->mkdir("rms");
713  jcdrms=drms->mkdir(cdrms_name);
714  }
715 
716  WSeries<double>* pws=NULL;
717  for(int i=0; i<nIFO; i++) {
718  pTF[i] = pD[i]->getTFmap();
719  pTF[i]->setlow(cfg.fLow);
720  pTF[i]->sethigh(cfg.fHigh);
721  pTF[i]->w_mode=1;
722 
724  LoadWaveforms(ifile, pD[i], ifactor); // read inj waveforms from input job file
725  SaveWaveforms(jfile, pD[i], ifactor); // write inj waveforms to output job file
726  }
727 
728  // restore ctrain
729  ifile->cd();
730  pws = (WSeries<double>*)ifile->Get(TString("cstrain/")+cdcstrain_name+"/"+pD[i]->Name);
731  if(pws==NULL)
732  {cout << "cwb2G::DataConditioning - Error : cstrain not present, job terminated!!!" << endl;EXIT(1);}
733  if(jobfOptions&CWB_JOBF_SAVE_CSTRAIN)
734  {jfile->cd();jcdcstrain->cd();pws->Write(ifo[i]);}
735  hot[i] = pD[i]->getHoT();
736  *hot[i]=*pws;
737  *pTF[i]=*hot[i];
738  delete pws;
739 
740  // restore strain rms
741  ifile->cd();
742  pws = (WSeries<double>*)ifile->Get(TString("rms/")+cdrms_name+"/"+pD[i]->Name);
743  if(pws==NULL)
744  {cout << "cwb2G::DataConditioning - Error : strain rms not present, job terminated!!!" << endl;EXIT(1);}
745  if(jobfOptions&CWB_JOBF_SAVE_CSTRAIN)
746  {jfile->cd();jcdcstrain->cd();jcdrms->cd();pws->Write(ifo[i]);}
747  pD[i]->nRMS=*pws;;
748  delete pws;
749 
750  if(singleDetector) {
751  *pD[1]=*pD[0];
752  // copy detector data not implemented in the copy operator
753  pD[1]->HRSS = pD[0]->HRSS;
754  pD[1]->ISNR = pD[0]->ISNR;
755  pD[1]->FREQ = pD[0]->FREQ;
756  pD[1]->BAND = pD[0]->BAND;
757  pD[1]->TIME = pD[0]->TIME;
758  pD[1]->TDUR = pD[0]->TDUR;
759  pD[1]->IWFID = pD[0]->IWFID;
760  pD[1]->IWFP = pD[0]->IWFP;
761  pD[1]->RWFID = pD[0]->RWFID;
762  pD[1]->RWFP = pD[0]->RWFP;
763  break;
764  }
765  }
766 
767  if(jfile!=NULL) jfile->Close();
768  ifile->Close();
769  return;
770 }
771 
772 void
774 //
775 // Select the significant pixels
776 //
777 // Loop over resolution levels (nRES)
778 // - Loop over detectors (cwb::nIFO)
779 // - Compute the maximum energy of TF pixels (WSeries<double>::maxEnergy)
780 // - Set pixel energy selection threshold (network::THRESHOLD)
781 // - Loop over time lags (network::nLag)
782 // - Select the significant pixels (network::getNetworkPixels)
783 // - Single resolution clustering (network::cluster)
784 // - Store selected pixels to job file (netcluster::write)
785 //
786 
787  // if 2G analysis istage>=CWB_STAGE_COHERENCE then skip coherence
788  if(istage>=CWB_STAGE_COHERENCE) return;
789 
790  if(bplugin) {
791  char sfactor[8];sprintf(sfactor,"%d",ifactor);
792  CWB_Plugin(jfile,&cfg,&NET,NULL,sfactor,CWB_PLUGIN_BCOHERENCE);
793  }
794 
795  PrintStageInfo(CWB_STAGE_COHERENCE,"cwb2G::Coherence");
796 
797  char info[256];
798  double Eo;
799  netcluster* pwc;
800  netcluster wc;
801 
802  // used to build sparse map
803  SSeries<double> SS;
804  std::vector<SSeries<double> > vSS;
805  for(int n=0; n<nIFO; n++) vSS.push_back(SS);
807 
808  // open temporary job file
809  jfile = new TFile(jname,"UPDATE");
810  if(jfile==NULL||!jfile->IsOpen())
811  {cout << "cwb2G::Coherence - Error : file " << jname << " not found" << endl;EXIT(1);}
812 
813  // create sparse directories in the job root file to store sparse TF maps
814  TDirectory* sdir = (TDirectory*)jfile->Get("csparse");
815  if(sdir==NULL) sdir = jfile->mkdir("csparse");
816 
817  if(bplugin) {
818  char sfactor[8];sprintf(sfactor,"%d",ifactor);
819  CWB_Plugin(jfile,&cfg,&NET,NULL,sfactor,CWB_PLUGIN_ICOHERENCE);
820  }
821 
822  if(!cfg.cohPlugin) { // built in coherence stage
823  TStopwatch watchCoh; // coherence benchmark
824  int upN = rateANA/1024; if(upN<1) upN=1; // calculate upsample factor
825 
826  for(int i=0; i<nRES; i++) { // loop over TF resolutions
827  watchCoh.Start();
828  // print level infos
829  int level=cfg.l_high-i;
830  int layers = level>0 ? 1<<level : 0;
831  int rate = rateANA>>level;
832  cout << "level : " << level << "\t rate(hz) : " << rate
833  << "\t layers : " << layers << "\t df(hz) : " << rateANA/2./double(1<<level)
834  << "\t dt(ms) : " << 1000./rate << endl;
835 
836  // produce TF maps with max over the sky energy
837  double alp=0;
838  for(int n=0; n<nIFO; n++) {
839  alp+=NET.getifo(n)->getTFmap()->maxEnergy(*hot[n],*pwdm[i],mTau,upN,NET.pattern);
840  // restore the frequency boundaries changed by the maxEnergy call
841  NET.getifo(n)->getTFmap()->setlow(cfg.fLow);
843  if(singleDetector) {
844  *(NET.getifo(1)->getTFmap()) = *(NET.getifo(0)->getTFmap());
845  break;
846  }
847  }
848 
849  if(bplugin) { // call user plugin
850  char cmd[128];
851  // export resolution level to CINT
852  sprintf(cmd,"gILEVEL = %lu;",level); EXPORT(size_t,gILEVEL,cmd);
854  }
855 
856  // threshold on pixel energy
857  alp /= nIFO;
858  if(NET.pattern!=0) {
859  Eo = NET.THRESHOLD(cfg.bpp,alp);
860  } else {
861  Eo = NET.THRESHOLD(cfg.bpp);
862  }
863  cout.precision(5);
864  cout<<"thresholds in units of noise variance: Eo="<<Eo<<" Emax="<<Eo*2<<endl;
865  // add infos to history
866  sprintf(info,"-RES:%d-THR:%g",i,Eo);
867  PrintAnalysisInfo(CWB_STAGE_COHERENCE,"cwb2G::Coherence",info,false);
868 
869  double TL = NET.setVeto(cfg.iwindow);
870  cout<<"live time in zero lag: "<<TL<<endl; // set veto array
871  if(TL <= 0.) {froot->Close();EXIT(1);} // exit if live time is zero
872 
873  if(bplugin) { // call user plugin
874  char cmd[128];
875  // export resolution level to CINT
876  sprintf(cmd,"gILEVEL = %lu;",level); EXPORT(size_t,gILEVEL,cmd);
878  }
879 
880 
881  // init sparse table (used in supercluster stage : set the TD filter size)
882  pwdm[i]->setTDFilter(cfg.TDSize, 1);
883  for(int n=0; n<nIFO; n++) {
884  WS[n].Forward(*hot[n],*pwdm[i]);
885  vSS[n].SetMap(&WS[n]);
886  vSS[n].SetHalo(mTau);
887  if(singleDetector) {
888  WS[1]=WS[0];
889  vSS[1].SetMap(&WS[1]);
890  vSS[1].SetHalo(mTau);
891  break;
892  }
893  }
894 
895  // select pixels
896  if(cfg.simulation) {cout<<"ifactor|clusters|pixels ";cout.flush();}
897  else {cout<<"lag|clusters|pixels "; cout.flush();}
898  int csize_tot=0;int psize_tot=0;
899  for(int j=0; j<(int)NET.nLag; j++) {
900 
901  NET.getNetworkPixels(j,Eo);
902  pwc = NET.getwc(j);
903  if(NET.pattern!=0) {
904  NET.cluster(2,3);
905  wc.cpf(*(pwc),false);
906  wc.select(const_cast<char*>("subrho"),SELECT_SUBRHO);
907  wc.select(const_cast<char*>("subnet"),SELECT_SUBNET);
908  pwc->cpf(wc,false);
909  } else NET.cluster(1,1);
910  // store cluster into temporary job file
911  int cycle = cfg.simulation ? ifactor : Long_t(pwc->shift);
912  pwc->write(jfile,"coherence","clusters",0,cycle);
913  pwc->write(jfile,"coherence","clusters",-1,cycle,-(rateANA>>(cfg.l_high-i)));
914  cout<<cycle<<"|"<<pwc->csize()<<"|"<<pwc->size()<<" ";cout.flush();
915  csize_tot+=pwc->csize(); psize_tot+=pwc->size();
916 
917  // add core pixels to sparse table
918  for(int n=0; n<nIFO; n++) vSS[n].AddCore(n,pwc);
919 
920  pwc->clear();
921  }
922  // write clusters to the job file & free trees memory
923  jfile->Write();
924  TList* fList = gDirectory->GetList();
925  TObject *obj;
926  TIter nextobj(fList);
927  while ((obj = (TObject *) nextobj())) {
928  if(TString(obj->GetName()).Contains("clusters")) delete obj;
929  }
930  // update sparse table (halo added to core pixels)
931  for(int n=0; n<nIFO; n++) {vSS[n].UpdateSparseTable();}
932  // write sparse table to job file
933  sdir->cd();
934  for(int n=0; n<nIFO; n++) {
935  vSS[n].Clean();
936  char wsname[32];
937  if(cfg.simulation) sprintf(wsname,"%s-level:%d:%d",ifo[n],ifactor,cfg.l_high-i);
938  else sprintf(wsname,"%s-level:%d",ifo[n],cfg.l_high-i);
939  vSS[n].Write(wsname,TObject::kWriteDelete);
940  }
941  // print benchmark infos
942  watchCoh.Stop(); cout<<endl;
943  PrintElapsedTime(watchCoh.RealTime(),"Coherence Elapsed Time for this level : ");
944  cout<<endl; watchCoh.Reset();
945  // add infos to history
946  sprintf(info,"-RES:%d-CSIZE:%d-PSIZE:%d",i,(int)csize_tot/NET.nLag,(int)psize_tot/NET.nLag);
947  PrintAnalysisInfo(CWB_STAGE_COHERENCE,"cwb2G::Coherence",info,false);
948  }
949  }
950 
951  if(bplugin) {
952  char sfactor[8];sprintf(sfactor,"%d",ifactor);
953  CWB_Plugin(jfile,&cfg,&NET,NULL,sfactor,CWB_PLUGIN_OCOHERENCE);
954  }
955 
956  jfile->Close();
957 
958  return;
959 }
960 
961 void
963 //
964 // Multi resolution clustering & Rejection of the sub-threshold clusters
965 //
966 // Loop over time lags
967 // - Read clusters from job file (netcluster::read)
968 // - Multi resolution clustering (netcluster::supercluster)
969 // - Compute for each pixel the time delay amplitudes (netcluster::loadTDampSSE)
970 // - Rejection of the sub-threshold clusters (network::subNetCut)
971 // - Defragment clusters (netcluster::defragment)
972 // - Store superclusters to job file (netcluster::write)
973 // Build & Write to job file the sparse TF maps (WriteSparseTFmap)
974 //
975 
976  // if 2G analysis istage>=CWB_STAGE_SUPERCLUSTER then skip super cluster analysis
977  if(istage>=CWB_STAGE_SUPERCLUSTER) return;
978 
979  PrintStageInfo(CWB_STAGE_SUPERCLUSTER,"cwb2G::SuperCluster");
980 
981  char info[256];
982 
983  // open input job file
984  TFile* ifile = ((istage==CWB_STAGE_COHERENCE)&&(iname!="")) ? new TFile(iname) : NULL;
985  // open temporary job file
986  jfile = new TFile(jname,"UPDATE");
987  if(jfile==NULL||!jfile->IsOpen())
988  {cout << "cwb2G::SuperCluster - Error opening : " << jname << endl;EXIT(1);}
989 
990  // create prod/sim/lag directories in the output root file to store dignostic histograms
991  if(froot==NULL) {cout << "cwb2G::SuperCluster - Error opening : " << froot->GetPath() << endl;EXIT(1);}
992  TDirectory* wdir = (TDirectory*)froot->Get("histogram");
993  if(wdir==NULL) wdir = froot->mkdir("histogram");
994 
995  // decrease skymap resolution to improve subNetCut performances
996  double skyres=0;
998  else skyres = cfg.angle<MIN_SKYRES_ANGLE ? MIN_SKYRES_ANGLE : 0;
999  if(skyres) {
1000  if(cfg.healpix) NET.setSkyMaps(int(skyres));
1001  else NET.setSkyMaps(skyres,cfg.Theta1,cfg.Theta2,cfg.Phi1,cfg.Phi2);
1002  NET.setAntenna();
1004  // the down resampling of the skymask works only for the built-in skymask
1005  if(strlen(cfg.skyMaskFile)>0) SetSkyMask(&NET,&cfg,cfg.skyMaskFile,'e',skyres);
1006  if(strlen(cfg.skyMaskCCFile)>0) SetSkyMask(&NET,&cfg,cfg.skyMaskCCFile,'c',skyres);
1007  }
1008 
1009  for(int i=0; i<nIFO; i++) pTF[i] = pD[i]->getTFmap();
1010  // set low-rate TD filters
1011  for(int k=0;k<nRES;k++) pwdm[k]->setTDFilter(cfg.TDSize, 1);
1012  // read sparse map from job file
1013  cout << "Loading sparse TF map ... " << endl;
1014  for(int n=0; n<nIFO; n++) {
1015  pD[n]->sclear(); // clear vector with sparse maps
1016  for(int i=0; i<nRES; i++) {
1017  char swname[32];
1018  if(cfg.simulation) sprintf(swname,"csparse/%s-level:%d:%d",ifo[n],ifactor,i+cfg.l_low);
1019  else sprintf(swname,"csparse/%s-level:%d",ifo[n],i+cfg.l_low);
1020  SSeries<double>* psw;
1021  if(ifile!=NULL) psw = (SSeries<double>*)ifile->Get(swname);
1022  else psw = (SSeries<double>*)jfile->Get(swname);
1023  if(psw==NULL) {
1024  cout << "cwb2G::SuperCluster : sparse map " << swname
1025  << " not exist in job file" << endl;EXIT(1);
1026  }
1027  SSeries<double> SS = *psw;
1028  pD[n]->vSS.push_back(SS);
1029  delete psw;
1030  }
1031  cout<<endl;
1032  }
1033 
1034  // export to CINT ifile (used by plugins)
1035  char cmd[128];
1036  if(ifile!=NULL) sprintf(cmd,"gIFILE = (void*)%p;",ifile);
1037  else sprintf(cmd,"gIFILE = NULL;");
1038  EXPORT(void*,gIFILE,cmd);
1039  // in supercluster plugin
1041 
1042  if(!cfg.scPlugin) { // built in supercluster function
1043 
1044  int nevt = 0;
1045  int nnn = 0;
1046  int mmm = 0;
1047  size_t count = 0;
1048  netcluster wc;
1049  netcluster* pwc;
1050 
1051  for(int j=0; j<(int)lags; j++) {
1052 
1053  int cycle = cfg.simulation ? ifactor : Long_t(NET.wc_List[j].shift);
1054 
1055  // read cluster metadata
1056  if(ifile!=NULL) wc.read(ifile,"coherence","clusters",0,cycle);
1057  else wc.read(jfile,"coherence","clusters",0,cycle);
1058  // read clusters from temporary job file, loop over TF resolutions
1059  if(ifile!=NULL) {
1060  for(int i=nRES-1; i>=0; i--) // reverse loop is faster loading cluster (?)
1061  wc.read(ifile,"coherence","clusters",-2,cycle,-(rateANA>>(i+cfg.l_low)));
1062  } else {
1063  for(int i=nRES-1; i>=0; i--) // reverse loop is faster loading cluster (?)
1064  wc.read(jfile,"coherence","clusters",-2,cycle,-(rateANA>>(i+cfg.l_low)));
1065  }
1066  cout<<"-----------------------------------------------------"<<endl;
1067  char cycle_name[32];
1068  if(cfg.simulation) sprintf(cycle_name," factor[%d]=%g",ifactor,cfg.factors[ifactor]);
1069  else sprintf(cycle_name," lag=%d",cycle);
1070  cout<<"-> Processing " <<cycle_name<<" ..."<<endl;
1071  cout<<" --------------------------------------------------"<<endl;
1072  cout<<" coher clusters|pixels : "
1073  <<setfill(' ')<<setw(6)<<wc.csize()<<"|"<<wc.size()<<endl;
1074 
1075  // supercluster analysis
1076  if(cfg.l_high==cfg.l_low) wc.pair=false; // if only one resolution is used pair is false
1077  if(NET.pattern!=0) wc.pair=false; // if other than pattern=0 - allow one resolution cluster
1078  wc.supercluster('L',NET.e2or,cfg.TFgap,false); // likehood2G
1079  cout<<" super clusters|pixels : "
1080  <<setfill(' ')<<setw(6)<<wc.esize(0)<<"|"<<wc.psize(0)<<endl;
1081 
1082  // defragmentation for pattern != 0
1083  if(NET.pattern!=0) {
1084  wc.defragment(cfg.Tgap,cfg.Fgap);
1085  cout<<" defrag clusters|pixels : "
1086  <<setfill(' ')<<setw(6)<<wc.esize(0)<<"|"<<wc.psize(0)<<"\n";
1087  }
1088 
1089  // copy selected clusters to network
1090  pwc = NET.getwc(j);
1091  pwc->cpf(wc, false);
1092 
1093  // apply subNetCut() only for pattern=0 || cfg.subnet>0 || cfg.subcut>0
1094  if(NET.pattern==0 || cfg.subnet>0 || cfg.subcut>0) {
1095  NET.setDelayIndex(hot[0]->rate());
1096  pwc->setcore(false);
1097  int psel = 0;
1098  while(1) {
1099  count = pwc->loadTDampSSE(NET, 'a', cfg.BATCH, cfg.LOUD);
1100  psel += NET.subNetCut((int)j,cfg.subnet,cfg.subcut,NULL);
1101  int ptot = pwc->psize(1)+pwc->psize(-1);
1102  double pfrac = ptot>0 ? double(psel)/double(ptot) : 0.;
1103  //cout<<"selected pixels: "<<psel<<", fraction: "<<pfrac<< endl;
1104  if(count<10000) break;
1105  }
1106  cout<<" subnet clusters|pixels : "
1107  <<setfill(' ')<<setw(6)<<NET.events()<<"|"<<pwc->psize(-1)<<"\n";
1108  }
1109  if(NET.pattern==0) {
1110  // defragmentation
1111  pwc->defragment(cfg.Tgap,cfg.Fgap);
1112  cout<<" defrag clusters|pixels : "
1113  <<setfill(' ')<<setw(6)<<NET.events()<<"|"<<pwc->psize(-1)<<"\n";
1114  }
1115 
1116  nevt += NET.events();
1117  nnn += pwc->psize(-1);
1118  mmm += pwc->psize(1)+pwc->psize(-1);
1119 
1120  // store cluster into temporary job file [NEWSS]
1121  pwc->write(jfile,"supercluster","clusters",0,cycle);
1122  pwc->write(jfile,"supercluster","clusters",-1,cycle);
1123  //cout<<cycle<<"|"<<pwc->csize()<<"|"<<pwc->size()<<" ";cout.flush();
1124 
1125  pwc->clear();
1126  cout<<endl;cout.flush();
1127  }
1128 
1129  // print final statistic
1130  cout<<endl<<"Supercluster done"<<endl;
1131  if(mmm) cout<<"total clusters|pixels|frac : "
1132  <<setfill(' ')<<setw(6)<<nevt<<"|"<<nnn<<"|"<<nnn/double(mmm)<<"\n";
1133  else cout<<"total clusters : "<<nevt<<"\n";
1134  cout<<endl;cout.flush();
1135 
1136  // add infos to history
1137  if(mmm) sprintf(info,"-NEVT:%d-PSIZE:%d-FRAC:%g",nevt,nnn,nnn/double(mmm));
1138  else sprintf(info,"-NEVT:%d",nevt);
1139  PrintAnalysisInfo(CWB_STAGE_SUPERCLUSTER,"cwb2G::SuperCluster",info,false);
1140  }
1141 
1142  // out supercluster plugin
1144 
1145  // restore skymap resolution
1146  if(skyres) {
1147  if(cfg.healpix) NET.setSkyMaps(int(cfg.healpix));
1149  NET.setAntenna();
1151  if(strlen(cfg.skyMaskFile)>0) SetSkyMask(&NET,&cfg,cfg.skyMaskFile,'e');
1152  if(strlen(cfg.skyMaskCCFile)>0) SetSkyMask(&NET,&cfg,cfg.skyMaskCCFile,'c');
1153  }
1154 
1155  if(ifile!=NULL) {
1156  for(int i=0; i<nIFO; i++) {
1158  LoadWaveforms(ifile, pD[i], ifactor); // read inj waveforms from input job file
1159  SaveWaveforms(jfile, pD[i], ifactor); // write inj waveforms to output job file
1160  }
1161  }
1162  }
1163 
1164  jfile->Write();
1165  jfile->Close();
1166 
1167  // clear vector with sparse maps
1168  for(int n=0; n<nIFO; n++) pD[n]->sclear();
1169  // write sparse table to job file
1170  jfile = new TFile(jname,"UPDATE","",9); // compression level 9
1171  WriteSparseTFmap(jfile, ifactor, "sparse", "supercluster");
1172  jfile->Close();
1173 
1174  // job root file garbage collector
1175  vector<TString> delObjList;
1176  // coherence clusters is removed if not set in the jobfOptions
1178  delObjList.push_back("coherence");
1179  }
1180  // cstrains and rms data are removed if not set in the jobfOptions
1182  char cdrms_name[32];sprintf(cdrms_name,"rms/rms-f%d",ifactor);
1183  char cdcstrain_name[32];sprintf(cdcstrain_name,"cstrain/cstrain-f%d",ifactor);
1184  delObjList.push_back(cdrms_name);
1185  delObjList.push_back(cdcstrain_name);
1186  }
1187  // sparse maps are removed if not set in the jobfOptions (only for the last factor)
1188  if(ifactor==cfg.nfactor-1) { // the last factor
1189  if(!(jobfOptions&CWB_JOBF_SAVE_CSPARSE)) delObjList.push_back("csparse");
1190  }
1191  FileGarbageCollector(jname,"",delObjList);
1192 
1193  return;
1194 }
1195 
1196 bool
1198 //
1199 // Event reconstruction & parameters estimation
1200 //
1201 // Read sparse map from job file
1202 // Loop over time lags
1203 // - Read cluster list from job file (netcluster::read)
1204 // - Loop over cluster list
1205 // - Read pixels (netcluster::read)
1206 // - Compute for each pixel the time delay amplitudes (netcluster::loadTDampSSE)
1207 // - Event reconstruction+parameter estimation (network::likelihood2G)
1208 // - Store event parameters to job file (netevent::output2G)
1209 // - If(config::cedDump>0) Generate Coherent Event Display (CWB::ced)
1210 //
1211 
1212  // if 2G analysis istage>=CWB_STAGE_LIKELIHOOD then skip likelihood analysis
1213  if(istage>=CWB_STAGE_LIKELIHOOD) return false;
1214 
1215  PrintStageInfo(CWB_STAGE_LIKELIHOOD,"cwb2G::Likelihood");
1216 
1217  char info[256];
1218  netcluster* pwc;
1219  int ceddir = 0; // flag if ced directory exists
1220  double factor = cfg.factors[ifactor];
1221 
1223  if(singleDetector) NET.setIndexMode(cfg.mode); // when nIFO=1 exclude duplicate delay configurations
1224 
1225  if(bplugin) {
1226  // export to CINT variables
1227  EXPORT(float*,gSLAGSHIFT,TString::Format("gSLAGSHIFT = (float*)%p;",slagShift).Data());
1228  }
1229 
1230  // open temporary job file
1231  TString xname = ((istage==CWB_STAGE_SUPERCLUSTER)&&(iname!="")) ? iname : jname;
1232  jfile = (xname==jname) ? new TFile(xname,"UPDATE") : new TFile(xname);
1233  if(jfile==NULL||!jfile->IsOpen())
1234  {cout << "cwb2G::Likelihood - Error : file " << xname.Data() << " not found" << endl;EXIT(1);}
1235 
1237  // restore in detector object the inj waveforms from input job file
1238  for(int i=0; i<nIFO; i++) LoadWaveforms(jfile, pD[i], ifactor);
1239  }
1241  // save to job file whitened inj/rec waveforms
1242  TFile* ifile = jfile;
1243  if(xname!=jname) {
1244  ifile = new TFile(jname,"UPDATE");
1245  if(ifile==NULL||!ifile->IsOpen()) {
1246  cout << "cwb2G::Likelihood - Error : file " << jname << " not found" << endl; EXIT(1); }
1247  }
1248  int wfSAVE = 0;
1249  if(jobfOptions&CWB_JOBF_SAVE_WFINJ) wfSAVE+=1;
1250  if(jobfOptions&CWB_JOBF_SAVE_WFREC) wfSAVE+=2;
1251  for(int i=0; i<nIFO; i++) SaveWaveforms(ifile, pD[i], ifactor, wfSAVE);
1252  ifile->Write();
1253  if(xname!=jname) ifile->Close();
1254  }
1255 
1256  // set low-rate TD filters
1257  for(int k=0; k<nRES; k++) pwdm[k]->setTDFilter(cfg.TDSize, cfg.upTDF);
1259 
1260  // read sparse map from job file
1261  cout << "Loading sparse TF map ... " << endl;
1262  for(int n=0; n<nIFO; n++) {
1263  pD[n]->sclear(); // clear vector with sparse maps
1264  for(int i=0; i<nRES; i++) {
1265  char swname[32];
1266  if(cfg.simulation) sprintf(swname,"sparse/%s-level:%d:%d",ifo[n],ifactor,i+cfg.l_low);
1267  else sprintf(swname,"sparse/%s-level:%d",ifo[n],i+cfg.l_low);
1268  SSeries<double>* psw = (SSeries<double>*)jfile->Get(swname);
1269  if(psw==NULL) {
1270  cout << "cwb2G::Likelihood : sparse map " << swname
1271  << " not exist in job file" << endl;EXIT(1);
1272  }
1273  SSeries<double> SS = *psw;
1274  pD[n]->vSS.push_back(SS);
1275  delete psw;
1276  }
1277  cout<<endl;
1278  }
1279 
1280  if(cfg.dump && (!cfg.outPlugin)) // write header to ascii dump file
1281  {netburst->dopen(outDump,const_cast<char*>("a"),true); netburst->dclose();}
1282 
1283  int LRETRY = 0; // number of likelihood2G retry for each event
1284  int NEVENTS = 0; // total recontructed events for all lags
1285  for(int j=0; j<(int)lags; j++) {
1286 
1287  // init likelihood plugin
1288  EXPORT(int,gLRETRY,"gLRETRY = 0;") // needed to avoid error msg when gLRETRY is not defined in the plugin
1289  if(bplugin) {
1290  EXPORT(int,gRET,"gRET = 0;") // needed to avoid error msg when gRET is not defined in the plugin
1291  char strcycle[8];sprintf(strcycle,"%d",cfg.simulation==0 ? j : ifactor);
1292  CWB_Plugin(jfile,&cfg,&NET,NULL,strcycle,CWB_PLUGIN_ILIKELIHOOD);
1293  // get gRET from plugin
1294  int gRET=0; IMPORT(int,gRET) if(gRET==-1) continue;
1295  }
1296  int gLRETRY=0; IMPORT(int,gLRETRY) LRETRY=gLRETRY;
1297 
1298  int cycle = cfg.simulation ? ifactor : Long_t(NET.wc_List[j].shift);
1299 
1300  pwc = NET.getwc(j);
1301  pwc->cData.clear();
1302 
1303  // read cluster list & metadata netcluster object
1304  vector<int> clist = pwc->read(jfile,"supercluster","clusters",0,cycle);
1305  //pwc->print();
1306 
1307  // print header
1308  char cycle_name[32];
1309  if(cfg.simulation) sprintf(cycle_name," factor[%d]=%g",ifactor,cfg.factors[ifactor]);
1310  else sprintf(cycle_name," lag=%d",cycle);
1311  cout<<"-------------------------------------------------------"<<endl;
1312  cout<<"-> Processing " << clist.size() << " clusters in" << cycle_name<<endl;
1313  cout<<" ----------------------------------------------------"<<endl;cout.flush();
1314 
1315  //int nmax = 1000; // load tdAmp associated to the first nmax loudest pixels
1316  int nmax = -1; // load all tdAmp
1317  int npixels = 0; // total loaded pixels per lag
1318  int nevents = 0; // total recontructed events per lag
1319  int nselected_core_pixels = 0;
1320  int nrejected_weak_pixels = 0; // remove weak glitches
1321  int nrejected_loud_pixels = 0; // remove loud glitches
1322  for(int k=0;k<(int)clist.size();k++) { // loop over the cluster list
1323 
1324  // read pixels & tdAmp into netcluster pwc
1325  pwc->read(jfile,"supercluster","clusters",nmax,cycle,0,clist[k]);
1326 
1327  // likelihood plugin
1329 
1330  wavearray<double> cid = pwc->get((char*)"ID", 0,'S',0); // get cluster ID
1331  int id = size_t(cid.data[cid.size()-1]+0.1);
1332  pwc->setcore(false,id);
1333  pwc->loadTDampSSE(NET, 'a', cfg.BATCH, cfg.BATCH); // attach TD amp to pixels
1334 
1335  int lag = j;
1336 
1337  NET.MRA = true;
1338  int ID = cfg.cedDump ? -id : 0;
1339  int selected_core_pixels = 0;
1340  if(NET.pattern>0) {
1341  selected_core_pixels = NET.likelihoodWP(cfg.search, lag, ID, NULL);
1342  } else {
1343  selected_core_pixels = NET.likelihood2G(cfg.search, lag, ID, NULL);
1344  }
1345  if(!cfg.outPlugin) { // if true then output to root file is provided by the user plugin
1346  double ofactor=0;
1347  if(cfg.simulation==4) ofactor=-factor;
1348  else if(cfg.simulation==3) ofactor=-ifactor;
1349  else ofactor=factor;
1350  if(cfg.dump) netburst->dopen(outDump,const_cast<char*>("a"),false);
1351  netburst->output2G(net_tree,&NET,id,lag,ofactor);
1352  if(cfg.dump) netburst->dclose();
1353  }
1354  int rejected_weak_pixels = 0;
1355  int rejected_loud_pixels = 0;
1356 
1357  bool detected = (bool)(NET.getwc(j)->sCuts[k] == -1);
1358 
1359  // print reconstructed event
1360  cout<<" cluster-id|pixels: "<<setfill(' ')<<setw(5)<<clist[k]<<"|"<<pwc->size()-npixels;
1361  if(detected) cout << "\t -> SELECTED !!!" << endl;
1362  else cout << "\t <- rejected " << endl;
1363  cout.flush();
1364 
1365  // likelihood clusters are stored into temporary job file if set in the jobfOptions
1366  if(((k==0)||detected)&&(jobfOptions&CWB_JOBF_SAVE_LIKELIHOOD)) {
1367  TFile* ifile = jfile;
1368  if(xname!=jname) {
1369  ifile = new TFile(jname,"UPDATE");
1370  if(ifile==NULL||!ifile->IsOpen()) {
1371  cout << "cwb2G::Likelihood - Error : file " << jname << " not found" << endl; EXIT(1); }
1372  }
1373  pwc->write(ifile,"likelihood","clusters",0,cycle);
1374  pwc->write(ifile,"likelihood","clusters",-1,cycle,0,k+1);
1375  if(detected) cout<<"saved"<<endl;cout.flush();
1376  ifile->Write();
1377  if(xname!=jname) ifile->Close();
1378  }
1379 
1380  if(detected) nevents++;
1381  if(gLRETRY==0) npixels=pwc->size();
1382  if(!detected && gLRETRY==LRETRY) {
1383  gLRETRY=0;
1384  EXPORT(int,gLRETRY,TString::Format("gLRETRY = %d;",gLRETRY).Data())
1385  }
1386 
1387  // likelihood plugin
1389 
1390  // ced dump
1391  if(detected&&(cfg.cedDump)) {
1392  TFile* ifile = jfile;
1393  if(xname!=jname) {
1394  ifile = new TFile(jname,"UPDATE");
1395  if(ifile==NULL||!ifile->IsOpen()) {
1396  cout << "cwb2G::Likelihood - Error : file " << jname << " not found" << endl; EXIT(1); }
1397  }
1398  ced = NULL;
1400  // save ced to temporary job file
1401  cout<<"dump ced to job file ..." <<endl;
1402  TDirectory* cdced = NULL;
1403  cdced = (TDirectory*)ifile->Get("ced");
1404  if(cdced == NULL) cdced = ifile->mkdir("ced");
1405  ced = new CWB::ced(&NET,netburst,cdced);
1406  } else {
1407  cout<<"dump ced to disk ..." <<endl;
1408  ced = new CWB::ced(&NET,netburst,ced_dir);
1409  }
1410  switch(istage) {
1411  case CWB_STAGE_FULL:
1412  case CWB_STAGE_INIT:
1413  // use TF map & inj signals
1415  for(int n=0; n<nIFO; n++) {pTF[n]->setlow(cfg.fLow); pTF[n]->sethigh(cfg.fHigh);}
1416  break;
1417  default:
1418  // use sparse map & inj signals
1420  }
1422  bool fullCED = singleDetector ? false : true;
1423  double ofactor=0;
1424  // CED plugin
1425  if(bplugin) CWB_Plugin(jfile,&cfg,&NET,NULL,"",CWB_PLUGIN_CED);
1426  if(cfg.simulation==4) ofactor=-factor;
1427  else if(cfg.simulation==3) ofactor=-ifactor;
1428  else ofactor=factor;
1429  if(ced->Write(ofactor,fullCED)) ceddir = 1;
1430  delete ced;
1431  if(xname!=jname) ifile->Close();
1432  }
1433 
1434  //pwc->print();
1435  pwc->clean(k+1); // clean time delay amplitudes
1436 
1437  IMPORT(int,gLRETRY)
1438  if(gLRETRY-- > 0) { // clear last cluster in pwc
1439  clist = pwc->read(jfile,"supercluster","clusters",0,cycle);
1440  if(pwc->cList.size()) {
1441  int npix = pwc->cList[id-1].size();
1442  for(int p=0;p<npix;p++) pwc->pList.pop_back();
1443  pwc->cList.pop_back(); pwc->cData.pop_back(); pwc->sCuts.pop_back();
1444  pwc->cRate.pop_back(); pwc->cTime.pop_back(); pwc->cFreq.pop_back();
1445  pwc->sArea.pop_back(); pwc->p_Map.pop_back(); pwc->p_Ind.pop_back();
1446  }
1447  k--;
1448  } else gLRETRY=LRETRY; // every new event gLRETRY is reinitialize to LRETRY
1449  EXPORT(int,gLRETRY,TString::Format("gLRETRY = %d;",gLRETRY).Data())
1450  } // end loop over the cluster list
1451  /*
1452  // print final lag statistic
1453  if(nevents==0)
1454  cout << "--------------------------------------------------------------" << endl;
1455  cout<<endl;
1456  cout<<" total pixels : "<<pwc->size()<<endl;
1457  cout<<" selected_core_pixels : "<<nselected_core_pixels<<endl;
1458  cout<<" rejected_weak_pixels : "<<nrejected_weak_pixels<<endl;
1459  cout<<" rejected_loud_pixels : "<<nrejected_loud_pixels <<endl;
1460  cout<<"-> reconstruct events : "<<nevents<<endl;
1461  cout<<endl;
1462  */
1463  cout<<endl<<endl;cout.flush();
1464  NEVENTS+=nevents;
1465  } // end loop over lags
1466 
1467  // add infos to history
1468  sprintf(info,"-NEVT:%d",NEVENTS);
1469  PrintAnalysisInfo(CWB_STAGE_LIKELIHOOD,"cwb2G::Likelihood",info,false);
1470 
1471  jfile->Close();
1472 
1473  // sparse maps & supercluster clusters are removed if not set in the jobfOptions (only for the last factor)
1474  if(ifactor==cfg.nfactor-1) { // the last factor
1475  vector<TString> delObjList;
1476  // supercluster clusters are removed if not set in the jobfOptions
1477  if(!(jobfOptions&CWB_JOBF_SAVE_SUPERCLUSTER)) delObjList.push_back("supercluster");
1478  if(!(jobfOptions&CWB_JOBF_SAVE_SPARSE)) delObjList.push_back("sparse");
1479  // If CED is saved to the job file than FileGarbageCollector is not applied
1480  // FileGarbageCollector do not preserve the style of the plots
1482  }
1483 
1484  return ceddir;
1485 }
1486 
1487 void
1489 //
1490 // Fill sparse map with cluster pixels and write to job file
1491 //
1492 // jfile : pointer to job file
1493 // ifactor : factor ID
1494 // tdir : directory name which contains the sparse map
1495 // tname : tree name : Ex "coherence" or "supercluster"
1496 //
1497 
1498  // create sparse directories in the job root file to store sparse TF maps
1499  TDirectory* sdir = (TDirectory*)jfile->Get(tdir);
1500  if(sdir==NULL) sdir = jfile->mkdir(tdir);
1501 
1502  int ncore_tot=0;
1503  int ncluster_tot=0;
1504  int ccluster_tot=0;
1505  netcluster wc;
1507  std::vector<SSeries<double> > vSS;
1508  for(int n=0; n<nIFO; n++) vSS.push_back(ss);
1509 
1510  for(int i=0; i<nRES; i++) {
1511  // init sparse table
1512  pwdm[i]->setTDFilter(cfg.TDSize, 1);
1513  for(int n=0; n<nIFO; n++) {
1514  pTF[n] = pD[n]->getTFmap();
1515  pTF[n]->Forward(*hot[n],*pwdm[i]);
1516  vSS[n].SetMap(pTF[n]);
1517  vSS[n].SetHalo(mTau);
1518  if(singleDetector) {
1519  pTF[1] = pD[1]->getTFmap();
1520  *pTF[1] = *pTF[0];
1521  vSS[1].SetMap(pTF[1]);
1522  vSS[1].SetHalo(mTau);
1523  break;
1524  }
1525  }
1526 
1527  for(int j=0; j<(int)lags; j++) {
1528 
1529  int cycle = cfg.simulation ? ifactor : Long_t(NET.wc_List[j].shift);
1530 
1531  // read cluster metadata
1532  wc.read(jfile,tname,"clusters",0,cycle);
1533 
1534  // read pixels & tdAmp into netcluster pwc
1535  if(tname=="coherence") {
1536  wc.read(jfile,tname,"clusters",-2,cycle,-vSS[0].wrate());
1537  } else {
1538  wc.read(jfile,tname,"clusters",-1,cycle,vSS[0].wrate());
1539  }
1540  //cout<<"loaded clusters|pixels: "<<wc.csize()<<"|"<<wc.size()<<endl;
1541 
1542  // add core pixels to sparse table
1543  for(int n=0; n<nIFO; n++) vSS[n].AddCore(n,&wc);
1544  wc.clear();
1545  }
1546 
1547  // update sparse table (halo added to core pixels)
1548  for(int n=0; n<nIFO; n++) {vSS[n].UpdateSparseTable();}
1549 
1550  // compute sparse statistic
1551  int ncore=0; // core pixels
1552  int ncluster=0; // core+halo pixels
1553  int ccluster=0; // core+halo pixels associated to each core pixels
1554  for(int n=0; n<nIFO; n++) {ncore+=vSS[n].GetSparseSize();ncluster+=vSS[n].GetSparseSize(0);}
1555  ccluster = 2*(vSS[0].GetHaloSlice()+vSS[0].GetHaloSlice(true))+1;
1556  ccluster*= 2*vSS[0].GetHaloLayer()+1;
1557  ncore_tot+=ncore; ncluster_tot+=ncluster; ccluster_tot+=2*ccluster*ncore;
1558  int rate = rateANA/(rateANA>>(cfg.l_high-i));
1559  /*
1560  // print sparse statistic (per resolution)
1561  cout << "----------------------- SPARSE STATISTIC / DETECTOR / RES -----------------------" << endl;
1562  cout << "rate|sSlice|sLayer|eSlice|ccluster : " << rate << "|" << vSS[0].GetHaloSlice() << "|"
1563  << vSS[0].GetHaloLayer() << "|" <<vSS[0].GetHaloSlice(1) << "|" << ccluster << endl;
1564  cout << "npix_core|npix_cluster|3*npix_cluster|2*ccluster*npix_core|ratio|sfilesize(byte) : " << endl;
1565  cout << ncore << " | " << ncluster << " | " << 3*ncluster
1566  << " | " << 2*ccluster*ncore << " | " << ((ccluster*ncore)?double(3*ncluster)/(2*ccluster*ncore):0)
1567  << " | " << 3*ncluster*4 << endl;
1568  cout<<endl;
1569  */
1570  // write sparse table to job file
1571  sdir->cd();
1572  for(int n=0; n<nIFO; n++) {
1573  if(!(cfg.cedDump && jstage==CWB_STAGE_FULL)) vSS[n].Clean(); // TF map not cleaned if used by ced
1574  char wsname[32];
1575  if(cfg.simulation) sprintf(wsname,"%s-level:%d:%d",ifo[n],ifactor,cfg.l_high-i);
1576  else sprintf(wsname,"%s-level:%d",ifo[n],cfg.l_high-i);
1577  vSS[n].Write(wsname,TObject::kWriteDelete);
1578  }
1579  }
1580 
1581  // print final sparse statistic
1582  cout << "----------------- FINAL SPARSE STATISTIC / DETECTOR -------------------" << endl;
1583  cout << "npix_core_tot|npix_cluster_tot|3*npix_cluster_tot|ccluster_tot|ratio : " << endl;
1584  cout << ncore_tot << " | " << ncluster_tot << " | " << 3*ncluster_tot
1585  << " | " << ccluster_tot << " | " << (ccluster_tot?double(3*ncluster_tot)/(ccluster_tot):0) << endl;
1586 
1587  return;
1588 }
1589 
1590 void
1592 //
1593 // Fill sparse map with cluster pixels
1594 //
1595 // jfile : pointer to job file
1596 // ifactor : factor ID
1597 // tname : tree name : Ex "coherence" or "supercluster"
1598 //
1599 
1600  netcluster wc;
1601 
1602  SSeries<double> SS;
1603  for(int n=0; n<nIFO; n++) pD[n]->sclear(); // clear vector with sparse maps
1604 
1605  cout << "cwb2G::FillSparseTFmap ..." << endl;
1606  for(int i=0; i<nRES; i++) {
1607 
1608  // init sparse table
1609  for(int n=0; n<nIFO; n++) {
1610  pTF[n] = pD[n]->getTFmap();
1611  pTF[n]->Forward(*hot[n],*pwdm[i]);
1612  pD[n]->vSS.push_back(SS);
1613  pD[n]->vSS[i].SetMap(pTF[n]);
1614  pD[n]->vSS[i].SetHalo(mTau);
1615  if(singleDetector) {
1616  pTF[1] = pD[1]->getTFmap();
1617  *pTF[1] = *pTF[0];
1618  pD[1]->vSS.push_back(SS);
1619  pD[1]->vSS[i].SetMap(pTF[1]);
1620  pD[1]->vSS[i].SetHalo(mTau);
1621  break;
1622  }
1623  }
1624 
1625  for(int j=0; j<(int)lags; j++) {
1626 
1627  int cycle = cfg.simulation ? ifactor : Long_t(NET.wc_List[j].shift);
1628 
1629  // read cluster metadata
1630  wc.read(jfile,"coherence","clusters",0,cycle);
1631 
1632  // read pixels & tdAmp into netcluster pwc
1633  if(tname=="coherence") {
1634  wc.read(jfile,tname,"clusters",-2,cycle,-pD[0]->vSS[i].wrate());
1635  } else {
1636  wc.read(jfile,tname,"clusters",-1,cycle,pD[0]->vSS[i].wrate());
1637  }
1638  //cout<<"loaded clusters|pixels: "<<wc.csize()<<"|"<<wc.size()<<endl;
1639 
1640  // add core pixels to sparse table
1641  for(int n=0; n<nIFO; n++) pD[n]->vSS[i].AddCore(n,&wc);
1642 
1643  wc.clear();
1644  }
1645  // update sparse table (halo added to core pixels)
1646  for(int n=0; n<nIFO; n++) {pD[n]->vSS[i].UpdateSparseTable();}
1647  }
1648  return;
1649 }
1650 
1651 void
1652 cwb2G::SaveWaveforms(TFile* jfile, detector* pD, int ifactor, int wfSAVE) {
1653 //
1654 // write injected waveforms to job file (only for simulation>0)
1655 //
1656 // jfile : pointer to job file
1657 // pD : pointer to detector object
1658 // ifactor : factor ID
1659 // wfSAVE : 1/2/3
1660 // 1 : save injected waveforms
1661 // 2 : save reconstructed waveforms
1662 // 3 : save injected & reconstructed waveforms
1663 //
1664 
1665  if(cfg.simulation==0) return;
1666 
1667  if(jfile==NULL)
1668  {cout << "cwb2G::SaveWaveforms - NULL input root file : " << jname << endl;EXIT(1);}
1669  if(wfSAVE==0)
1670  {cout << "cwb2G::SaveWaveforms - save mode must be 1/2/3 " << endl;EXIT(1);}
1671 
1672  // create wf directory
1673  TDirectory* wf = (TDirectory*)jfile->Get("waveform");
1674  if(wf==NULL) wf=jfile->mkdir("waveform");
1675  char cdwf_name[32];sprintf(cdwf_name,"waveform-f%d",ifactor);
1676  TDirectory* cdwf = (TDirectory*)wf->Get(cdwf_name);
1677  if(cdwf==NULL) cdwf=wf->mkdir(cdwf_name);
1678  cdwf->cd();
1679 
1680  pD->wfsave(wfSAVE); // save waveforms
1681  pD->Write(pD->Name,TObject::kSingleKey);
1682  pD->wfsave(0); // restore default value
1683 
1684  return;
1685 }
1686 
1687 void
1688 cwb2G::LoadWaveforms(TFile* ifile, detector* pD, int ifactor, int wfSAVE) {
1689 //
1690 // restore in detector object the wf waveforms from input job file (only for simulation>0)
1691 //
1692 // ifile : pointer to job file
1693 // pD : pointer to detector object
1694 // ifactor : factor ID
1695 // wfSAVE : 1/2/3
1696 // 1 : restore injected waveforms
1697 // 2 : restore reconstructed waveforms
1698 // 3 : restore injected & reconstructed waveforms
1699 //
1700 
1701  if(cfg.simulation==0) return;
1702 
1703  if(ifile==NULL)
1704  {cout << "cwb2G::LoadWaveforms - NULL input root file : " << iname << endl;EXIT(1);}
1705  if(wfSAVE==0)
1706  {cout << "cwb2G::LoadWaveforms - save mode must be 1/2/3 " << endl;EXIT(1);}
1707 
1708  // read wf waveforms
1709  ifile->cd();
1710  char cdwf_name[32];sprintf(cdwf_name,"waveform/waveform-f%d",ifactor);
1711  detector* pd = (detector*)ifile->Get(TString(cdwf_name)+TString("/")+pD->Name);
1712  if(pd==NULL) return;
1713 
1714  pd->wfsave(wfSAVE); // load waveforms
1715  *pd >> *pD;
1716 
1717  delete pd;
1718 }
1719 
char channelNamesMDC[NIFO_MAX][50]
Definition: config.hh:311
monster wdmMRA
list of pixel pointers for MRA
Definition: network.hh:651
void SetOptions(int simulation, double rho, double inRate, bool useSparse=false, char *gtype=const_cast< char *>("png"), int paletteId=0)
Definition: ced.hh:88
void sethigh(double f)
Definition: wseries.hh:132
detector * getifo(size_t n)
param: detector index
Definition: network.hh:436
size_t rateANA
Definition: cwb.hh:279
char analysis[8]
Definition: config.hh:117
virtual size_t defragment(double T, double F, TH2F *=NULL)
T - maximum time gap in seconds F - maximum frequency gap in Hz.
Definition: netcluster.cc:1294
TString outDump
double iwindow
Definition: config.hh:200
virtual void resize(unsigned int)
Definition: wseries.cc:901
CWB::frame fr[2 *NIFO_MAX]
Definition: cwb.hh:251
void PrintElapsedTime(int job_elapsed_time, TString info)
Definition: cwb.cc:2000
void LoadWaveforms(TFile *ifile, detector *pD, int ifactor, int wfSAVE=1)
Definition: cwb2G.cc:1688
size_t write(const char *file, int app=0)
Definition: netcluster.cc:3008
size_t nLag
Definition: network.hh:573
#define REGRESSION_FILTER_LENGTH
Definition: cwb2G.cc:32
void white(double dT=0, int wtype=1, double offset=0., double stride=0.)
what it does: see algorithm description in wseries.hh
Definition: detector.hh:258
static size_t GetProcInfo(bool mvirtual=true)
Definition: cwb.cc:1847
std::vector< vector_int > cRate
Definition: netcluster.hh:398
void WriteSparseTFmap(TFile *jfile, int ifactor, TString tdir, TString tname)
Definition: cwb2G.cc:1488
void setMatrix(double edge=0., double f=1.)
Definition: regression.cc:425
size_t TDSize
Definition: config.hh:217
void SaveWaveforms(TFile *jfile, detector *pD, int ifactor, int wfSAVE=1)
Definition: cwb2G.cc:1652
void FileGarbageCollector(TString ifName, TString ofName="", vector< TString > delObjList=vector< TString >())
Definition: cwb.cc:2338
std::vector< netcluster > wc_List
Definition: network.hh:610
int stop
Definition: Toolbox.hh:94
size_t add(detector *)
param: detector structure return number of detectors in the network
Definition: network.cc:2559
bool singleDetector
used for the stage stuff
Definition: cwb.hh:277
WSeries< double > * pTF[NIFO_MAX]
pointers to detectors
Definition: cwb.hh:264
void setAntenna(detector *)
param: detector (use theta, phi index array)
Definition: network.cc:2846
std::vector< double > * getmdcTime()
Definition: network.hh:422
bool mdcPlugin
Definition: config.hh:365
wavearray< double > HoT
Definition: detector.hh:350
size_t setsnr(wavearray< double > &, std::vector< double > *, std::vector< double > *, double=5., double=8.)
Definition: detector.cc:1670
WDM< double > * pwdm[NRES_MAX]
Definition: cwb2G.hh:80
virtual void rate(double r)
Definition: wavearray.hh:141
double cedRHO
Definition: config.hh:298
float factor
char skyMaskFile[1024]
Definition: config.hh:308
size_t upTDF
Definition: config.hh:216
#define REGRESSION_MATRIX_FRACTION
Definition: cwb2G.cc:33
bool dataPlugin
Definition: config.hh:364
bool Likelihood(int ifactor, char *ced_dir, netevent *output=NULL, TTree *net_tree=NULL, char *outDump=NULL)
Definition: cwb2G.cc:1197
long getNetworkPixels(int LAG, double Eo, double DD=1., TH1F *hist=NULL)
Definition: network.cc:78
#define B
int n
Definition: cwb_net.C:28
TString iname
stage benchmark
Definition: cwb.hh:241
void CWB_Plugin(TFile *jfile, CWB::config *, network *, WSeries< double > *, TString, int)
COHERENCE.
wavearray< double > z
Definition: Test10.C:32
double mTau
Definition: cwb.hh:285
detector * pD[NIFO_MAX]
noise variability
Definition: cwb.hh:263
TString("c")
bool bplugin
Definition: cwb.hh:297
size_t setIndexMode(size_t=0)
Definition: network.cc:8105
bool cohPlugin
Definition: config.hh:367
#define EXIT(ERR)
Definition: cwb2G.cc:47
void bandpass(wavearray< DataType_t > &ts, double flow, double fhigh, int n=-1)
Definition: wseries.cc:313
TFile * jfile
output root file
Definition: cwb.hh:259
int ID
Definition: TestMDC.C:70
int start
Definition: Toolbox.hh:93
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1) *itime+ifreq;double time=itime *dT;double freq=(ifreq >0) ? ifreq *dF :dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
wavearray< double > get(char *name, size_t index=0, char atype='R', int type=1, bool=true)
param: string with parameter name param: index in the amplitude array, which define detector param: c...
Definition: netcluster.cc:2207
std::vector< wavearray< double > * > RWFP
Definition: detector.hh:382
wavearray< double > HRSS
Definition: detector.hh:371
double whiteWindow
Definition: config.hh:189
double fLow
Definition: config.hh:140
void add(const wavearray< DataType_t > &, int=0, int=0, int=0)
Definition: wavearray.cc:746
TFile * froot
Definition: cwb.hh:256
std::vector< vector_float > sArea
Definition: netcluster.hh:401
netcluster * pwc
Definition: cwb_job_obj.C:38
bool dcPlugin
Definition: config.hh:366
CWB_STAGE istage
Definition: cwb.hh:247
CWB_STAGE jstage
Definition: cwb.hh:248
bool cedDump
Definition: config.hh:297
void setlow(double f)
Definition: wseries.hh:125
std::vector< SSeries< double > > vSS[NIFO_MAX]
int layers
char ifo[NIFO_MAX][8]
Definition: cwb.hh:245
virtual double ReadData(double mdcShift, int ifactor)
Definition: cwb.hh:201
void Coherence(int ifactor)
Definition: cwb2G.cc:773
ss AddCore(ifoID,&wc)
waveform wf
size_t setcore(bool core, int id=0)
Definition: netcluster.cc:243
int nevt
virtual size_t supercluster(char atype, double S, bool core)
param: statistic: E - excess power, L - likelihood param: selection threshold T for likelihood cluste...
Definition: netcluster.cc:808
std::vector< vector_int > cList
Definition: netcluster.hh:397
size_t add(WSeries< double > &target, char *name, double fL=0., double fH=0.)
Definition: regression.cc:91
virtual void start(double s)
Definition: wavearray.hh:137
int j
Definition: cwb_net.C:28
i drho i
size_t lags
Definition: cwb.hh:293
std::vector< double > mdcTime
Definition: network.hh:614
void setTDFilter(int nCoeffs, int L=1)
Definition: WDM.cc:639
long subNetCut(int lag, float subnet=0.6, float subcut=0.33, TH2F *hist=NULL)
Definition: network.cc:1014
network NET
pointers to WSeries
Definition: cwb.hh:266
size_t cluster(int kt, int kf)
param: time gap in pixels return: number of reconstructed clusters
Definition: network.hh:321
int ncore
static bool checkFile(TString fName, bool question=false, TString message="")
Definition: Toolbox.cc:4670
double Phi1
Definition: config.hh:279
long likelihood2G(char mode, int lag, int ID, TH2F *hist=NULL)
Definition: network.cc:1387
bool outPlugin
Definition: config.hh:369
std::vector< vector_float > p_Map
Definition: netcluster.hh:402
void dopen(const char *fname, char *mode, bool header=true)
Definition: netevent.hh:312
double Edge
Definition: network.hh:578
Definition: cwb2G.hh:33
cout<< "Selected Pixels : "<< nPix<< endl;wc.cluster(1, 1);SSeries< double > ss
char ced_dir[512]
Definition: test_config1.C:154
double mdcShift
Definition: cwb_net.C:196
CWB::mdc * gMDC
fprintf(stdout,"start=%f duration=%f rate=%f\, x.start(), x.size()/x.rate(), x.rate())
void clean(int cID=0)
Definition: netcluster.hh:451
wavearray< double > w
Definition: Test1.C:27
double segEdge
Definition: config.hh:164
size_t mode
Definition: config.hh:275
virtual size_t size() const
Definition: wavearray.hh:145
float slagShift[20]
Definition: cwb.hh:292
double dcCal[NIFO_MAX]
Definition: config.hh:196
int ccluster
double Theta2
Definition: config.hh:278
size_t esize(int k=2)
Definition: netcluster.hh:153
size_t fResample
Definition: config.hh:142
long likelihoodWP(char mode, int lag, int ID, TH2F *hist=NULL)
Definition: network.cc:293
void output2G(TTree *, network *, size_t, int, double)
Definition: netevent.cc:711
TTree * net_tree
double setVeto(double=5.)
param: time window around injections
Definition: network.cc:3487
char filter_dir[1024]
Definition: config.hh:305
static TString GetMDCLog(TString log, int pos)
Definition: Toolbox.cc:2278
double shift
Definition: netcluster.hh:382
size_t w_mode
Definition: wseries.hh:458
int simulation
Definition: config.hh:199
int l_high
Definition: config.hh:156
double Fgap
Definition: config.hh:136
void PrintAnalysisInfo(CWB_STAGE stage, TString comment, TString info, bool out=true, bool log=true)
Definition: cwb.cc:2129
#define IMPORT(TYPE, VAR)
Definition: cwb.hh:69
WSeries< double > nRMS
Definition: detector.hh:358
int segID[20]
Definition: cwb.hh:291
char search
Definition: config.hh:121
wavearray< double > * px
static void resampleToPowerOfTwo(wavearray< double > &w)
Definition: Toolbox.cc:6602
double maxEnergy(wavearray< DataType_t > &ts, Wavelet &w, double=0, int=1, int=0, TH1F *=NULL)
Definition: wseries.cc:504
wavearray< double > TIME
Definition: detector.hh:375
void SuperCluster(int ifactor)
Definition: cwb2G.cc:962
i() int(T_cor *100))
std::vector< netpixel > pList
Definition: netcluster.hh:390
double TDRate
Definition: cwb2G.hh:78
wavearray< double > * getHoT()
param: no parameters
Definition: detector.hh:175
gwavearray< double > * gx
void setDelayIndex(double rate)
param: MDC log file
Definition: network.cc:2896
int ifactor
std::vector< std::string > mdcList
Definition: network.hh:612
const int NIFO_MAX
Definition: wat.hh:22
void Resample(const wavearray< DataType_t > &, double, int=6)
Definition: wavearray.cc:640
#define MIN_SKYRES_HEALPIX
Definition: cwb2G.cc:28
int BATCH
Definition: config.hh:245
void apply(double threshold=0., char c='a')
Definition: regression.cc:709
void FillSparseTFmap(TFile *jfile, int ifactor, TString tname)
Definition: cwb2G.cc:1591
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
Definition: wseries.cc:193
printf("total live time: non-zero lags = %10.1f \, liveTot)
int l
bool detected
std::vector< vector_int > p_Ind
Definition: netcluster.hh:403
double dataShift[NIFO_MAX]
Definition: config.hh:205
std::vector< float > cFreq
Definition: netcluster.hh:400
std::vector< int > RWFID
Definition: detector.hh:381
int Write(double factor, bool fullCED=true)
Definition: ced.cc:607
char channelNamesRaw[NIFO_MAX][50]
Definition: config.hh:310
double precision
int gRATEANA
void SetChannelName(char *chName)
Definition: ced.hh:101
void PrintStageInfo(CWB_STAGE stage, TString comment, bool out=true, bool log=true, TString fname="")
Definition: cwb.cc:2020
int k
int pattern
Definition: network.hh:601
void solve(double th, int nE=0, char c='s')
Definition: regression.cc:610
char jname[1024]
job file object
Definition: cwb.hh:260
#define REGRESSION_SOLVE_EIGEN_THR
Definition: cwb2G.cc:34
int nfactor
Definition: config.hh:201
std::vector< wavearray< double > * > IWFP
Definition: detector.hh:379
size_t healpix
Definition: config.hh:282
size_t cpf(const netcluster &, bool=false, int=0)
Definition: netcluster.cc:117
char refIFO[4]
Definition: config.hh:125
double THRESHOLD(double bpp)
param: selected fraction of LTF pixels assuming Gaussian noise
Definition: network.cc:2615
size_t setFilter(size_t)
Definition: regression.cc:276
void setDelay(const char *="L1")
Definition: network.cc:2767
size_t loadTDampSSE(network &net, char c, size_t BATCH=10000, size_t LOUD=0)
Definition: netcluster.cc:3317
const int NRES_MAX
Definition: wat.hh:23
wavearray< double > TDUR
Definition: detector.hh:376
double e2or
Definition: network.hh:584
std::vector< int > IWFID
Definition: detector.hh:378
wavearray< double > getClean()
Definition: regression.hh:135
WSeries< double > wM
Definition: cwb_job_obj.C:41
TFile * ifile
frfile FRF[2 *NIFO_MAX]
Definition: cwb.hh:253
CWB::Toolbox TB
Definition: cwb.hh:243
size_t size()
Definition: netcluster.hh:147
int nevents
wavearray< double > BAND
Definition: detector.hh:374
int ncluster
int npix
int * xstart
double subnet
Definition: config.hh:247
size_t events()
Definition: network.hh:329
static int shiftBurstMDCLog(std::vector< std::string > &mdcList, vector< TString > ifos, double mdc_shift)
Definition: Toolbox.cc:2174
CWB::config cfg
Definition: cwb.hh:192
pWDM setTDFilter(12)
#define REGRESSION_APPLY_THR
Definition: cwb2G.cc:37
int l_low
Definition: config.hh:155
#define EXPORT(TYPE, VAR, CMD)
Definition: config.cc:92
#define SELECT_SUBRHO
Definition: cwb2G.cc:40
void sclear()
Definition: detector.hh:191
#define WDM_PRECISION
Definition: cwb2G.cc:45
double whiteStride
Definition: config.hh:190
std::vector< int > sCuts
Definition: netcluster.hh:392
void dclose()
Definition: netevent.hh:321
std::vector< float > cTime
Definition: netcluster.hh:399
double fHigh
Definition: config.hh:141
virtual double mean() const
Definition: wavearray.cc:1071
WSeries< double > * getTFmap()
param: no parameters
Definition: detector.hh:179
int LOUD
Definition: config.hh:246
#define SELECT_SUBNET
Definition: cwb2G.cc:41
void Init()
Definition: cwb2G.cc:61
Definition: Meyer.hh:36
char Name[16]
Definition: detector.hh:327
size_t csize()
Definition: netcluster.hh:151
unsigned int jobfOptions
history object
Definition: cwb.hh:276
int * layers
Definition: monster.hh:121
CWB::ced * ced
temporary time series
Definition: cwb2G.hh:83
static TString SetMDCLog(TString log, int pos, TString val)
Definition: Toolbox.cc:2238
int SetSkyMask(network *net, CWB::config *cfg, char *options, char skycoord, double skyres=-1)
Definition: cwb.cc:2484
double segLen
Definition: config.hh:161
void readFrames(char *filename, char *channel, wavearray< double > &w)
Definition: frame.cc:828
int tag
Definition: monster.hh:116
vector< TString > fList
Definition: LoopChirpMass.C:28
char cmd[1024]
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Definition: wseries.cc:246
netcluster * getwc(size_t n)
param: delay index
Definition: network.hh:439
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double subcut
Definition: config.hh:248
#define NEVENTS
int nRes
Definition: monster.hh:117
std::vector< clusterdata > cData
Definition: netcluster.hh:391
#define MIN_SKYRES_ANGLE
Definition: cwb2G.cc:29
double Phi2
Definition: config.hh:280
netevent * netburst
livetime object
Definition: cwb.hh:270
char skyMaskCCFile[1024]
Definition: config.hh:309
virtual wavearray< double > select(char *, double)
Definition: netcluster.cc:263
DataType_t * data
Definition: wavearray.hh:319
double Tgap
Definition: config.hh:134
netcluster wc
wavearray< double > FREQ
Definition: detector.hh:373
Long_t id
Definition: ced.hh:44
void setMRAcatalog(char *fn)
Definition: network.hh:315
double factors[FACTORS_MAX]
Definition: config.hh:202
void DataConditioning(int ifactor)
Definition: cwb2G.cc:450
int BetaOrder
Definition: monster.hh:118
wavearray< double > * hot[NIFO_MAX]
wavelet pointers: pwdm[0] - l_high, wdm[nRES-1] l_low
Definition: cwb2G.hh:81
cout<< "live time after cat 2 : "<< detSegs_ctime<< endl;if(detSegs_ctime< segTHR) {cout<< "job segment live time after cat2 < "<< segTHR<< " sec, job terminated !!!"<< endl;exit(1);} double Tb=detSegs[0].start;double Te=detSegs[0].stop;double dT=Te-Tb;char file[512], tdf00[512], tdf90[512], buFFer[1024];int rnID=int(gRandom->Rndm(13) *1.e9);if(simulation) { i=NET.readMDClog(injectionList, double(long(Tb)) -mdcShift);printf("GPS: %16.6f saved, injections: %d\", double(long(Tb)), i);frTB[nIFO].shiftBurstMDCLog(NET.mdcList, ifos, mdcShift);for(int i=0;i< NET.mdcTime.size();i++) NET.mdcTime[i]+=mdcShift;vector< waveSegment > mdcSegs(NET.mdcTime.size());for(int k=0;k< NET.mdcTime.size();k++) {mdcSegs[k].start=NET.mdcTime[k]-gap;mdcSegs[k].stop=NET.mdcTime[k]+gap;} vector< waveSegment > mdcSegs_dq2=slagTB.mergeSegLists(detSegs_dq2, mdcSegs);double mdcSegs_ctime=slagTB.getTimeSegList(mdcSegs_dq2);cout<< "live time in zero lag after cat2+inj : "<< mdcSegs_ctime<< endl;if(mdcSegs_ctime==0) {cout<< "job segment with zero cat2+inj live time in zero lag, job terminated !!!"<< endl;exit(1);} } if(dump_infos_and_exit) exit(0);if(mask >0.) NET.setSkyMask(mask, skyMaskFile);for(i=0;i< nIFO;i++) { frTB[i].readFrames(FRF[i], channelNamesRaw[i], x);x.start(x.start()+dataShift[i]);x.start(x.start() -segLen *(segID[i]-segID[0]));if(singleDetector) TB.resampleToPowerOfTwo(x);sprintf(file,"%s/%s_%d_%s_%d_%d.dat", nodedir, ifo[i], int(Tb), data_label, runID, rnID);if(dump_sensitivity_and_exit) { sprintf(file,"%s/sensitivity_%s_%d_%s_job%d.txt", dump_dir, ifo[i], int(Tb), data_label, runID);cout<< endl<< "Dump Sensitivity : "<< file<< endl<< endl;TB.makeSpectrum(file, x);continue;} if(dcCal[i]>0.) x *=dcCal[i];if(fResample >0) { x.FFT(1);x.resize(fResample/x.rate() *x.size());x.FFT(-1);x.rate(fResample);} pTF[i]=pD[i]-> getTFmap()
int nIFO
Toolbox.
Definition: cwb.hh:244
char wdmXTalk[1024]
Definition: config.hh:215
std::vector< SSeries< double > > vSS
Definition: detector.hh:352
double lagStep
Definition: config.hh:169
size_t read(const char *)
Definition: netcluster.cc:3115
double ReadData(double mdcShift, int ifactor)
Definition: cwb2G.cc:202
#define REGRESSION_SOLVE_EIGEN_NUM
Definition: cwb2G.cc:35
char fName[256]
#define REGRESSION_SOLVE_REGULATOR
Definition: cwb2G.cc:36
virtual void resize(unsigned int)
Definition: wavearray.cc:463
double TFgap
Definition: config.hh:138
bool MRA
Definition: network.hh:599
double bpp
Definition: config.hh:133
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
Definition: wseries.cc:291
wavearray< double > y
Definition: Test10.C:31
double Theta1
Definition: config.hh:277
void putLayer(wavearray< DataType_t > &, double n)
param: n - layer number
Definition: wseries.cc:219
size_t psize(int k=2)
Definition: netcluster.hh:163
size_t setsim(WSeries< double > &, std::vector< double > *, double=5., double=8., bool saveWF=false)
Definition: detector.cc:1377
size_t inRate
Definition: config.hh:132
#define WDM_BETAORDER
Definition: cwb2G.cc:44
double frequency(int l)
Definition: wseries.cc:117
double segMLS
Definition: config.hh:162
wavearray< double > ISNR
Definition: detector.hh:372
int maxLayer()
Definition: wseries.hh:139
int m_H
Definition: Wavelet.hh:121
int levelR
Definition: config.hh:152
int nRES
Definition: cwb2G.hh:77
int wfsave()
Definition: detector.hh:316
bool scPlugin
Definition: config.hh:368
TString ifos[60]
double angle
Definition: config.hh:276
void setSkyMaps(double, double=0., double=180., double=0., double=360.)
param: sky map granularity step, degrees param: theta begin, degrees param: theta end...
Definition: network.cc:2687
int precision
Definition: monster.hh:119
virtual WSeries< double > white(double, int, double=0., double=0.)
what it does: each wavelet layer is devided into k intervals.
Definition: wseries.cc:1146
void clear()
Definition: netcluster.hh:427
bool dump
Definition: config.hh:295