Logo Coherent WaveBurst  
Reference Guide 		Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PCA_Benchmark.C
Go to the documentation of this file.
1 #include <vector>
2 
3 // --------------------------------------------------------
4 // PRINCIPAL COMPONENT ANALYSIS
5 // --------------------------------------------------------
6 
7 #define PCA_NRES 7 // number of resolutions
8 #define PCA_IRES 4 // initial resolution
9 
10 #define PCA_MAX 60 // number of max PC extractions
11 #define PCA_PREC 0.01 // precision = Eresidual/Esignal
12 
13 #define ACORE 1.7 // set selection threshold : Ethr = 2*ACORE*ACORE
14 
15 //#define APPLY_TSHIFT // enable/disable (uncomment/comment) find best wavelet time shift
16 
17 //#define OUTPUT_FILE "pca_benchmark.txt"
18  // enable/disable (uncomment/comment) write precision vs nPC to ascii file
19 
20 // --------------------------------------------------------
21 // SIGNAL
22 // --------------------------------------------------------
23 
24 #define WAVE_TYPE "INSPIRAL" // inspiral
25 //#define WAVE_TYPE "SGE" // sine gaussian
26 //#define WAVE_TYPE "WNB" // white noise burst
27 //#define WAVE_TYPE "GA" // gaussian
28 
29 #define WAVE_LENGTH 6 // Is effective only for WAVE_TYPE = SGE,WNB,GA (def=1sec)
30 
31 #define WAVE_POL "hp" // select hp,hx component
32 
33 #define WAVE_SNR 30 // signal SNR
34 
35 // --------------------------------------------------------
36 // WDM
37 // --------------------------------------------------------
38 
39 #define WDM_NU 6 // WDM mayer order
40  // NOTE : wp_xtalk coefficients are defined only for WDM_NU=2,4,6
41 #define WDM_PREC 10 // WDM precision
42 #define WDM_TDSIZE 12 // Time Delay filter size
43 //#define WDM_SCRATCH 14.0 // WDM TF segment scratch (sec)
44 #define WDM_SCRATCH 27.0 // WDM TF segment scratch (sec)
45 //#define WDM_SCRATCH 54.0 // WDM TF segment scratch (sec)
46 
47 //#define LOAD_XTALK_CATALOG // enable/disable (uncomment/comment) xcatalog
48 //#define XTALK_CATALOG "wdmXTalk/OverlapCatalog16-1024.bin" // WDM_NU = 4
49 #define XTALK_CATALOG "wdmXTalk/OverlapCatalog-ilLev4-hLev10-iNu6-P10.bin" // WDM_NU = 6
50 
51 // --------------------------------------------------------
52 // WAVELET PACKET
53 // --------------------------------------------------------
54 
55 #define NWP 4 // number of WP types used in the PCA
56 
57 #define WP_SINGLE // enable/disable (uncomment/comment) single wavelet
58 //#define WP_DIAG // enable/disable (uncomment/comment) diagonal WP
59 //#define WP_HORIZ // enable/disable (uncomment/comment) horizontal WP
60 //#define WP_VERT // enable/disable (uncomment/comment) vertical WP
61 //#define WP_BDIAG // enable/disable (uncomment/comment) back-diagonal WP
62 
63  // enable/disable (uncomment/comment) wavelet packet analysis
64  // NOTE : for more infos see InitXTALK function
65  // ----------------------------------------------------------
66  // Phase 00 Phase 90
67  // - - c*w - - C*W
68  // - b*v - - B*V -
69  // a*u - - A*U - -
70  //
71  // WAVELET PACKET BASE
72  // parity = odd -> z=U+v+W , Z=u-V+w
73  // parity = even -> z=U-v+W , Z=u+V+w
74  // (v,U) = (v,W) = -(V,u) = -(V,w) = WP_XTALK
75  // (W,u) = (U,w) = 0
76  // (z,z) = (Z,Z) = 3+4*WP_XTALK
77  // (z,Z) = 0
78  //
79  // DATA
80  // x=a*u+b*v+c*w
81  // A00 = (x,z)/(z,z) // A00 is the amplitude of x respect to z
82  // A90 = (x,Z)/(Z,Z) // A90 is the amplitude of x respect to Z
83  // parity = odd
84  // A00 = (x,z) = (x,U+v+W) = A+b+C
85  // A90 = (x,Z) = (x,u-V+w) = a-B+c
86  // parity = even
87  // A00 = (x,z) = (x,U-v+W) = A-b+C
88  // A90 = (x,Z) = (x,u+V+w) = a+B+c
89  // ----------------------------------------------------------
90 
91 // --------------------------------------------------------
92 // PLOTS
93 // --------------------------------------------------------
94 
95 //#define PLOT_MONSTER // enable/disable (uncomment/comment) show likelihood TFmap
96 #define PLOT_REC_VS_INJ // enable/disable (uncomment/comment) show time rec vs inj waveform
97 
98 //#define PLOT_TFMAP // enable/disable (uncomment/comment) show TF map @ PLOT_TFRES resolution
99 #define PLOT_TFRES 6 // TF frequency resolution showed
100 #define PLOT_MINFREQ 0 // min freq shoed in the TF plots
101 #define PLOT_MAXFREQ 256 // max freq shoed in the TF plots
102 #define PLOT_TFTYPE 2
103  // 0 : sqrt((E00+E90)/2)
104  // 1 : (E00+E90)/2
105  // 2 : sqrt((E00+E90)/2)
106  // 3 : amplitude:00
107  // 4 : energy:00
108  // 5 : |amplitude:00|
109  // 6 : amplitude:90
110  // 7 : energy:90
111  // 8 : |amplitude:90|
112 
113 // --------------------------------------------------------
114 // DECLARATIONS
115 // --------------------------------------------------------
116 
117 WDM<double> wdm[PCA_NRES]; // define a WDM transform containers
118 WSeries<double> wsig[PCA_NRES]; // waveform TF map containers
119 wavearray<double> sig; // signal
120 gwavearray<double>* grec; // used to plot rec
121 gwavearray<double>* gsig; // used to plot sig
122 watplot *WTS[3]; // plot objects
123 double wp_xtalk[4];
124 bool wp_mask[5];
125 
126 #ifdef LOAD_XTALK_CATALOG
127 monster wdmMRA; // wdm multi-resolution analysis
128 void CheckXTalkCatalog(int l_low, int l_high);
129 double GetPrincipalFactor(int nmax, int mmax, int rmax);
130 #endif
131 
132 int GetPrincipalComponent(double &amax, int &nmax, int &mmax, int &rmax);
133 
134 void InitXTALK();
135 
136 void PCA_Benchmark() {
137 
138  //
139  // Principal Component Analysis : Test Macro
140  // Extraction of Pincipal Components from multiresolution TF maps
141  // Author : Gabriele Vedovato
142 
143  if(TString("hp")!=WAVE_POL && TString("hx")!=WAVE_POL) {
144  cout << "Error - Bad WAVE_POL value : must be hp or hx" << endl;
145  exit(1);
146  }
147 #ifdef PLOT_TFMAP
148  if(PLOT_TFRES<0 || PLOT_TFRES>=PCA_NRES) {
149  cout << "Error - Bad PLOT_TFRES value : must be >=0 && < " << PCA_NRES << endl;
150  exit(1);
151  }
152 #endif
153 #ifdef OUTPUT_FILE
154  ofstream out;
155  out.open(OUTPUT_FILE,ios::out);
156 #endif
157 
158  for(int l=0;l<=NWP;l++) wp_mask[l] = false;
159 #ifdef WP_SINGLE
160  wp_mask[0] = true;
161 #endif
162 #ifdef WP_DIAG
163  wp_mask[1] = true;
164 #endif
165 #ifdef WP_HORIZ
166  wp_mask[2] = true;
167 #endif
168 #ifdef WP_VERT
169  wp_mask[3] = true;
170 #endif
171 #ifdef WP_BDIAG
172  wp_mask[4] = true;
173 #endif
174 
175  InitXTALK(); // init wp_xtalk coefficients
176 
177  CWB::mdc MDC;
178  MDC.SetInjLength(WAVE_LENGTH);
179 
180 #if (WAVE_TYPE == "INSPIRAL")
181  // ---------------------------------
182  // set inspiral parameters
183  // ---------------------------------
184  TString inspOptions="";
185  inspOptions = "--time-step 60.0 --time-interval 0 ";
186  inspOptions+= "--l-distr random ";
187  inspOptions+= "--gps-start-time 931072160 --gps-end-time 931072235 ";
188  inspOptions+= "--d-distr volume --m-distr totalMassRatio --i-distr uniform ";
189 
190 // inspOptions+= "--f-lower 80.000000 "; // 0.5 sec
191 // inspOptions+= "--f-lower 50.000000 "; // 1.5 sec
192 // inspOptions+= "--f-lower 40.000000 "; // 3 sec
193  inspOptions+= "--f-lower 30.000000 "; // 6.5 sec
194 // inspOptions+= "--f-lower 20.000000 "; // 19 sec
195 // inspOptions+= "--f-lower 10.000000 "; // 122 sec
196  inspOptions+= "--min-mass1 5.0 --max-mass1 5.0 ";
197  inspOptions+= "--min-mass2 5.0 --max-mass2 5.0 ";
198  inspOptions+= "--min-mtotal 10.0 --max-mtotal 10.0 ";
199 
200 /*
201  inspOptions+= "--f-lower 40.000000 ";
202  inspOptions+= "--min-mass1 10.0 --max-mass1 10.0 ";
203  inspOptions+= "--min-mass2 10.0 --max-mass2 10.0 ";
204  inspOptions+= "--min-mtotal 20.0 --max-mtotal 20.0 ";
205 */
206 /*
207  inspOptions+= "--f-lower 24.000000 ";
208  inspOptions+= "--min-mass1 42.1 --max-mass1 42.1 ";
209  inspOptions+= "--min-mass2 29.6 --max-mass2 29.6 ";
210  inspOptions+= "--min-mtotal 71.7 --max-mtotal 71.7 ";
211 */
212 /*
213  inspOptions+= "--f-lower 10.000000 ";
214  inspOptions+= "--min-mass1 75.000000 --max-mass1 75.000000 ";
215  inspOptions+= "--min-mass2 75.000000 --max-mass2 75.000000 ";
216  inspOptions+= "--min-mtotal 150.000000 --max-mtotal 150.000000 ";
217 */
218 
219  inspOptions+= "--m-distr componentMass ";
220  inspOptions+= "--min-distance 1000000.0 --max-distance 1500000.0 ";
221  inspOptions+= "--approximant EOBNRv2pseudoFourPN --disable-spin ";
222  inspOptions+= "--taper-injection start --seed 123456789 ";
223  inspOptions+= "--dir ./ ";
224  inspOptions+= "--output inspirals.xml "; // set output xml file
225 
226  MDC.SetInspiral("EOBNRv2",inspOptions);
227 
228  // Get the first waveform hp,hx components starting from gps = 931072160
229  sig = MDC.GetInspiral(WAVE_POL,931072160,931072235);
230 #endif
231 
232 #if (WAVE_TYPE == "SGE")
233  vector<mdcpar> par(2);
234  par[0].name="frequency"; par[0].value=100.;
235  par[1].name="Q"; par[1].value=8.9;
236  MDC.AddWaveform(MDC_SGE, par);
237  MDC.Print();
238  waveform wf = MDC.GetWaveform(0,0);
239  if(TString(WAVE_POL)=="hp") sig = wf.hp;
240  if(TString(WAVE_POL)=="hx") sig = wf.hx;
241 #endif
242 
243 #if (WAVE_TYPE == "WNB")
244  vector<mdcpar> par(6);
245  par[0].name="frequency"; par[0].value=250;
246  par[1].name="bandwidth"; par[1].value=100;
247  par[2].name="duration"; par[2].value=0.1;
248  par[3].name="pseed"; par[3].value=111;
249  par[4].name="xseed"; par[4].value=222;
250  par[5].name="mode"; par[5].value=1; // simmetric
251  MDC.AddWaveform(MDC_WNB, par);
252  MDC.Print();
253  waveform wf = MDC.GetWaveform(0,0);
254  if(TString(WAVE_POL)=="hp") sig = wf.hp;
255  if(TString(WAVE_POL)=="hx") sig = wf.hx;
256 #endif
257 
258 #if (WAVE_TYPE == "GA")
259  vector<mdcpar> par(1);
260  par[0].name="duration"; par[0].value=0.004;
261  MDC.AddWaveform(MDC_GA, par);
262  MDC.Print();
263  waveform wf = MDC.GetWaveform(0,0);
264  if(TString(WAVE_POL)=="hp") sig = wf.hp;
265  if(TString(WAVE_POL)=="hx") {
266  cout << "Error - WAVE=GA -> hx=0" << endl;
267  exit(1);
268  }
269 #endif
270 
271  cout << endl << "sig size : " << sig.size() << " sig rate : "
272  << sig.rate() << " sig start : " << (int)sig.start()
273  << " sig length : " << sig.size()/sig.rate() << " sec" << endl;
274  sig.start(0); // set start to 0 (needed by draw method)
275 
276  // add scratch array
277  int iscratch = WDM_SCRATCH*sig.rate();
278  wavearray<double> X(sig.size()+2*iscratch);
279  X.rate(sig.rate());
280  X=0.; for(int i=0;i<sig.size();i++) X[i+iscratch]=sig[i];
281  // apply time shift
282  //int jshift = 0.05*sig.rate();
283  //for(int i=jshift;i<sig.size();i++) X[i+iscratch]=sig[i-jshift];
284  sig=X;
285 
286  // resample data 16384 Hz -> 2048 Hz
287  sig.FFTW(1);
288  sig.resize(sig.size()/8);
289  sig.FFTW(-1);
290  sig.rate(sig.rate()/8.);
291 
292  // normalize sig energy to WAVE_SNR^2
293  double Esig=0.; for(int i=0;i<sig.size();i++) Esig+=sig[i]*sig[i];
294  for(int i=0;i<sig.size();i++) sig[i]/=sqrt(Esig)/sqrt(WAVE_SNR*WAVE_SNR);
295  Esig=WAVE_SNR*WAVE_SNR;
296 
297  // print pixel resolutions
298  cout << endl;
299  for(int level=PCA_NRES+PCA_IRES-1; level>=PCA_IRES; level--) {
300  int rateANA = int(sig.rate());
301  int layers = level>0 ? 1<<level : 0;
302  int rate = rateANA>>level;
303  cout << "level : " << level << "\t rate(hz) : " << rate << "\t layers : " << layers
304  << "\t df(hz) : " << rateANA/2./double(1<<level)
305  << "\t dt(ms) : " << 1000./rate << endl;
306  }
307  cout << endl;
308 
309 #ifdef APPLY_TSHIFT
310  cout << endl << "find best wavelet time shift ENABLED" << endl << endl;
311 #else
312  cout << endl << "find best wavelet time shift DISABLED" << endl << endl;
313 #endif
314  bool wpEnabled=false;
315  for(int l=0;l<=NWP;l++) if(wp_mask[l]) wpEnabled=true;
316  if(wpEnabled) {
317  cout << endl;
318  if(wp_mask[0]) cout << "wavelet packet single ENABLED" << endl;
319  if(wp_mask[1]) cout << "wavelet packet diagonal ENABLED" << endl;
320  if(wp_mask[2]) cout << "wavelet packet horizontal ENABLED" << endl;
321  if(wp_mask[3]) cout << "wavelet packet vertical ENABLED" << endl;
322  if(wp_mask[4]) cout << "wavelet packet back-diagonal ENABLED" << endl;
323  cout << endl;
324  } else {
325  cout << endl << "Error : all WP options are disabled !!!" << endl << endl;
326  exit(1);
327  }
328 #endif
329 
330  // compute WDM with PCA_NRES TF resolutions
331  cout << "Init WDM ..." << endl;
332  for(int i=PCA_NRES-1;i>=0;i--) {
333  int j=i+PCA_IRES;
334  wdm[i] = new WDM<double>(1<<j, 1<<j, WDM_NU, WDM_PREC); // define a WDM transform
335 
336  // check if filter lenght is less than cwb scratch length
337  double wdmFLen = double(wdm[i].m_H)/sig.rate(); // sec
338  if(wdmFLen > WDM_SCRATCH) {
339  cout << endl;
340  cout << "Error - filter length must be <= WDM_SCRATCH !!!" << " RES = " << j << endl;
341  cout << "filter length : " << wdmFLen << " sec" << endl;
342  cout << "cwb scratch : " << WDM_SCRATCH << " sec" << endl;
343  exit(1);
344  }
345  }
346 
347 #ifdef LOAD_XTALK_CATALOG
348 
349  char filter_dir[1024];
350  if(gSystem->Getenv("HOME_WAT_FILTERS")==NULL) {
351  cout << "Error : environment HOME_WAT_FILTERS is not defined!!!" << endl;exit(1);
352  } else {
353  strcpy(filter_dir,TString(gSystem->Getenv("HOME_WAT_FILTERS")).Data());
354  }
355 
356  char MRAcatalog[1024];
357  sprintf(MRAcatalog,"%s/%s",filter_dir,XTALK_CATALOG);
358  cout << "cwb2G::Init - Loading catalog of WDM cross-talk coefficients ... " << endl;
359  cout << MRAcatalog << endl;
360  CWB::Toolbox::checkFile(MRAcatalog);
361 
362  wdmMRA.read(MRAcatalog);
363 
364  // check if analysis layers are contained in the MRAcatalog
365  // level : is the decomposition level
366  // layes : are the number of layers along the frequency axis rateANA/(rateANA>>level)
367 
368  int l_low = PCA_IRES;
369  int l_high = PCA_IRES+PCA_NRES-1;
370 
371  int check_layers=0;
372  for(int level=l_high; level>=l_low; level--) {
373  int layers = level>0 ? 1<<level : 0;
374  for(int j=0;j<wdmMRA.nRes;j++) if(layers==wdmMRA.layers[j]) check_layers++;
375  }
376 
377  if(check_layers!=PCA_NRES) {
378  cout << "Error : analysis layers do not match the MRA catalog" << endl;
379  cout << endl << "analysis layers : " << endl;
380  for(int level=l_high; level>=l_low; level--) {
381  int layers = level>0 ? 1<<level : 0;
382  cout << "level : " << level << " layers : " << layers << endl;
383  }
384  cout << endl << "MRA catalog layers : " << endl;
385  for(int i=0;i<wdmMRA.nRes;i++)
386  cout << "layers : " << wdmMRA.layers[i] << endl;
387  exit(1);
388  }
389 
390  //CheckXTalkCatalog(l_low, l_high);
391 #endif
392 
393  // compute TF of sig signal
394  cout << "Apply WDM Transform to signal ..." << endl;
395  for(int i=0;i<PCA_NRES;i++) {
396  wsig[i].Forward(sig, wdm[i]); // apply the WDM to the time series
397  }
398 
399  // display TF sig
400 #ifdef PLOT_TFMAP
401  double dt0 = 1./wsig[PLOT_TFRES].rate();
402  double tbeg0 = WDM_SCRATCH-2*dt0;
403  double tend0 = dt0*wsig[PLOT_TFRES].size()/2-WDM_SCRATCH;
404  WTS[0] = new watplot(const_cast<char*>("wts0"));
405  WTS[0]->plot(&wsig[PLOT_TFRES], PLOT_TFTYPE, tbeg0, tend0, const_cast<char*>("COLZ"));
406  WTS[0]->hist2D->GetYaxis()->SetRangeUser(PLOT_MINFREQ, PLOT_MAXFREQ);
407  double hist2D_MAX = WTS[0]->hist2D->GetMaximum();
408 #endif
409 
410  // define time array of the PC reconstructed waveform
411  wavearray<double> rec(sig.size());
412  rec.rate(sig.rate());
413  rec=0;
414  double Erec=0.;
415 
416  // define netcluster for monster display
417  netcluster* pwc = new netcluster;
418  netpixel* pix = new netpixel;
419  pixdata pdata;
420  float crate = sig.rate();
421  float ctime;
422  float cfreq;
423  std::vector<int> list; // cluster list
424  std::vector<int> vtof(1); // time configuration array
425  std::vector<int> vtmp; // sky index array
426  std::vector<float> varea; // sky error regions array
427  std::vector<float> vpmap; // sky pixel map array
428  clusterdata cd; // dummy cluster data
429 
430  // start PC extraction
431  cout << "Start PC extraction ..." << endl;
432  double Ethr = 2*ACORE*ACORE; // selection pixel threshold
433  double precision=1;
434  int npix=0;
435  int nn,mm;
436  int nWP[NWP+1]; for(int l=0;l<=NWP;l++) nWP[l]=0;
437  int nPC=0;
438  for(int k=0;k<PCA_MAX;k++) {
439 
440  //if(fabs(Esig-Erec)/Esig<PCA_PREC) break;
441  if(precision<PCA_PREC) break;
442 
443  // find pixel with max energy -> principal component
444  double amax; // max pixel amplitude
445  int nmax; // time index @ amax
446  int mmax; // freq index @ amax
447  int rmax; // resolution @ amax
448  int isWP=GetPrincipalComponent(amax, nmax, mmax, rmax);
449  if(amax*amax<Ethr) break; // break loop if the PC energy is below the threshold
450  double factor = 1.;
451 #ifdef LOAD_XTALK_CATALOG
452  factor = GetPrincipalFactor(nmax, mmax, rmax);
453 #endif
454  nPC++; nWP[isWP]++;
455 
456  cout << endl << "--------------------------------------------------------" << endl;
457  cout << "Principal Component # " << k << endl;
458  cout << "--------------------------------------------------------" << endl;
459 
460  if(isWP) cout << "---> WAVELET PACKET : " << isWP << endl;
461  else cout << "---> wavelet pixel" << endl;
462  cout << " amax : " << amax << " nmax : " << nmax
463  << " mmax : " << mmax << " rmax : " << rmax << endl;
464 
465  double a00[3];
466  double a90[3];
467  a00[0] = wsig[rmax].getSample(nmax,mmax);
468  a90[0] = wsig[rmax].getSample(nmax,-mmax);
469 
470  if(isWP) { // wavelet packet : select 2 more pixels
471 
472  if(isWP==1) {nn=nmax+1; mm=mmax+1;}
473  if(isWP==2) {nn=nmax-1; mm=mmax; }
474  if(isWP==3) {nn=nmax; mm=mmax+1;}
475  if(isWP==4) {nn=nmax-1; mm=mmax+1;}
476 
477  a00[1] = wsig[rmax].getSample(nn,mm);
478  a90[1] = wsig[rmax].getSample(nn,-mm);
479 
480  if(isWP==1) {nn=nmax-1; mm=mmax-1;}
481  if(isWP==2) {nn=nmax+1; mm=mmax; }
482  if(isWP==3) {nn=nmax; mm=mmax-1;}
483  if(isWP==4) {nn=nmax+1; mm=mmax-1;}
484 
485  a00[2] = wsig[rmax].getSample(nn,mm);
486  a90[2] = wsig[rmax].getSample(nn,-mm);
487  }
488 
489  double ishift=0;
490 #ifdef APPLY_TSHIFT
491  // get time delayed amplitutes A00,A90
492  WDM<double>* pWDM = (WDM<double>*)wsig[rmax].pWavelet;
493  pWDM->setTDFilter(WDM_TDSIZE); // computes TD filters
494  int nsmp = (1<<rmax+PCA_IRES)/2;
495  double Emax=a00[0]*a00[0]+a90[0]*a90[0];
496  if(isWP) { // wavelet packet
497  Emax+=(a00[1]*a00[1]+a90[1]*a90[1]);
498  Emax+=(a00[2]*a00[2]+a90[2]*a90[2]);
499  }
500  double A00[3];
501  double A90[3];
502  A00[0]=a00[0]; A90[0]=a90[0];
503  if(isWP) {
504  A00[1]=a00[1]; A90[1]=a90[1];
505  A00[2]=a00[2]; A90[2]=a90[2];
506  }
507  double d00[3];
508  double d90[3];
509  double EMAX=0;
510  for(int i=-nsmp;i<=nsmp;i++) {
511  d00[0] = (double)pWDM->getPixelAmplitude(mmax, nmax, i, false); // delay a00 by i samples
512  d90[0] = (double)pWDM->getPixelAmplitude(mmax, nmax, i, true); // delay a90 by i samples
513  EMAX = d00[0]*d00[0]+d90[0]*d90[0];
514  if(EMAX>Emax) {Emax=EMAX;ishift=i;A00[0]=d00[0];A90[0]=d90[0];}
515 
516  if(isWP) { // wavelet packet : select 2 more pixels
517 
518  if(isWP==1) {nn=nmax+1; mm=mmax+1;}
519  if(isWP==2) {nn=nmax-1; mm=mmax; }
520  if(isWP==3) {nn=nmax; mm=mmax+1;}
521  if(isWP==4) {nn=nmax-1; mm=mmax+1;}
522 
523  d00[1] = (double)pWDM->getPixelAmplitude(mm, nn, i, false); // delay a00 by i samples
524  d90[1] = (double)pWDM->getPixelAmplitude(mm, nn, i, true); // delay a90 by i samples
525 
526  if(isWP==1) {nn=nmax-1; mm=mmax-1;}
527  if(isWP==2) {nn=nmax+1; mm=mmax; }
528  if(isWP==3) {nn=nmax; mm=mmax-1;}
529  if(isWP==4) {nn=nmax+1; mm=mmax-1;}
530 
531  d00[2] = (double)pWDM->getPixelAmplitude(mm, nn, i, false); // delay a00 by i samples
532  d90[2] = (double)pWDM->getPixelAmplitude(mm, nn, i, true); // delay a90 by i samples
533 
534  // parity
535  int sign;
536  if(isWP==1) sign = mmax%2 ? 1 :-1;
537  if(isWP==2) sign = mmax%2 ? -1 : 1;
538  if(isWP==3) sign = mmax%2 ? -1 : 1;
539  if(isWP==4) sign = mmax%2 ? -1 : 1;
540 
541  // get D00,D90 with respect to z,Z
542  // x = a90[1]*U +/- a00[0]*v + a90[2]*W
543  // X = a00[1]*u -/+ a90[0]*V + a00[2]*w
544  double D00 = 0;
545  double D90 = 0;
546  if(sign==1) {
547  D00 = a00[0]+a90[1]+a90[2]; // D00 = (x,z)
548  D90 = -a90[0]+a00[1]+a00[2]; // D00 = (x,z)
549  } else {
550  D00 = -a00[0]+a90[1]+a90[2]; // D00 = (x,z)
551  D90 = a90[0]+a00[1]+a00[2]; // D00 = (x,z)
552  }
553 
554  // D00 = (x,z)/(z,z)
555  // D90 = (X,Z)/(Z,Z)
556  // (z,z) = (Z,Z) = 3+4*WP_XTALK
557  D00/=(3+4*wp_xtalk[isWP-1]);
558  D90/=(3+4*wp_xtalk[isWP-1]);
559 
560  // z = z/sqrt(3+4*WP_XTALK)
561  // Z = Z/sqrt(3+4*WP_XTALK)
562  // (z,z) = (Z,Z) = 1
563  D00*=sqrt(3+4*wp_xtalk[isWP-1]);
564  D90*=sqrt(3+4*wp_xtalk[isWP-1]);
565 
566  EMAX = (D00*D00+D90*D90);
567 
568  if(EMAX>Emax) {
569  Emax=EMAX;ishift=i;
570  A00[0]=d00[0];A90[0]=d90[0];
571  A00[1]=d00[1];A90[1]=d90[1];
572  A00[2]=d00[2];A90[2]=d90[2];
573  }
574  }
575  }
576  cout << 0 << " a00[0] : " << a00[0] << " a90[0] : " << a90[0] << endl;
577  cout << ishift << " A00[0] : " << A00[0] << " A90[0] : " << A90[0] << endl;
578  cout << "tshift : " << 1000*ishift/wsig[rmax].rate() << " msec" << endl;
579  cout << "tmax : " << 1000*nsmp/wsig[rmax].rate() << " msec" << endl;
580  a00[0]=A00[0];a90[0]=A90[0];
581  if(isWP) {
582  a00[1]=A00[1];a90[1]=A90[1];
583  a00[2]=A00[2];a90[2]=A90[2];
584  }
585 #endif
586 
587  // get max pixels time representation
588 
589  int layers = wsig[rmax].maxLayer()+1; // numbers of frequency bins (first & last bins have df/2)
590  int slices = wsig[rmax].sizeZero(); // number of time bins
591  //cout << "layers " << layers << " slices " << slices << endl;
592 
593  int itime = nmax*layers + mmax;
594 
595  wavearray<double> w00[3]; //time series container
596  wavearray<double> w90[3]; //time series container
597  int j00[3];
598  int j90[3];
599 
600  j00[0] = wdm[rmax].getBaseWave(itime,w00[0],false)-ishift;
601  j90[0] = wdm[rmax].getBaseWave(itime,w90[0],true)-ishift;
602 
603  wavearray<double> x(sig.size()); // PC time series container
604  x.rate(sig.rate());
605  x=0;
606 
607  double D00 = 0;
608  double D90 = 0;
609  if(isWP) { // wavelet packet : select 2 more pixels
610 
611  if(isWP==1) {nn=nmax+1; mm=mmax+1;}
612  if(isWP==2) {nn=nmax-1; mm=mmax; }
613  if(isWP==3) {nn=nmax; mm=mmax+1;}
614  if(isWP==4) {nn=nmax-1; mm=mmax+1;}
615 
616  itime = nn*layers + mm;
617  j00[1] = wdm[rmax].getBaseWave(itime,w00[1],false)-ishift;
618  j90[1] = wdm[rmax].getBaseWave(itime,w90[1],true)-ishift;
619 
620  if(isWP==1) {nn=nmax-1; mm=mmax-1;}
621  if(isWP==2) {nn=nmax+1; mm=mmax; }
622  if(isWP==3) {nn=nmax; mm=mmax-1;}
623  if(isWP==4) {nn=nmax+1; mm=mmax-1;}
624 
625  itime = nn*layers + mm;
626  j00[2] = wdm[rmax].getBaseWave(itime,w00[2],false)-ishift;
627  j90[2] = wdm[rmax].getBaseWave(itime,w90[2],true)-ishift;
628 
629  // parity
630  // sign= 1 -> z=U+v+W , Z=u-V+w
631  // sign=-1 -> z=U-v+W , Z=u+V+w
632  // v=w00[0], u=w00[1], w=w00[2]
633  // V=w90[0], U=w90[1], W=w90[2]
634  // (z,z) = (Z,Z) = 3+4*WP_XTALK
635  // (z,Z) = 0
636  int sign;
637  if(isWP==1) sign = mmax%2 ? 1 :-1;
638  if(isWP==2) sign = mmax%2 ? -1 : 1;
639  if(isWP==3) sign = mmax%2 ? -1 : 1;
640  if(isWP==4) sign = mmax%2 ? -1 : 1;
641  cout << " parity : " << sign << endl;
642 
643  // get D00,D90 with respect to z,Z
644  // x = a90[1]*U +/- a00[0]*v + a90[2]*W
645  // X = a00[1]*u -/+ a90[0]*V + a00[2]*w
646  // (x,X)=0
647  if(sign==1) {
648  D00 = a00[0]+a90[1]+a90[2]; // D00 = (x,z)
649  D90 = -a90[0]+a00[1]+a00[2]; // D00 = (x,z)
650  } else {
651  D00 = -a00[0]+a90[1]+a90[2]; // D00 = (x,z)
652  D90 = a90[0]+a00[1]+a00[2]; // D00 = (x,z)
653  }
654  D00/=(3+4*wp_xtalk[isWP-1]); // D00 = (x,z)/(z,z)
655  D90/=(3+4*wp_xtalk[isWP-1]); // D90 = (x,Z)/(Z,Z)
656 
657  // D00*z
658  // sign= 1 -> z=U+v+W
659  // sign=-1 -> z=U-v+W
660  for(int i=0; i<w90[1].size(); i++){ // U=w90[1]
661  if(j90[1]+i<0 || j90[1]+i>=x.size()) continue;
662  x.data[j90[1]+i] += D00*w90[1][i];
663  }
664  for(int i=0; i<w00[0].size(); i++){ // v=w00[0]
665  if(j00[0]+i<0 || j00[0]+i>=x.size()) continue;
666  x.data[j00[0]+i] += sign*D00*w00[0][i];
667  }
668  for(int i=0; i<w90[2].size(); i++){ // v=w00[0]
669  if(j90[2]+i<0 || j90[2]+i>=x.size()) continue;
670  x.data[j90[2]+i] += D00*w90[2][i];
671  }
672 
673  // D90*Z
674  // sign= 1 -> Z=u-V+w
675  // sign=-1 -> Z=u+V+w
676  for(int i=0; i<w00[1].size(); i++){ // u=w00[1]
677  if(j00[1]+i<0 || j00[1]+i>=x.size()) continue;
678  x.data[j00[1]+i] += D90*w00[1][i];
679  }
680  for(int i=0; i<w90[0].size(); i++){ // V=w90[0]
681  if(j90[0]+i<0 || j90[0]+i>=x.size()) continue;
682  x.data[j90[0]+i] -= sign*D90*w90[0][i];
683  }
684  for(int i=0; i<w00[2].size(); i++){ // w=w00[2]
685  if(j00[2]+i<0 || j00[2]+i>=x.size()) continue;
686  x.data[j00[2]+i] += D90*w00[2][i];
687  }
688 
689  } else { // single pixel
690 
691  for(int i=0; i<w00[0].size(); i++){
692  if(j00[0]+i<0 || j00[0]+i>=x.size()) continue;
693  x.data[j00[0]+i] += factor*a00[0]*w00[0][i];
694  }
695  for(int i=0; i<w90[0].size(); i++){
696  if(j90[0]+i<0 || j90[0]+i>=x.size()) continue;
697  x.data[j90[0]+i] += factor*a90[0]*w90[0][i];
698  }
699 
700  D00 = a00[0];
701  D90 = a90[0];
702  }
703 
704  // fill pixel infos for monster display
705  npix++;
706  pix->core = 1;
707  pix->clusterID = 0; // setup new ID
708  pix->rate = x.rate()/(1<<(rmax+PCA_IRES));
709  pix->layers = layers;
710  pix->time = itime;
711  pix->frequency = mmax;
712  pix->data.clear();
713  pdata.asnr = D00;
714  pdata.a_90 = D90;
715  pix->data.push_back(pdata);
716  pwc->pList.push_back(*pix); // update pList and counter
717 
718  // PC waveform
719  for(int i=0; i<x.size(); i++) rec.data[i] += x[i];
720  Erec=0.; for(int i=0; i<rec.size(); i++) Erec+=rec[i]*rec[i];
721 
722  // subtract max pixels from TF maps
723  WSeries<double> wx; // TF map container
724  for(int i=0;i<PCA_NRES;i++) {
725  wx.Forward(x, wdm[i]); // apply the WDM to the max pixel time series
726  wsig[i]-=wx; // subtract max pixels from TF maps
727  }
728 
729  // compute precision = Eresidual/Esignal
730  wavearray<double> y = rec;
731  y-=sig; y*=y;
732  wavearray<double> y2 = sig;
733  y2*=y2;
734  precision = y.mean()/y2.mean();
735  cout << " Erec " << Erec << " Esig " << Esig
736  << " fabs(Esig-Erec)/Esig : " << fabs(Esig-Erec)/Esig
737  << " precision " << precision << endl;
738 #ifdef OUTPUT_FILE
739  out << Erec << " " << Esig << " " << fabs(Esig-Erec)/Esig << " " << precision << endl;
740 #endif
741  }
742 #ifdef OUTPUT_FILE
743  out.close();
744 #endif
745 
746  cout << endl;
747  cout << "-----------------------------------------------------" << endl;
748  cout << "number of SW/nPC used in the PCA : " << nWP[0] << "/" << nPC << endl;
749  cout << "number of diagonal_WP/nPC used in the PCA : " << nWP[1] << "/" << nPC << endl;
750  cout << "number of horizontal_WP/nPC used in the PCA : " << nWP[2] << "/" << nPC << endl;
751  cout << "number of vertical_WP/nPC used in the PCA : " << nWP[3] << "/" << nPC << endl;
752  cout << "number of back-diagonal_WP/nPC used in the PCA : " << nWP[4] << "/" << nPC << endl;
753 #ifdef OUTPUT_FILE
754  cout << "output file : " << OUTPUT_FILE << endl;
755 #endif
756  cout << "-----------------------------------------------------" << endl;
757  cout << endl;
758 
759  // display tf (residuals)
760 #ifdef PLOT_TFMAP
761  double dt1 = 1./wsig[PLOT_TFRES].rate();
762  double tbeg1 = WDM_SCRATCH-2*dt1;
763  double tend1 = dt1*wsig[PLOT_TFRES].size()/2-WDM_SCRATCH;
764  WTS[1] = new watplot(const_cast<char*>("wts1"));
765  WTS[1]->plot(&wsig[PLOT_TFRES], PLOT_TFTYPE, tbeg1, tend1, const_cast<char*>("COLZ"));
766  WTS[1]->hist2D->GetYaxis()->SetRangeUser(PLOT_MINFREQ, PLOT_MAXFREQ);
767  WTS[1]->hist2D->GetZaxis()->SetRangeUser(0, hist2D_MAX);
768 #endif
769 
770 #ifdef PLOT_MONSTER
771  // monster display
772  for(int i=0;i<npix;i++) list.push_back(i);
773  // fill pixel infos for monster display
774  pwc->rate = sig.rate();
775  pwc->cList.push_back(list);
776  pwc->cRate.push_back(crate);
777  pwc->cTime.push_back(ctime);
778  pwc->cFreq.push_back(cfreq);
779  pwc->sArea.push_back(varea); // recreate sky error regions array
780  pwc->p_Map.push_back(vpmap); // recreate sky pixel map array
781  pwc->nTofF.push_back(vtof); // recreate time configuration array
782  pwc->p_Ind.push_back(vtmp); // recreate sky index array
783  pwc->cData.push_back(cd); // recreate dummy cluster data
784 
785  WTS[2] = new watplot(const_cast<char*>("wts2"));
786  WTS[2]->plot(pwc, 1, 1, 'L', 0, const_cast<char*>("COLZ"));
787 #endif
788 
789 #ifdef PLOT_REC_VS_INJ
790  // display rec vs sig waveform in time domain
791  grec = new gwavearray<double>;
792  *grec = rec;
793  gsig = new gwavearray<double>;
794  *gsig = sig;
795  gsig->SetTimeRange(WDM_SCRATCH, gsig->size()/gsig->rate()-WDM_SCRATCH);
796  gsig->Draw(GWAT_TIME,"ALPCUSTOM");
797  gsig->Draw(grec,GWAT_TIME,"SAME",kRed);
798 #endif
799 }
800 
801 int GetPrincipalComponent(double &amax, int &nmax, int &mmax, int &rmax) {
802 
803  rmax=0;
804  nmax=0;
805  mmax=0;
806  amax=0;
807 
808  int Rmax=0;
809  double Nmax=0;
810  double Mmax=0;
811  double Amax=0;
812 
813  double Ethr = 2*ACORE*ACORE; // selection pixel threshold
814  int isWP=0;
815  int nn,mm;
816  double a00[3];
817  double a90[3];
818  double A00 = 0;
819  double A90 = 0;
820 
821  for(int i=0;i<PCA_NRES;i++) {
822 
823  int layers = wsig[i].maxLayer()+1; // numbers of frequency bins (first & last bins have df/2)
824  int slices = wsig[i].sizeZero(); // number of time bins
825  //cout << "layers " << layers << " slices " << slices << endl;
826 
827  int iscratch = WDM_SCRATCH*(wsig[i].rate()/(1<<(i+PCA_IRES)));
828  for(int n=iscratch;n<slices-iscratch;n++) { // exclude scratch data
829  for(int m=2;m<layers-2;m++) {
830  a00[0] = wsig[i].getSample(n,m);
831  a90[0] = wsig[i].getSample(n,-m);
832 
833  double E = (a00[0]*a00[0]+a90[0]*a90[0]);
834 
835  double A = sqrt(E);
836  if(wp_mask[0]) {
837  if(A>amax) {amax=A; nmax=n; mmax=m; rmax=i; isWP=0;}
838  if(A>Amax) {Amax=A; Nmax=n; Mmax=m; Rmax=i; A00=a00[0]; A90=a90[0];}
839  }
840 
841  for(int l=0;l<NWP;l++) {
842 
843  if(!wp_mask[l+1]) continue;
844 
845  if(l==0) {nn=n+1; mm=m+1;} // diagonal WP
846  if(l==1) {nn=n-1; mm=m; } // horizontal WP
847  if(l==2) {nn=n; mm=m+1;} // vertical WP
848  if(l==3) {nn=n-1; mm=m+1;} // back-diagonal WP
849 
850  a00[1] = wsig[i].getSample(nn,mm);
851  a90[1] = wsig[i].getSample(nn,-mm);
852 
853  if(l==0) {nn=n-1; mm=m-1;} // diagonal WP
854  if(l==1) {nn=n+1; mm=m; } // horizontal WP
855  if(l==2) {nn=n; mm=m-1;} // vertical WP
856  if(l==3) {nn=n+1; mm=m-1;} // back-diagonal WP
857 
858  a00[2] = wsig[i].getSample(nn,mm);
859  a90[2] = wsig[i].getSample(nn,-mm);
860 
861  // parity
862  int sign;
863  if(l==0) sign = m%2 ? 1 :-1;
864  if(l==1) sign = m%2 ? -1 : 1;
865  if(l==2) sign = m%2 ? -1 : 1;
866  if(l==3) sign = m%2 ? -1 : 1;
867 
868  // sign= 1 -> z=U+v+W , Z=u-V+w
869  // sign=-1 -> z=U-v+W , Z=u+V+w
870  // v=w00[0], u=w00[1], w=w00[2]
871  // V=w90[0], U=w90[1], W=w90[2]
872  // (v,U) = (v,W) = -(V,u) = -(V,w) = WP_XTALK
873  // (z,Z) = 0
874 
875  // get D00,D90 with respect to z,Z
876  // x = a90[1]*U +/- a00[0]*v + a90[2]*W
877  // X = a00[1]*u -/+ a90[0]*V + a00[2]*w
878  double D00 = 0;
879  double D90 = 0;
880  if(sign==1) {
881  D00 = a00[0]+a90[1]+a90[2]; // D00 = (x,z)
882  D90 = -a90[0]+a00[1]+a00[2]; // D00 = (x,z)
883  } else {
884  D00 = -a00[0]+a90[1]+a90[2]; // D00 = (x,z)
885  D90 = a90[0]+a00[1]+a00[2]; // D00 = (x,z)
886  }
887  // (z,z) = (Z,Z) = 3+4*WP_XTALK
888  D00/=(3+4*wp_xtalk[l]);
889  D90/=(3+4*wp_xtalk[l]);
890 
891  // get D00,D90 with respect to the normalized z,Z
892  // z = z/sqrt(3+4*WP_XTALK)
893  // Z = Z/sqrt(3+4*WP_XTALK)
894  // (z,z) = (Z,Z) = 1
895  D00*=sqrt(3+4*wp_xtalk[l]);
896  D90*=sqrt(3+4*wp_xtalk[l]);
897 
898  double E = (D00*D00+D90*D90);
899 
900  A = sqrt(E);
901 
902  if(A>amax) {amax=A; nmax=n; mmax=m; rmax=i; isWP=l+1;}
903  }
904  // set isWP=0 if the difference between 1 pixels and WP is below the noise threshold
905  // 2*Ethr is the noise of the 2 WP extra pixels
906 // if(fabs(amax*amax-Amax*Amax)<2*Ethr) {amax=Amax; nmax=Nmax; mmax=Mmax; rmax=Rmax; isWP=0;}
907  }
908  }
909  }
910  return isWP;
911 }
912 
913 #ifdef LOAD_XTALK_CATALOG
914 void CheckXTalkCatalog(int l_low, int l_high) {
915 
916  for(int level=l_high; level>=l_low; level--) {
917  int layers = level>0 ? (1<<level)+1 : 0;
918  cout << "layers : " << layers << endl;
919 
920  // v,V
921  int _n = 100;
922  int _m = 9; // parity = m%2 ? even : odd; m=9 -> parity=even
923  // w,W
924  int _nn=_n+1;
925  int _mm=_m+1;
926 
927 /*
928  for(int ii=-1;ii<=1;ii++) {
929  for(int jj=-1;jj<=1;jj++) {
930  // CC[0] CC[2]
931  // CC[1] CC[3]
932  _nn=_n+ii;
933  _mm=_m+jj;
934  struct xtalk xt = wdmMRA.getXTalk(layers, _n*layers+_m, layers, _nn*layers+_mm);
935  cout << ii << " " << jj << " (v,w) : " << xt.CC[0] << " (v,W) : " << xt.CC[1]
936  << " (V,w) : " << xt.CC[2] << " (V,W) : " << xt.CC[3] << endl;
937  }
938  }
939 */
940 
941  struct xtalk xt = wdmMRA.getXTalk(layers, _n*layers+_m, layers, _nn*layers+_mm);
942  cout << "xtalk catalog : " << " (v,w) : " << xt.CC[0] << " (v,W) : " << xt.CC[1]
943  << " (V,w) : " << xt.CC[2] << " (V,W) : " << xt.CC[3] << endl;
944 
945  int j00[2];
946  int j90[2];
947  wavearray<double> w00[2]; //time series container
948  wavearray<double> w90[2]; //time series container
949 
950  j00[0] = wdm[level-l_low].getBaseWave( _n*layers+_m, w00[0],false); // v
951  j90[0] = wdm[level-l_low].getBaseWave( _n*layers+_m, w90[0],true); // V
952  j00[1] = wdm[level-l_low].getBaseWave(_nn*layers+_mm,w00[1],false); // w
953  j90[1] = wdm[level-l_low].getBaseWave(_nn*layers+_mm,w90[1],true); // W
954 
955  wavearray<double> W00[2]; //time series container
956  wavearray<double> W90[2]; //time series container
957  for(int k=0;k<2;k++) {
958  W00[k]=sig; W00[k]=0;
959  for(int i=0;i<w00[k].size();i++) {
960  if(i+j00[k]<0 || i+j00[k]>W00[k].size()) continue;
961  W00[k][i+j00[k]] = w00[k][i];
962  }
963  W90[k]=sig; W90[k]=0;
964  for(int i=0;i<w90[k].size();i++) {
965  if(i+j90[k]<0 || i+j90[k]>W90[k].size()) continue;
966  W90[k][i+j90[k]] = w90[k][i];
967  }
968  }
969 
970  double vw=0.;
971  double vW=0.;
972  double Vw=0.;
973  double VW=0.;
974  for(int i=0;i<W90[0].size();i++) {
975  vw+=W00[0][i]*W00[1][i]; // (v,w)
976  vW+=W00[0][i]*W90[1][i]; // (v,W)
977  Vw+=W90[0][i]*W00[1][i]; // (V,w)
978  VW+=W90[0][i]*W90[1][i]; // (V,W)
979  }
980  // NOTE : for more infos see InitXTALK function
981  cout << "xtalk base waves : " << " (v,w) : " << vw << " (v,W) : " << vW
982  << " (V,w) : " << Vw << " (V,W) : " << VW << endl << endl;
983  cout << "wp_xtalk : " << " (v,w) : " << 0 << " (v,W) : " << wp_xtalk[0]
984  << " (V,w) : " << -wp_xtalk[0] << " (V,W) : " << 0 << endl << endl;
985  }
986 }
987 #endif
988 
989 double GetPrincipalFactor(int nmax, int mmax, int rmax) {
990 
991  double Ethr = 2*ACORE*ACORE; // selection pixel threshold
992  int nn,mm;
993  vector<double> vb;
994  vector<double> vB;
995  vector<int> vN;
996  vector<int> vM;
997  vector<int> vL;
998 
999  struct xtalk xt;
1000 
1001  for(int i=0;i<PCA_NRES;i++) {
1002 
1003  int layers = wsig[i].maxLayer()+1; // numbers of frequency bins (first & last bins have df/2)
1004  int slices = wsig[i].sizeZero(); // number of time bins
1005  //cout << "layers " << layers << " slices " << slices << endl;
1006 
1007  int iscratch = WDM_SCRATCH*(wsig[i].rate()/(1<<(i+PCA_IRES)));
1008  for(int n=iscratch;n<slices-iscratch;n++) { // exclude scratch data
1009  for(int m=2;m<layers-2;m++) {
1010  double a00 = wsig[i].getSample(n,m);
1011  double a90 = wsig[i].getSample(n,-m);
1012 
1013  double E = (a00*a00+a90*a90);
1014 
1015  if(E<Ethr) continue;
1016 
1017  //cout << "r:n:m " << i << ":" << n << ":" << m << " E : " << E << endl;
1018  //cout << "lmax:nmax:mmax " << lmax << ":" << nmax << ":" << mmax << endl;
1019 
1020  vb.push_back(a00);
1021  vB.push_back(a90);
1022  vN.push_back(n);
1023  vM.push_back(m);
1024  vL.push_back(layers);
1025  }
1026  }
1027  }
1028 
1029  // PC is X = a*Y00 + A*Y90
1030  // x - q*X
1031  // sum_i((x,Yi)-q*(X,Yi))^2 = sum_i(bi-q*(X,Yi))^2 + sum_i(Bi-q*(X,Yi))^2
1032  // q = (sum_i(bi*(X,Yi)) + sum_i(Bi*(X,Yi))) / (sum_i((X,Yi)^2) + sum_i((X,Yi)^2))
1033 
1034  int lmax = wsig[rmax].maxLayer()+1;
1035 
1036  double a = wsig[rmax].getSample(nmax,mmax);
1037  double A = wsig[rmax].getSample(nmax,-mmax);
1038 
1039  double num=0;
1040  double den=0;
1041  for(int k=0;k<vb.size();k++) {
1042  double b = vb[k];
1043  double B = vB[k];
1044  int n = vN[k];
1045  int m = vM[k];
1046  int l = vL[k];
1047  xt = wdmMRA.getXTalk(lmax, nmax*lmax+mmax, l, n*l+m);
1048  if(xt.CC[0]>2) continue;
1049  //cout << "r:n:m " << l << ":" << n << ":" << m << " b : " << b << " B : " << B << endl;
1050  //cout << "lmax:nmax:mmax " << lmax << ":" << nmax << ":" << mmax << endl;
1051  //cout << k << " (u,w) : " << xt.CC[0] << " (u,W) : " << xt.CC[1]
1052  // << " (U,w) : " << xt.CC[2] << " (U,W) : " << xt.CC[3] << endl;
1053 
1054  num+=b*(xt.CC[0]*a+xt.CC[2]*A)+B*(xt.CC[1]*a+xt.CC[3]*A);
1055  den+=pow(xt.CC[0]*a+xt.CC[2]*A,2)+pow(xt.CC[1]*a+xt.CC[3]*A,2);
1056  }
1057 
1058  double q = num/den;
1059 
1060  cout << "q = " << q << endl;
1061 
1062  return q;
1063 }
1064 
1065 void InitXTALK() {
1066 
1067 // -------------------------------------------------------------------------
1068 // XTALK DEFINITIONS
1069 // NOTE : All xtalk values are computed with 'v,V' parity = even (see below)
1070 // -------------------------------------------------------------------------
1071 //
1072 // WAVELET PACKET BASE
1073 // parity = m%2 ? even : odd; (n,m) are the 'v,V' time,frequency index
1074 // parity = odd -> z=U+v+W , Z=u-V+w
1075 // parity = even -> z=U-v+W , Z=u+V+w
1076 // (z,z) = (Z,Z) = 3+4*WP_XTALK
1077 // (z,Z) = 0
1078 //
1079 // diagonal WP
1080 // - - w - - W
1081 // - v - - V -
1082 // u - - U - -
1083 //
1084 // horizontal WP
1085 // - - - - - -
1086 // u v w U V W
1087 // - - - - - -
1088 //
1089 // vertical WP
1090 // - w - - W -
1091 // - v - - V -
1092 // - u - - U -
1093 //
1094 // back-diagonal WP
1095 // w - - W - -
1096 // - v - - V -
1097 // - - u - - U
1098 //
1099 // wp_xtalk = (v,W) = (W,v) = -(V,w) = -(w,V)
1100 // wp_xtalk = (v,U) = (U,v) = -(V,u) = -(u,V)
1101 //
1102 
1103 #if WDM_NU == 2
1104 
1105  wp_xtalk[0] = 0.207345; // diagonal WP (v,W)=(n:m, n+1:m+1) (v,U)=(n:m, n-1:m-1)
1106  wp_xtalk[1] = 0.561945; // horizontal WP (v,U)=(n:m, n-1:m ) (v,W)=(n:m, n+1:m )
1107  wp_xtalk[2] = 0.243038; // vertical WP (v,W)=(n:m, n :m+1) (v,U)=(n:m, n :m-1)
1108  wp_xtalk[3] = -0.207345; // back-diagonal WP (v,W)=(n:m, n-1:m+1) (v,U)=(n:m, n+1:m-1)
1109 
1110 #elif WDM_NU == 4
1111 
1112  wp_xtalk[0] = 0.162725; // diagonal WP (v,W)=(n:m, n+1:m+1) (v,U)=(n:m, n-1:m-1)
1113  wp_xtalk[1] = 0.597494; // horizontal WP (v,U)=(n:m, n-1:m ) (v,W)=(n:m, n+1:m )
1114  wp_xtalk[2] = 0.178856; // vertical WP (v,W)=(n:m, n :m+1) (v,U)=(n:m, n :m-1)
1115  wp_xtalk[3] = -0.162725; // back-diagonal WP (v,W)=(n:m, n-1:m+1) (v,U)=(n:m, n+1:m-1)
1116 
1117 #elif WDM_NU == 6
1118 
1119  wp_xtalk[0] = 0.138369; // diagonal WP (v,W)=(n:m, n+1:m+1) (v,U)=(n:m, n-1:m-1)
1120  wp_xtalk[1] = 0.610111; // horizontal WP (v,U)=(n:m, n-1:m ) (v,W)=(n:m, n+1:m )
1121  wp_xtalk[2] = 0.148011; // vertical WP (v,W)=(n:m, n :m+1) (v,U)=(n:m, n :m-1)
1122  wp_xtalk[3] = -0.138369; // back-diagonal WP (v,W)=(n:m, n-1:m+1) (v,U)=(n:m, n+1:m-1)
1123 
1124 #else
1125 
1126  bool wpCheck=false;
1127  for(int l=1;l<=NWP;l++) if(wp_mask[l]) wpCheck=true;
1128  if(wpCheck) {
1129  if(WDM_NU!=2) {
1130  cout << endl << "Error : wp_xtalk coefficients are defined only for WDM_NU=2,4,6 !!!" << endl << endl;
1131  exit(1);
1132  }
1133  }
1134 
1135 #endif
1136 
1137 }
1138 
wavearray::size
virtual size_t size() const
Definition: wavearray.hh:127
WDM::getBaseWave
int getBaseWave(int m, int n, SymmArray< double > &w)
Definition: WDM.cc:305
gsig
gwavearray< double > * gsig
Definition: PCA_Benchmark.C:121
netcluster::cRate
std::vector< vector_int > cRate
Definition: netcluster.hh:380
slices
int slices
Definition: DrawClusterWithSkyplot.C:30
filter_dir
char filter_dir[1024]
Definition: cwb1G_parameters.C:163
waveform
Definition: mdc.hh:189
j00
int j00
Definition: BaseFunctionWDM.C:52
vN
std::vector< WSeries< double > > vN[NIFO_MAX]
Definition: CWB_Plugin_xWRC.C:84
PCA_NRES
#define PCA_NRES
Definition: PCA_Benchmark.C:7
netpixel::clusterID
size_t clusterID
Definition: netpixel.hh:91
wavearray::rate
virtual void rate(double r)
Definition: wavearray.hh:123
PCA_IRES
#define PCA_IRES
Definition: PCA_Benchmark.C:8
a
wavearray< double > a(hp.size())
netpixel::frequency
size_t frequency
Definition: netpixel.hh:93
MDC_GA
Definition: mdc.hh:128
B
#define B
wdm
WDM< double > wdm[PCA_NRES]
Definition: PCA_Benchmark.C:117
n
int n
Definition: cwb_net.C:10
const_cast< char * >
cout<< "skymap size : "<< L<< endl;for(int l=0;l< L;l++) sm.set(l, l);sm > const_cast< char * >("skymap.dat")
XTALK_CATALOG
#define XTALK_CATALOG
Definition: PCA_Benchmark.C:49
TString
TString("c")
Definition: cwb_report_skymap.C:111
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::sArea
std::vector< vector_float > sArea
Definition: netcluster.hh:383
pix
netpixel pix(nifo)
pwc
netcluster * pwc
Definition: cwb_job_obj.C:20
CWB::mdc::SetInjLength
void SetInjLength(double inj_length=MDC_INJ_LENGTH)
Definition: mdc.hh:272
layers
int layers
Definition: DrawClusterWithSkyplot.C:29
waveform::hp
wavearray< double > hp
Definition: mdc.hh:196
WAVE_POL
#define WAVE_POL
Definition: PCA_Benchmark.C:31
MDC
CWB::mdc * MDC
Definition: CWB_Plugin_HEN_SIM_Config.C:78
WDM_TDSIZE
#define WDM_TDSIZE
Definition: PCA_Benchmark.C:42
wf
waveform wf
Definition: CWB_Plugin_HEN_SIM_Config.C:82
PCA_Benchmark
void PCA_Benchmark()
Definition: PCA_Benchmark.C:136
size
Long_t size
Definition: cwb_condor_check.C:48
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
wavearray::start
virtual void start(double s)
Definition: wavearray.hh:119
j
int j
Definition: cwb_net.C:10
cd
canvas cd(1)
WTS
watplot * WTS[3]
Definition: PCA_Benchmark.C:122
i
i drho i
Definition: cwb_epparameters.C:85
l_low
int l_low
Definition: test_config1.C:40
netcluster::rate
double rate
Definition: netcluster.hh:360
NWP
#define NWP
Definition: PCA_Benchmark.C:55
WDM::setTDFilter
void setTDFilter(int nCoeffs, int L=1)
Definition: WDM.cc:621
ACORE
#define ACORE
Definition: PCA_Benchmark.C:13
list
TList * list
Definition: cwb_dump_events.C:31
PCA_PREC
#define PCA_PREC
Definition: PCA_Benchmark.C:11
PLOT_MINFREQ
#define PLOT_MINFREQ
Definition: PCA_Benchmark.C:100
CWB::Toolbox::checkFile
static bool checkFile(TString fName, bool question=false, TString message="")
Definition: Toolbox.cc:3956
watplot::plot
void plot(wavearray< double > &, char *=NULL, int=1, double=0., double=0., bool=false, float=0., float=0., bool=false, float=0., bool=false)
Definition: watplot.cc:132
wavearray::mean
virtual double mean() const
Definition: wavearray.cc:1053
pixdata::asnr
double asnr
Definition: netpixel.hh:20
InitXTALK
void InitXTALK()
Definition: PCA_Benchmark.C:1065
CWB::mdc::AddWaveform
mdcid AddWaveform(MDC_TYPE mdc_type, vector< mdcpar > par, TString uname="")
Definition: mdc.cc:445
PLOT_TFRES
#define PLOT_TFRES
Definition: PCA_Benchmark.C:99
CWB::mdc::Print
void Print(int level=0)
Definition: mdc.cc:2707
netpixel::core
bool core
Definition: netpixel.hh:102
WDM::getPixelAmplitude
double getPixelAmplitude(int, int, int, bool=false)
Definition: WDM.cc:730
netcluster::p_Map
std::vector< vector_float > p_Map
Definition: netcluster.hh:384
watplot::hist2D
TH2F * hist2D
Definition: watplot.hh:175
WDM_SCRATCH
#define WDM_SCRATCH
Definition: PCA_Benchmark.C:44
out
ofstream out
Definition: cwb_merge.C:196
factor
double factor
Definition: ReadCStrainFromJobFile.C:18
PCA_MAX
#define PCA_MAX
Definition: PCA_Benchmark.C:10
CWB::mdc
Definition: mdc.hh:216
GetPrincipalComponent
int GetPrincipalComponent(double &amax, int &nmax, int &mmax, int &rmax)
Definition: PCA_Benchmark.C:801
xtalk
Definition: monster.hh:12
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))
netcluster::pList
std::vector< netpixel > pList
Definition: netcluster.hh:372
wp_xtalk
double wp_xtalk[4]
Definition: PCA_Benchmark.C:123
PLOT_MAXFREQ
#define PLOT_MAXFREQ
Definition: PCA_Benchmark.C:101
netcluster::p_Ind
std::vector< vector_int > p_Ind
Definition: netcluster.hh:385
GetPrincipalFactor
double GetPrincipalFactor(int nmax, int mmax, int rmax)
Definition: PCA_Benchmark.C:989
netcluster::cFreq
std::vector< float > cFreq
Definition: netcluster.hh:382
monster::read
void read(char *filename)
param: file name
Definition: monster.cc:244
precision
double precision
Definition: cwb1G_parameters.C:145
inspOptions
TString inspOptions
Definition: CWB_Plugin_InjectCBC_Config.C:12
k
int k
Definition: cwb_report_prod_2.C:150
MDC_SGE
Definition: mdc.hh:122
A
static double A
Definition: geodesics.cc:8
WAVE_SNR
#define WAVE_SNR
Definition: PCA_Benchmark.C:33
gwavearray::SetTimeRange
void SetTimeRange(double tMin, double tMax)
Definition: gwavearray.hh:57
netpixel::time
size_t time
Definition: netpixel.hh:92
par
vector< mdcpar > par
Definition: CWB_Plugin_HEN_SIM_Config.C:83
num
int num
Definition: check_slag_bkg.C:71
wsig
WSeries< double > wsig[PCA_NRES]
Definition: PCA_Benchmark.C:118
rateANA
size_t rateANA
Definition: test_config1.C:21
npix
int npix
Definition: ReadHistFromJobFile.C:66
WDM_NU
#define WDM_NU
Definition: PCA_Benchmark.C:39
wavearray::FFTW
virtual void FFTW(int=1)
Definition: wavearray.cc:878
waveform::hx
wavearray< double > hx
Definition: mdc.hh:197
PLOT_TFTYPE
#define PLOT_TFTYPE
Definition: PCA_Benchmark.C:102
netcluster::cTime
std::vector< float > cTime
Definition: netcluster.hh:381
pWDM
WDM< double > * pWDM
Definition: testWDM_5.C:28
fabs
double fabs(const Complex &x)
Definition: numpy.cc:37
wp_mask
bool wp_mask[5]
Definition: PCA_Benchmark.C:124
monster::layers
int * layers
Definition: monster.hh:103
CWB::mdc::GetWaveform
waveform GetWaveform(int ID, int id=0)
Definition: mdc.cc:2512
WSeries::Forward
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Definition: wseries.cc:228
strcpy
strcpy(RunLabel, RUN_LABEL)
l
int l
Definition: cbc_plots.C:434
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
monster::nRes
int nRes
Definition: monster.hh:99
netcluster::cData
std::vector< clusterdata > cData
Definition: netcluster.hh:373
xtalk::CC
float CC[4]
Definition: monster.hh:12
mmax
double mmax
Definition: CreateListEBBH.C:19
l_high
int l_high
Definition: test_config1.C:41
itime
double * itime
Definition: cwb_condor_create_ced.C:53
WDM_PREC
#define WDM_PREC
Definition: PCA_Benchmark.C:41
WSeries::getSample
DataType_t getSample(int n, double m)
Definition: wseries.hh:167
ctime
double ctime
Definition: TestCreateFrList.C:32
netpixel::rate
float rate
Definition: netpixel.hh:95
wavearray::resize
virtual void resize(unsigned int)
Definition: wavearray.cc:445
MDC_WNB
Definition: mdc.hh:127
WSeries::sizeZero
size_t sizeZero()
Definition: wseries.hh:126
y
wavearray< double > y
Definition: Test10.C:31
grec
gwavearray< double > * grec
Definition: PCA_Benchmark.C:120
sig
wavearray< double > sig
Definition: PCA_Benchmark.C:119
WSeries::maxLayer
int maxLayer()
Definition: wseries.hh:121
gwavearray::Draw
void Draw(GWAT_DRAW type=GWAT_TIME, TString options="ALP", Color_t color=kBlack)
Definition: gwavearray.cc:71
pixdata::a_90
double a_90
Definition: netpixel.hh:21
netcluster::nTofF
std::vector< vector_int > nTofF
Definition: netcluster.hh:386
exit
exit(0)
WAVE_LENGTH
#define WAVE_LENGTH
Definition: PCA_Benchmark.C:29