Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CWB_Plugin_supercluster_bench.C
Go to the documentation of this file.
1 #define XIFO 4
2 
3 #pragma GCC system_header
4 
5 #include "cwb.hh"
6 #include "cwb2G.hh"
7 #include "config.hh"
8 #include "network.hh"
9 #include "wavearray.hh"
10 #include "TString.h"
11 #include "TObjArray.h"
12 #include "TObjString.h"
13 #include "TRandom.h"
14 #include "TComplex.h"
15 
16 //!SUPERCLUSTER
17 
18 long subNetCut(network* net, int lag, float snc, TH2F* hist);
19 inline int _sse_MRA_ps(network* net, float* amp, float* AMP, float Eo, int K);
20 void PrintElapsedTime(int job_elapsed_time, double cpu_time, TString info);
21 
22 #define USE_LOCAL_SUBNETCUT // comment to use the builtin implementation of subNetCut
23 
24 void
25 CWB_Plugin(TFile* jfile, CWB::config* cfg, network* net, WSeries<double>* x, TString ifo, int type) {
26 
27 // This plugin implements the standard supercluster stage & use a local implementation of subNetCut (only 2G)
28 
29  cout << endl;
30  cout << "-----> CWB_Plugin_supercluster_bench.C" << endl;
31  cout << "ifo " << ifo.Data() << endl;
32  cout << "type " << type << endl;
33  cout << endl;
34 
35  if(type==CWB_PLUGIN_CONFIG) {
36  cfg->scPlugin=true; // disable built-in supercluster function
37  }
38 
39  if(type==CWB_PLUGIN_ISUPERCLUSTER) {
40 
41  cout << "type==CWB_PLUGIN_ISUPERCLUSTER" << endl;
42 
43  TStopwatch bench;
44  bench.Stop();
45 
46  // import ifile
47  void* gIFILE=NULL; IMPORT(void*,gIFILE)
48  cout << "-----> CWB_Plugin_wavegraph.C -> " << " gIFILE : " << gIFILE << endl;
49  TFile* ifile = (TFile*)gIFILE;
50 
51  // import ifactor
52  int gIFACTOR=-1; IMPORT(int,gIFACTOR)
53  cout << "-----> CWB_Plugin_wavegraph.C -> " << " gIFACTOR : " << gIFACTOR << endl;
54  int ifactor = gIFACTOR;
55 
56  int nIFO = net->ifoListSize(); // number of detectors
57  int rateANA=cfg->inRate>>cfg->levelR;
58  int nRES = net->wdmListSize(); // number of resolution levels
59  int lags = net->getifo(0)->lagShift.size();
60 
61  wavearray<double>* hot[NIFO_MAX]; // temporary time series
62  for(int i=0; i<nIFO; i++) hot[i] = net->getifo(i)->getHoT();
63 
64  int nevt = 0;
65  int nnn = 0;
66  int mmm = 0;
67  size_t count = 0;
68  netcluster wc;
69  netcluster* pwc;
70 
71  for(int j=0; j<(int)lags; j++) {
72 
73  int cycle = cfg->simulation ? ifactor : Long_t(net->wc_List[j].shift);
74 
75  // read cluster metadata
76  if(ifile!=NULL) wc.read(ifile,"coherence","clusters",0,cycle);
77  else wc.read(jfile,"coherence","clusters",0,cycle);
78  // read clusters from temporary job file, loop over TF resolutions
79  if(ifile!=NULL) {
80  for(int i=nRES-1; i>=0; i--) // reverse loop is faster loading cluster (?)
81  wc.read(ifile,"coherence","clusters",-1,cycle,rateANA>>(i+cfg->l_low));
82  } else {
83  for(int i=nRES-1; i>=0; i--) // reverse loop is faster loading cluster (?)
84  wc.read(jfile,"coherence","clusters",-1,cycle,rateANA>>(i+cfg->l_low));
85  }
86  if(!cfg->simulation) cout<<"process lag " <<cycle<<" ..."<<endl;
87  cout<<"loaded clusters|pixels: "<<wc.csize()<<"|"<<wc.size()<<endl;
88 
89  // supercluster analysis
90  wc.supercluster('L',net->e2or,cfg->TFgap,false); //likehood2G
91  cout<<"super clusters|pixels: "<<wc.esize(0)<<"|"<<wc.psize(0)<<endl;
92 
93  // release all pixels
94  pwc = net->getwc(j);
95  pwc->cpf(wc, false);
96 
97  net->setDelayIndex(hot[0]->rate());
98  pwc->setcore(false);
99 
100  // apply cuts
101  int psel = 0;
102  while(1) {
103  count = pwc->loadTDampSSE(*net, 'a', cfg->BATCH, cfg->LOUD);
104  bench.Continue();
105 #ifdef USE_LOCAL_SUBNETCUT
106  psel += subNetCut(net,(int)j,cfg->subnet,NULL);
107 #else
108  psel += net->subNetCut((int)j,cfg->subnet,NULL);
109 #endif
110  bench.Stop();
111  PrintElapsedTime(bench.RealTime(),bench.CpuTime(),"subNetCut : Processing Time - ");
112  int ptot = pwc->psize(1)+pwc->psize(-1);
113  double pfrac = ptot>0 ? double(psel)/double(ptot) : 0.;
114  cout<<"selected pixels: "<<psel<<", fraction: "<<pfrac<< endl;
115  if(count<10000) break;
116  }
117 
118  pwc->defragment(cfg->Tgap,cfg->Fgap); // SK added defragmentation
119 
120  nevt = net->events();
121  nnn += pwc->psize(-1);
122  mmm += pwc->psize(1)+pwc->psize(-1);
123 
124  if(mmm) cout<<"events in the buffer: "<<net->events()<<"|"<<nnn<<"|"<<nnn/double(mmm)<<"\n";
125  else cout<<"events in the buffer: "<<net->events()<<"\n";
126 
127  // store cluster into temporary job file [NEWSS]
128  pwc->write(jfile,"supercluster","clusters",0,cycle);
129  pwc->write(jfile,"supercluster","clusters",-1,cycle);
130  cout<<cycle<<"|"<<pwc->csize()<<"|"<<pwc->size()<<" ";cout.flush();
131 
132  pwc->clear();
133  cout<<endl;
134  }
135  }
136 
137  return;
138 }
139 
140 void
141 PrintElapsedTime(int job_elapsed_time, double cpu_time, TString info) {
142 //
143 // convert job_elapsed_time to (hh:mm:ss) format and print it
144 //
145 // job_elapsed_time : time (seconds)
146 //
147 // info : info string added to (hh:mm:ss)
148 //
149 
150  int job_elapsed_hour = int(job_elapsed_time/3600);
151  int job_elapsed_min = int((job_elapsed_time-3600*job_elapsed_hour)/60);
152  int job_elapsed_sec = int(job_elapsed_time-3600*job_elapsed_hour-60*job_elapsed_min);
153  char buf[1024];
154  sprintf(buf,"%s %02d:%02d:%02d (hh:mm:ss) : cpu time : %f (sec)\n",info.Data(),job_elapsed_hour,job_elapsed_min,job_elapsed_sec,cpu_time);
155  cout << buf;
156 
157  return;
158 }
159 
160 long subNetCut(network* net, int lag, float snc, TH2F* hist)
161 {
162 // sub-network cut with dsp regulator
163 // lag: lag index
164 // snc: sub network threshold, if snc<0 use weak constraint
165 // hist: diagnostic histogram
166 // return number of processed pixels
167 
168  if(!net->wc_List[lag].size()) return 0;
169 
170  size_t nIFO = net->ifoList.size();
171 
172  if(nIFO>NIFO) {
173  cout<<"network::subNetCut(): invalid network.\n";
174  exit(0);
175  }
176 
177  float En = 2*net->acor*net->acor*nIFO; // network energy threshold in the sky loop
178  float Es = 2*net->e2or; // subnet energy threshold in the sky loop
179  float TH = fabs(snc); // sub network threshold
180 
181  __m128 _En = _mm_set1_ps(En);
182  __m128 _Es = _mm_set1_ps(Es);
183  __m128 _oo = _mm_set1_ps(1.e-12);
184  __m128 _0 = _mm_set1_ps(0.);
185  __m128 _05 = _mm_set1_ps(0.5);
186  __m128 _1 = _mm_set1_ps(1.);
187  __m128* _pe[NIFO];
188 
189  int f_ = NIFO/4;
190  int l,lm,Vm;
191  float Lm,Em,Am,Lo,Eo,Co,Lr,Er,ee,em,To;
192  float cc,aa,AA,rHo,stat,Ls,Ln,EE;
193 
194  size_t i,j,k,m,V,V4,id,K,M;
195  int Lsky = int(net->index.size()); // total number of source locations
196  short* mm = net->skyMask.data;
197 
198  float vvv[NIFO];
199  float* v00[NIFO];
200  float* v90[NIFO];
201  float* pe[NIFO];
202  float* pa[NIFO];
203  float* pA[NIFO];
204  short* ml[NIFO];
205  double* FP[NIFO];
206  double* FX[NIFO];
207  double xx[NIFO];
208 
209  for(i=0; i<NIFO; i++) {
210  if(i<nIFO) {
211  ml[i] = net->getifo(i)->index.data;
212  FP[i] = net->getifo(i)->fp.data;
213  FX[i] = net->getifo(i)->fx.data;
214  }
215  else {
216  ml[i] = net->getifo(0)->index.data;
217  FP[i] = net->getifo(0)->fp.data;
218  FX[i] = net->getifo(0)->fx.data;
219  }
220  }
221 
222  // allocate buffers
223  std::vector<int> pI; // buffer for pixel IDs
224  wavearray<double> cid; // buffers for cluster ID
225  netpixel* pix;
226  std::vector<int>* vint;
227  netcluster* pwc = &net->wc_List[lag];
228 
229  size_t count = 0;
230  size_t tsize = 0;
231 
232 //+++++++++++++++++++++++++++++++++++++++
233 // loop over clusters
234 //+++++++++++++++++++++++++++++++++++++++
235 
236  cid = pwc->get((char*)"ID", 0,'S',0); // get cluster ID
237 
238  K = cid.size();
239  for(k=0; k<K; k++) { // loop over clusters
240  id = size_t(cid.data[k]+0.1);
241  if(pwc->sCuts[id-1] != -2) continue; // skip rejected/processed clusters
242  vint = &(pwc->cList[id-1]); // pixel list
243  V = vint->size(); // pixel list size
244  if(!V) continue;
245 
246  pI = net->wdmMRA.getXTalk(pwc, id);
247 
248  V = pI.size(); // number of loaded pixels
249  if(!V) continue;
250 
251  pix = pwc->getPixel(id,pI[0]);
252  tsize = pix->tdAmp[0].size();
253  if(!tsize || tsize&1) { // tsize%1 = 1/0 = power/amplitude
254  cout<<"network::subNetCut() error: wrong pixel TD data\n";
255  exit(1);
256  }
257  tsize /= 2;
258  V4 = V + (V%4 ? 4 - V%4 : 0);
259 
260  //cout<<En<<" "<<Es<<" "<<lag<<" "<<id<<" "<<V4<<" "<<" "<<tsize<<endl;
261 
262  std::vector<wavearray<float> > vtd; // vectors of TD amplitudes
263  std::vector<wavearray<float> > vTD; // vectors of TD amplitudes
264  std::vector<wavearray<float> > eTD; // vectors of TD energies
265 
266  wavearray<float> tmp(tsize*V4); tmp=0; // aligned array for TD amplitudes
267  wavearray<float> fp(NIFO*V4); fp=0; // aligned array for + antenna pattern
268  wavearray<float> fx(NIFO*V4); fx=0; // aligned array for x antenna pattern
269  wavearray<float> nr(NIFO*V4); nr=0; // aligned array for inverse rms
270  wavearray<float> Fp(NIFO*V4); Fp=0; // aligned array for pattern
271  wavearray<float> Fx(NIFO*V4); Fx=0; // aligned array for patterns
272  wavearray<float> am(NIFO*V4); am=0; // aligned array for TD amplitudes
273  wavearray<float> AM(NIFO*V4); AM=0; // aligned array for TD amplitudes
274  wavearray<float> bb(NIFO*V4); bb=0; // temporary array for MRA amplitudes
275  wavearray<float> BB(NIFO*V4); BB=0; // temporary array for MRA amplitudes
276  wavearray<float> xi(NIFO*V4); xi=0; // 00 array for reconctructed responses
277  wavearray<float> XI(NIFO*V4); XI=0; // 90 array for reconstructed responses
278  wavearray<float> ww(NIFO*V4); ww=0; // 00 array for phase-shifted data vectors
279  wavearray<float> WW(NIFO*V4); WW=0; // 90 array for phase-shifted data vectors
280  wavearray<float> u4(NIFO*4); u4=0; // temp array
281  wavearray<float> U4(NIFO*4); U4=0; // temp array
282 
283  __m128* _Fp = (__m128*) Fp.data;
284  __m128* _Fx = (__m128*) Fx.data;
285  __m128* _am = (__m128*) am.data;
286  __m128* _AM = (__m128*) AM.data;
287  __m128* _xi = (__m128*) xi.data;
288  __m128* _XI = (__m128*) XI.data;
289  __m128* _fp = (__m128*) fp.data;
290  __m128* _fx = (__m128*) fx.data;
291  __m128* _nr = (__m128*) nr.data;
292  __m128* _ww = (__m128*) ww.data;
293  __m128* _WW = (__m128*) WW.data;
294  __m128* _bb = (__m128*) bb.data;
295  __m128* _BB = (__m128*) BB.data;
296 
297  for(i=0; i<NIFO; i++) {
298  vtd.push_back(tmp); // array of aligned energy vectors
299  vTD.push_back(tmp); // array of aligned energy vectors
300  eTD.push_back(tmp); // array of aligned energy vectors
301  }
302 
303  for(i=0; i<NIFO; i++) { // set up zero deley pointers
304  pa[i] = vtd[i].data + (tsize/2)*V4;
305  pA[i] = vTD[i].data + (tsize/2)*V4;
306  pe[i] = eTD[i].data + (tsize/2)*V4;
307  }
308 
309  net->a_00.resize(NIFO*V4); net->a_00=0.;
310  net->a_90.resize(NIFO*V4); net->a_90=0.;
311  net->rNRG.resize(V4); net->rNRG=0.;
312  net->pNRG.resize(V4); net->pNRG=0.;
313 
314  __m128* _aa = (__m128*) net->a_00.data; // set pointer to 00 array
315  __m128* _AA = (__m128*) net->a_90.data; // set pointer to 90 array
316 
317  net->pList.clear();
318  for(j=0; j<V; j++) { // loop over selected pixels
319  pix = pwc->getPixel(id,pI[j]); // get pixel pointer
320  net->pList.push_back(pix); // store pixel pointers for MRA
321 
322  double rms = 0.;
323  for(i=0; i<nIFO; i++) {
324  xx[i] = 1./pix->data[i].noiserms;
325  rms += xx[i]*xx[i]; // total inverse variance
326  }
327 
328  for(i=0; i<nIFO; i++) {
329  nr.data[j*NIFO+i]=(float)xx[i]/sqrt(rms); // normalized 1/rms
330  for(l=0; l<tsize; l++) {
331  aa = pix->tdAmp[i].data[l]; // copy TD 00 data
332  AA = pix->tdAmp[i].data[l+tsize]; // copy TD 90 data
333  vtd[i].data[l*V4+j] = aa; // copy 00 data
334  vTD[i].data[l*V4+j] = AA; // copy 90 data
335  eTD[i].data[l*V4+j] = aa*aa+AA*AA; // copy power
336  }
337  }
338  }
339 
340 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
341 // first sky loop
342 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
343 
344  int lb = 0;
345  int le = Lsky-1;
346  bool mra = false;
347  double suball=0;
348  double submra=0;
349 
350  stat=Lm=Em=Am=EE=0.; lm=Vm= -1;
351 
352  skyloop:
353 
354  for(l=lb; l<=le; l++) { // loop over sky locations
355  if(!mm[l] || l<0) continue; // skip delay configurations
356 
357  _sse_point_ps(_pe, pe, ml, int(l), (int)V4); // point _pe to energy vectors
358 
359  __m128 _msk;
360  __m128 _E_o = _mm_setzero_ps(); // total network energy
361  __m128 _E_n = _mm_setzero_ps(); // network energy above the threshold
362  __m128 _E_s = _mm_setzero_ps(); // subnet energy above the threshold
363  __m128 _M_m = _mm_setzero_ps(); // # of pixels above threshold
364  __m128* _rE = (__m128*) net->rNRG.data; // m128 pointer to energy array
365  __m128* _pE = (__m128*) net->pNRG.data; // m128 pointer to energy array
366 
367  for(j=0; j<V4; j+=4) { // loop over selected pixels
368  *_rE = _sse_sum_ps(_pe); // get pixel energy
369  _msk = _mm_and_ps(_1,_mm_cmpge_ps(*_rE,_En)); // E>En 0/1 mask
370  _M_m = _mm_add_ps(_M_m,_msk); // count pixels above threshold
371  *_pE = _mm_mul_ps(*_rE,_msk); // zero sub-threshold pixels
372  _E_o = _mm_add_ps(_E_o,*_pE); // network energy
373  _sse_minSNE_ps(_rE,_pe,_pE); // subnetwork energy with _pe increment
374  _E_s = _mm_add_ps(_E_s,*_pE); // subnetwork energy
375  _msk = _mm_and_ps(_1,_mm_cmpge_ps(*_pE++,_Es)); // subnet energy > Es 0/1 mask
376  _E_n = _mm_add_ps(_E_n,_mm_mul_ps(*_rE++,_msk)); // network energy
377  }
378 
379  _mm_storeu_ps(vvv,_E_n);
380  Ln = vvv[0]+vvv[1]+vvv[2]+vvv[3]; // network energy above subnet threshold
381  _mm_storeu_ps(vvv,_E_o);
382  Eo = vvv[0]+vvv[1]+vvv[2]+vvv[3]+0.01; // total network energy
383  _mm_storeu_ps(vvv,_E_s);
384  Ls = vvv[0]+vvv[1]+vvv[2]+vvv[3]; // subnetwork energy
385  _mm_storeu_ps(vvv,_M_m);
386  m = 2*(vvv[0]+vvv[1]+vvv[2]+vvv[3])+0.01; // pixels above threshold
387 
388  aa = Ls*Ln/(Eo-Ls);
389  if((aa-m)/(aa+m)<0.33) continue;
390 
391  net->pnt_(v00, pa, ml, (int)l, (int)V4); // pointers to first pixel 00 data
392  net->pnt_(v90, pA, ml, (int)l, (int)V4); // pointers to first pixel 90 data
393  float* pfp = fp.data; // set pointer to fp
394  float* pfx = fx.data; // set pointer tp fx
395  float* p00 = net->a_00.data; // set pointer for 00 array
396  float* p90 = net->a_90.data; // set pointer for 90 array
397 
398  m = 0;
399  for(j=0; j<V; j++) {
400  int jf = j*f_; // source sse pointer increment
401  net->cpp_(p00,v00); net->cpp_(p90,v90); // copy amplitudes with target increment
402  net->cpf_(pfp,FP,l); net->cpf_(pfx,FX,l); // copy antenna with target increment
403  _sse_zero_ps(_xi+jf); // zero MRA amplitudes
404  _sse_zero_ps(_XI+jf); // zero MRA amplitudes
405  _sse_cpf_ps(_am+jf,_aa+jf); // duplicate 00
406  _sse_cpf_ps(_AM+jf,_AA+jf); // duplicate 90
407  if(net->rNRG.data[j]>En) m++; // count superthreshold pixels
408  }
409 
410  __m128* _pp = (__m128*) am.data; // point to multi-res amplitudes
411  __m128* _PP = (__m128*) AM.data; // point to multi-res amplitudes
412 
413  if(mra) { // do MRA
414  _sse_MRA_ps(net,xi.data,XI.data,En,m); // get principal components
415  _pp = (__m128*) xi.data; // point to PC amplitudes
416  _PP = (__m128*) XI.data; // point to PC amplitudes
417  }
418 
419  m = 0; Ls=Ln=Eo=0;
420  for(j=0; j<V; j++) {
421  int jf = j*f_; // source sse pointer increment
422  int mf = m*f_; // target sse pointer increment
423  _sse_zero_ps(_bb+jf); // reset array for MRA amplitudes
424  _sse_zero_ps(_BB+jf); // reset array for MRA amplitudes
425  ee = _sse_abs_ps(_pp+jf,_PP+jf); // total pixel energy
426  if(ee<En) continue;
427  _sse_cpf_ps(_bb+mf,_pp+jf); // copy 00 amplitude/PC
428  _sse_cpf_ps(_BB+mf,_PP+jf); // copy 90 amplitude/PC
429  _sse_cpf_ps(_Fp+mf,_fp+jf); // copy F+
430  _sse_cpf_ps(_Fx+mf,_fx+jf); // copy Fx
431  _sse_mul_ps(_Fp+mf,_nr+jf); // normalize f+ by rms
432  _sse_mul_ps(_Fx+mf,_nr+jf); // normalize fx by rms
433  m++;
434  em = _sse_maxE_ps(_pp+jf,_PP+jf); // dominant pixel energy
435  Ls += ee-em; Eo += ee; // subnetwork energy, network energy
436  if(ee-em>Es) Ln += ee; // network energy above subnet threshold
437  }
438 
439  size_t m4 = m + (m%4 ? 4 - m%4 : 0);
440  _E_n = _mm_setzero_ps(); // + likelihood
441 
442  for(j=0; j<m4; j+=4) {
443  int jf = j*f_;
444  _sse_dpf4_ps(_Fp+jf,_Fx+jf,_fp+jf,_fx+jf); // go to DPF
445  _E_s = _sse_like4_ps(_fp+jf,_fx+jf,_bb+jf,_BB+jf); // std likelihood
446  _E_n = _mm_add_ps(_E_n,_E_s); // total likelihood
447  }
448  _mm_storeu_ps(vvv,_E_n);
449 
450  Lo = vvv[0]+vvv[1]+vvv[2]+vvv[3];
451  AA = aa/(fabs(aa)+fabs(Eo-Lo)+2*m*(Eo-Ln)/Eo); // subnet stat with threshold
452  ee = Ls*Eo/(Eo-Ls);
453  em = fabs(Eo-Lo)+2*m; // suball NULL
454  ee = ee/(ee+em); // subnet stat without threshold
455  aa = (aa-m)/(aa+m);
456 
457  if(AA>stat && !mra) {
458  stat=AA; Lm=Lo; Em=Eo; Am=aa; lm=l; Vm=m; suball=ee; EE=em;
459  }
460  }
461 
462  if(!mra && lm>=0) {mra=true; le=lb=lm; goto skyloop;} // get MRA principle components
463 
464  pwc->sCuts[id-1] = -1;
465  pwc->cData[id-1].likenet = Lm;
466  pwc->cData[id-1].energy = Em;
467  pwc->cData[id-1].theta = net->nLikelihood.getTheta(lm);
468  pwc->cData[id-1].phi = net->nLikelihood.getPhi(lm);
469  pwc->cData[id-1].skyIndex = lm;
470 
471  rHo = 0.;
472  if(mra) {
473  submra = Ls*Eo/(Eo-Ls); // MRA subnet statistic
474  submra/= fabs(submra)+fabs(Eo-Lo)+2*(m+6); // MRA subnet coefficient
475  To = 0;
476  pwc->p_Ind[id-1].push_back(lm);
477  for(j=0; j<vint->size(); j++) {
478  pix = pwc->getPixel(id,j);
479  pix->theta = net->nLikelihood.getTheta(lm);
480  pix->phi = net->nLikelihood.getPhi(lm);
481  To += pix->time/pix->rate/pix->layers;
482  if(j==0&&mra) pix->ellipticity = submra; // subnet MRA propagated to L-stage
483  if(j==0&&mra) pix->polarisation = fabs(Eo-Lo)+2*(m+6); // submra NULL propagated to L-stage
484  if(j==1&&mra) pix->ellipticity = suball; // subnet all-sky propagated to L-stage
485  if(j==1&&mra) pix->polarisation = EE; // suball NULL propagated to L-stage
486  }
487  To /= vint->size();
488  rHo = sqrt(Lo*Lo/(Eo+2*m)/nIFO); // estimator of coherent amplitude
489  }
490 
491  if(hist && rHo>net->netRHO)
492  for(j=0;j<vint->size();j++) hist->Fill(suball,submra);
493 
494  if(fmin(suball,submra)>TH && rHo>net->netRHO) {
495  count += vint->size();
496  if(hist) {
497  printf("lag|id %3d|%3d rho=%5.2f To=%5.1f stat: %5.3f|%5.3f|%5.3f ",
498  int(lag),int(id),rHo,To,suball,submra,stat);
499  printf("E: %6.1f|%6.1f L: %6.1f|%6.1f|%6.1f pix: %4d|%4d|%3d|%2d \n",
500  Em,Eo,Lm,Lo,Ls,int(vint->size()),int(V),Vm,int(m));
501  }
502  }
503  else pwc->sCuts[id-1]=1;
504 
505 // clean time delay data
506 
507  V = vint->size();
508  for(j=0; j<V; j++) { // loop over pixels
509  pix = pwc->getPixel(id,j);
510  pix->core = true;
511  if(pix->tdAmp.size()) pix->clean();
512  }
513  } // end of loop over clusters
514  return count;
515 }
516 
517 inline int _sse_MRA_ps(network* net, float* amp, float* AMP, float Eo, int K) {
518 // fast multi-resolution analysis inside sky loop
519 // select max E pixel and either scale or skip it based on the value of residual
520 // pointer to 00 phase amplitude of monster pixels
521 // pointer to 90 phase amplitude of monster pixels
522 // Eo - energy threshold
523 // K - number of principle components to extract
524 // returns number of MRA pixels
525  int j,n,mm;
526  int k = 0;
527  int m = 0;
528  int f = NIFO/4;
529  int V = (int)net->rNRG.size();
530  float* ee = net->rNRG.data; // residual energy
531  float* pp = net->pNRG.data; // residual energy
532  float EE = 0.; // extracted energy
533  float E;
534  float mam[NIFO];
535  float mAM[NIFO];
536  net->pNRG=-1;
537  for(j=0; j<V; ++j) if(ee[j]>Eo) pp[j]=0;
538 
539  __m128* _m00 = (__m128*) mam;
540  __m128* _m90 = (__m128*) mAM;
541  __m128* _amp = (__m128*) amp;
542  __m128* _AMP = (__m128*) AMP;
543  __m128* _a00 = (__m128*) net->a_00.data;
544  __m128* _a90 = (__m128*) net->a_90.data;
545 
546  while(k<K){
547 
548  for(j=0; j<V; ++j) if(ee[j]>ee[m]) m=j; // find max pixel
549  if(ee[m]<=Eo) break; mm = m*f;
550 
551  //cout<<" V= "<<V<<" m="<<m<<" ee[m]="<<ee[m];
552 
553  E = _sse_abs_ps(_a00+mm,_a90+mm); EE += E; // get PC energy
554  int J = net->wdmMRA.getXTalk(m)->size()/7;
555  float* c = net->wdmMRA.getXTalk(m)->data; // c1*c2+c3*c4=c1*c3+c2*c4=0
556 
557  if(E/EE < 0.01) break; // ignore small PC
558 
559  _sse_cpf_ps(mam,_a00+mm); // store a00 for max pixel
560  _sse_cpf_ps(mAM,_a90+mm); // store a90 for max pixel
561  _sse_add_ps(_amp+mm,_m00); // update 00 PC
562  _sse_add_ps(_AMP+mm,_m90); // update 90 PC
563 
564  for(j=0; j<J; j++) {
565  n = int(c[0]+0.1);
566  if(ee[n]>Eo) {
567  ee[n] = _sse_rotsub_ps(_m00,c[1],_m90,c[2],_a00+n*f); // subtract PC from a00
568  ee[n]+= _sse_rotsub_ps(_m00,c[3],_m90,c[4],_a90+n*f); // subtract PC from a90
569  }
570  c += 7;
571  }
572  //cout<<" "<<ee[m]<<" "<<k<<" "<<E<<" "<<EE<<" "<<endl;
573  pp[m] = _sse_abs_ps(_amp+mm,_AMP+mm); // store PC energy
574  k++;
575  }
576 /*
577  cout<<"EE="<<EE<<endl;
578  EE = 0;
579  for(j=0; j<V; ++j) {
580  if(pp[j]>=0) EE += ee[j];
581  if(pp[j]>=0.) cout<<j<<"|"<<pp[j]<<"|"<<ee[j]<<" "; // find max pixel
582  }
583  cout<<"EE="<<EE<<endl;
584 */
585  return k;
586 }
587 
network::wdmMRA
monster wdmMRA
list of pixel pointers for MRA
Definition: network.hh:633
network::getifo
detector * getifo(size_t n)
param: detector index
Definition: network.hh:418
cfg
CWB::config * cfg
Definition: TestCWB_Plugin.C:5
netcluster::defragment
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:1275
_sse_abs_ps
static float _sse_abs_ps(__m128 *_a)
Definition: watsse.hh:119
wavearray::size
virtual size_t size() const
Definition: wavearray.hh:127
NIFO
#define NIFO
Definition: wat.hh:56
netcluster::write
size_t write(const char *file, int app=0)
Definition: netcluster.cc:2989
netpixel::phi
float phi
Definition: netpixel.hh:99
job_elapsed_min
int job_elapsed_min
Definition: cwb_net.C:720
network::wc_List
std::vector< netcluster > wc_List
Definition: network.hh:592
cwb_online.f
tuple f
Definition: cwb_online.py:91
printf
printf("total live time: non-zero lags = %10.1f \n", liveTot)
network::cpp_
static void cpp_(float *&, float **)
Definition: network.hh:1100
netpixel::tdAmp
std::vector< wavearray< float > > tdAmp
Definition: netpixel.hh:105
n
int n
Definition: cwb_net.C:10
_sse_zero_ps
static void _sse_zero_ps(__m128 *_p)
Definition: watsse.hh:26
network::a_90
wavearray< float > a_90
buffer for cluster sky 00 amplitude
Definition: network.hh:635
TString
TString("c")
Definition: cwb_report_skymap.C:111
detector::fx
wavearray< double > fx
Definition: detector.hh:346
count
int count
Definition: compare_bkg.C:373
x
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1)*itime+ifreq;double time=itime *dT;double freq=(ifreq >0)?ifreq *dF:dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
Definition: BaseFunctionWDM.C:51
netpixel::data
std::vector< pixdata > data
Definition: netpixel.hh:104
netcluster::get
wavearray< double > get(char *name, size_t index=0, char atype='R', int type=1, bool=true)
param: string with parameter name param: index in the amplitude array, which define detector param: c...
Definition: netcluster.cc:2188
pix
netpixel pix(nifo)
network::pList
std::vector< netpixel * > pList
Definition: network.hh:632
pwc
netcluster * pwc
Definition: cwb_job_obj.C:20
detector::index
wavearray< short > index
Definition: detector.hh:350
_sse_rotsub_ps
static float _sse_rotsub_ps(__m128 *_u, float c, __m128 *_v, float s, __m128 *_a)
Definition: watsse.hh:381
skymap::getTheta
double getTheta(size_t i)
Definition: skymap.hh:206
netcluster::setcore
size_t setcore(bool core, int id=0)
Definition: netcluster.cc:224
_sse_cpf_ps
static void _sse_cpf_ps(float *a, __m128 *_p)
Definition: watsse.hh:289
nevt
int nevt
Definition: RootSkyMap2Fits.C:12
M
#define M
Definition: UniqSLagsList.C:3
netcluster::supercluster
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:789
netpixel::layers
size_t layers
Definition: netpixel.hh:94
m
int m
Definition: cwb_net.C:10
netcluster::cList
std::vector< vector_int > cList
Definition: netcluster.hh:379
j
int j
Definition: cwb_net.C:10
i
i drho i
Definition: cwb_epparameters.C:85
network::subNetCut
long subNetCut(int lag, float subnet=0.6, float subcut=0.33, TH2F *hist=NULL)
Definition: network.cc:983
network::ifoList
std::vector< detector * > ifoList
Definition: network.hh:590
subNetCut
long subNetCut(network *net, int lag, float snc, TH2F *hist)
SUPERCLUSTER.
Definition: CWB_Plugin_supercluster_bench.C:160
network::wdmListSize
size_t wdmListSize()
Definition: network.hh:428
netpixel::core
bool core
Definition: netpixel.hh:102
ifo
char ifo[NIFO_MAX][8]
Definition: cwb1G_parameters.C:17
net
network ** net
NOISE_MDC_SIMULATION.
Definition: CWB_Plugin_HEN_BKG_Config.C:4
network::ifoListSize
size_t ifoListSize()
Definition: network.hh:413
network::nLikelihood
skymap nLikelihood
Definition: network.hh:604
nIFO
#define nIFO
Definition: DrawSensitivitiesS5.C:11
network::pNRG
wavearray< float > pNRG
buffers for cluster residual energy
Definition: network.hh:637
netcluster::esize
size_t esize(int k=2)
Definition: netcluster.hh:135
jfile
jfile
Definition: cwb_job_obj.C:25
xx
wavearray< double > xx
Definition: TestFrame1.C:11
_sse_MRA_ps
int _sse_MRA_ps(network *net, float *amp, float *AMP, float Eo, int K)
Definition: CWB_Plugin_supercluster_bench.C:517
CWB::config::simulation
int simulation
Definition: config.hh:181
CWB::config::Fgap
double Fgap
Definition: config.hh:118
hot
wavearray< double > hot[2]
Definition: ReadCStrainFromJobFile.C:15
IMPORT
#define IMPORT(TYPE, VAR)
Definition: cwb.hh:51
monster::getXTalk
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
int
i() int(T_cor *100))
netpixel::polarisation
float polarisation
Definition: netpixel.hh:101
detector::getHoT
wavearray< double > * getHoT()
param: no parameters
Definition: detector.hh:157
network::setDelayIndex
void setDelayIndex(double rate)
param: MDC log file
Definition: network.cc:2865
NIFO_MAX
const int NIFO_MAX
Definition: wat.hh:4
CWB::config::BATCH
int BATCH
Definition: config.hh:227
tmp
double * tmp
Definition: testWDM_5.C:31
network::a_00
wavearray< float > a_00
wdm multi-resolution analysis
Definition: network.hh:634
gIFACTOR
int gIFACTOR
Definition: cbc_plots_sim4.C:115
detector::fp
wavearray< double > fp
Definition: detector.hh:345
netcluster::p_Ind
std::vector< vector_int > p_Ind
Definition: netcluster.hh:385
network::cpf_
static void cpf_(float *&a, double **p)
Definition: network.hh:1086
CWB_Plugin
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *net, WSeries< double > *x, TString ifo, int type)
COHERENCE.
Definition: CWB_Plugin_supercluster_bench.C:25
network::pnt_
static void pnt_(float **, float **, short **, int, int)
Definition: network.hh:1073
k
int k
Definition: cwb_report_prod_2.C:150
PrintElapsedTime
void PrintElapsedTime(int job_elapsed_time, double cpu_time, TString info)
Definition: CWB_Plugin_supercluster_bench.C:141
netcluster::cpf
size_t cpf(const netcluster &, bool=false, int=0)
Definition: netcluster.cc:99
network::index
wavearray< int > index
Definition: network.hh:622
network::acor
double acor
Definition: network.hh:567
netpixel::time
size_t time
Definition: netpixel.hh:92
netcluster::loadTDampSSE
size_t loadTDampSSE(network &net, char c, size_t BATCH=10000, size_t LOUD=0)
Definition: netcluster.cc:3298
detector::lagShift
wavearray< double > lagShift
Definition: detector.hh:351
network::e2or
double e2or
Definition: network.hh:566
e
double e
Definition: CreateListEBBH.C:35
ifile
TFile * ifile
Definition: cwb_dump_config.C:23
rateANA
size_t rateANA
Definition: test_config1.C:21
netcluster::size
size_t size()
Definition: netcluster.hh:129
ww
WSeries< double > ww
Definition: Regression_H1.C:33
netcluster::getPixel
netpixel * getPixel(size_t n, size_t i)
Definition: netcluster.hh:395
CWB::config::subnet
double subnet
Definition: config.hh:229
network::events
size_t events()
Definition: network.hh:311
CWB::config::l_low
int l_low
Definition: config.hh:137
skymap::getPhi
double getPhi(size_t i)
Definition: skymap.hh:146
job_elapsed_time
int job_elapsed_time
Definition: cwb_net.C:718
netpixel::ellipticity
float ellipticity
Definition: netpixel.hh:100
_sse_add_ps
static void _sse_add_ps(__m128 *_a, __m128 *_b)
Definition: watsse.hh:230
netcluster::sCuts
std::vector< int > sCuts
Definition: netcluster.hh:374
CWB::config::LOUD
int LOUD
Definition: config.hh:228
netpixel::theta
float theta
Definition: netpixel.hh:98
fabs
double fabs(const Complex &x)
Definition: numpy.cc:37
netcluster::csize
size_t csize()
Definition: netcluster.hh:133
ifactor
int ifactor
Definition: cbc_plots_sim4.C:679
network::getwc
netcluster * getwc(size_t n)
param: delay index
Definition: network.hh:421
l
int l
Definition: cbc_plots.C:434
_sse_sum_ps
static __m128 _sse_sum_ps(__m128 **_p)
Definition: watsse.hh:497
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
netcluster::cData
std::vector< clusterdata > cData
Definition: netcluster.hh:373
_sse_point_ps
static void _sse_point_ps(__m128 **_p, float **p, short **m, int l, int n)
Definition: watsse.hh:484
_sse_mul_ps
static void _sse_mul_ps(__m128 *_a, float b)
Definition: watsse.hh:38
wavearray::data
DataType_t * data
Definition: wavearray.hh:301
CWB::config::Tgap
double Tgap
Definition: config.hh:116
network::rNRG
wavearray< float > rNRG
buffer for cluster sky 90 amplitudes
Definition: network.hh:636
wc
netcluster wc
Definition: DrawClusterWithSkyplot.C:39
K
int K
Definition: CWB_Plugin_Sim4_BRST_Config.C:87
nRES
const int nRES
Definition: revMonster.cc:7
id
Long_t id
Definition: cwb_condor_check.C:48
job_elapsed_hour
int job_elapsed_hour
Definition: cwb_net.C:719
network::skyMask
wavearray< short > skyMask
Definition: network.hh:623
_sse_maxE_ps
static float _sse_maxE_ps(__m128 *_a, __m128 *_A)
Definition: watsse.hh:536
netpixel::rate
float rate
Definition: netpixel.hh:95
_sse_minSNE_ps
static void _sse_minSNE_ps(__m128 *_pE, __m128 **_pe, __m128 *_es)
Definition: watsse.hh:522
netcluster::read
size_t read(const char *)
Definition: netcluster.cc:3096
_sse_like4_ps
static __m128 _sse_like4_ps(__m128 *_f, __m128 *_a, __m128 *_A)
Definition: watsse.hh:980
job_elapsed_sec
int job_elapsed_sec
Definition: cwb_net.C:721
wavearray::resize
virtual void resize(unsigned int)
Definition: wavearray.cc:445
CWB::config::TFgap
double TFgap
Definition: config.hh:120
netcluster::psize
size_t psize(int k=2)
Definition: netcluster.hh:145
network::netRHO
double netRHO
Definition: network.hh:579
CWB::config::inRate
size_t inRate
Definition: config.hh:114
_sse_dpf4_ps
static void _sse_dpf4_ps(__m128 *_Fp, __m128 *_Fx, __m128 *_fp, __m128 *_fx)
Definition: watsse.hh:613
CWB::config::levelR
int levelR
Definition: config.hh:134
CWB::config::scPlugin
bool scPlugin
Definition: config.hh:350
exit
exit(0)
netpixel::clean
void clean()
Definition: netpixel.hh:81
netcluster::clear
void clear()
Definition: netcluster.hh:409