Logo coherent WaveBurst  
Library Reference Guide
Logo
watplot.cc
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Sergey Klimenko, Gabriele Vedovato
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 
19 #define WATPLOT_CC
20 #include <iostream>
21 #include <stdexcept>
22 #include <fstream>
23 #include "watplot.hh"
24 #include "WDM.hh"
25 #include "TSystem.h"
26 #include "TObjArray.h"
27 #include "TObjString.h"
28 #include "TPaletteAxis.h"
29 
30 ClassImp(watplot) // used by THtml doc
31 
32 
33 //______________________________________________________________________________
34 /* Begin_Html
35 <center><h2>watlpot class</h2></center>
36 The class watlpot is used to plot different objects types (wavearray, WSeries, skymap, netcluster, ...)
37 
38 <p>
39 <h3><a name="example">Example</a></h3>
40 <p>
41 The macro <a href="./tutorials/wat/testWDM_3.C.html">testWDM_3.C</a> is an example which shown how to plot wavearray object.<br>
42 <p>
43 
44 End_Html
45 
46 Begin_Macro
47 testWDM_3.C
48 End_Macro
49 
50 Begin_Html
51 The macro <a href="./tutorials/wat/testWDM_1.C.html">testWDM_1.C</a> is an example which shown how to plot WSeries object.<br>
52 <p>
53 
54 End_Html
55 
56 Begin_Macro
57 testWDM_1.C
58 End_Macro
59 
60 Begin_Html
61 The macro <a href="./tutorials/wat/DrawSkymapWithSkyplot.C.html">DrawSkymapWithSkyplot.C</a> is an example which shown how to plot skymap object.<br>
62 <p>
63 
64 End_Html
65 
66 Begin_Macro
67 DrawSkymapWithSkyplot.C
68 End_Macro
69 
70 Begin_Html
71 The macro <a href="./tutorials/wat/DrawClusterWithSkyplot.C.html">DrawClusterWithSkyplot.C</a> is an example which shown how to plot netcluster object.<br>
72 <p>
73 
74 End_Html
75 
76 Begin_Macro
77 DrawClusterWithSkyplot.C
78 End_Macro */
79 
80 
81 // ++++++++++++++++++++++++++++++++++++++++++++++
82 // S. Klimenko, University of Florida
83 // WAT plot class
84 // ++++++++++++++++++++++++++++++++++++++++++++++
85 
86 watplot::watplot(char* name, int i1, int i2, int i3, int i4) {
87 //
88 // Default constructor
89 //
90 // Input parameters are used to create a new canvas.
91 //
92 // i1,i2 are the pixel coordinates of the top left corner of
93 // the canvas (if i1 < 0) the menubar is not shown)
94 // i3 is the canvas size in pixels along X
95 // i4 is the canvas size in pixels along Y
96 //
97 // default values are : name=NULL, i1=200, i2=20, i3=800, i4=600
98 //
99 
100  this->title = "";
101  this->xtitle = "";
102  this->ytitle = "";
103  this->ncol = 0;
104  this->opt = "";
105  this->col = 1;
106  this->t1 = 0.;
107  this->t2 = 0.;
108  this->fft = false;
109  this->f1 = 0.;
110  this->f2 = 0.;
111  this->hist2D = NULL;
112 
113  char defn[16] = "watplot";
114  if (name && gROOT->FindObject(name)!=NULL) {
115  printf("watplot : Error Canvas Name %s already exist",name);
116  exit(1);
117  }
118 
119  null();
120  if(!name) name = defn;
121  canvas= new TCanvas(name, name, i1, i2, i3, i4);
122  canvas->Clear();
123  canvas->ToggleEventStatus();
124  canvas->SetGridx();
125  canvas->SetGridy();
126  canvas->SetFillColor(kWhite);
127  canvas->SetRightMargin(0.10);
128  canvas->SetLeftMargin(0.10);
129  canvas->SetBottomMargin(0.13);
130  canvas->SetBorderMode(0);
131 
132  // remove the red box around canvas
133  gStyle->SetFrameBorderMode(0);
134  gROOT->ForceStyle();
135 
136  gStyle->SetTitleH(0.050);
137  gStyle->SetTitleW(0.95);
138  gStyle->SetTitleY(0.98);
139  gStyle->SetTitleFont(12,"D");
140  gStyle->SetTitleColor(kBlue,"D");
141  gStyle->SetTextFont(12);
142  gStyle->SetTitleFillColor(kWhite);
143  gStyle->SetLineColor(kWhite);
144  gStyle->SetNumberContours(256);
145  gStyle->SetCanvasColor(kWhite);
146  gStyle->SetStatBorderSize(1);
147 
148 }
149 
150 void watplot::plot(wavearray<double> &td, char* opt, int col,
151  double t1, double t2, bool fft, float f1, float f2, bool psd, float t3, bool oneside) {
152 //
153 // Plot wavearray<double> time series
154 //
155 // td : wavearray time series
156 // opt : TGraph::Draw options
157 // col : set TGraph line color
158 // t1 : start of time interval in seconds
159 // t2 : end of time interval in seconds
160 // fft : true -> plot fft
161 // f1 : set begin frequency (Hz)
162 // f2 : set end frequency (Hz)
163 // psd : true -> plot psd using blackmanharris window
164 // t3 : is the chunk length (sec) used to produce the psd
165 // oneside : true/false -> oneside/doubleside
166 //
167 
168  if ( t2 < t1) {
169  printf(" Plot error: t2 must be greater then t1.\n");
170  return;
171  }
172  if ( f2 < f1) {
173  printf(" Plot error: f2 must be greater then f1.\n");
174  return;
175  }
176 
177  if(t2==0.) { t1 = td.start(); t2 = t1+td.size()/td.rate(); }
178 
179  double Ts = (td.rate() == 0)? 1.: 1./td.rate();
180  int i1 = int((t1-td.start())/Ts);
181  int i2 = (t2 == 0.)? td.size() : int((t2-td.start())/Ts);
182  //int nmax = i2-i1;
183  int nmax = (i2-i1)-(i2-i1)%2; // nmax even
184  wavearray<double> _x(nmax);
185  wavearray<double> _y(nmax);
186  double xmin, xmax;
187  double ymin, ymax;
188  xmin=0.;
189  xmax=(nmax-1)*Ts;
190  ymin=1.e40;
191  ymax=0.;
192  double x0 = i1*Ts+td.start();
193  TGraph* G;
194 
195  for (int i=0; i < nmax; i++) {
196  _x.data[i] = x0 + Ts*i - 0.;
197  _y.data[i] = td.data[i+i1];
198  if (ymin > _y.data[i]) ymin=_y.data[i];
199  if (ymax < _y.data[i]) ymax=_y.data[i];
200  }
201 
202  if(fft) {
203  if(f2==0.) { f1 = 0.; f2 = td.rate()/2.; }
204  if(psd) { // power spectrum density
205  t3 = t3<0. ? 0. : t3;
206  t3 = t3<(_y.size()/td.rate()) ? t3 : (_y.size()/td.rate());
207  int bsize = t3*td.rate();
208  bsize-=bsize%2;
209  int loops = int(_y.size()/bsize)==0 ? 1 : int(_y.size()/bsize);
210  double Fs=(double)td.rate()/(double)bsize;
211  wavearray<double> w(bsize);
212  wavearray<double> u(bsize);
213  wavearray<double> psd(bsize); psd=0;
214  blackmanharris(w);
215  for (int n=0;n<loops;n++) {
216  int shift=n*bsize;
217  //cout << "shift: " << shift << endl;
218  for (int i=0;i<bsize;i++) u[i]=_y[i+shift];
219  for (int i=0;i<bsize;i++) u[i]*=w[i];
220  u.FFTW(1);
221  for (int i=0;i<bsize;i+=2) psd[i/2]+=pow(u[i],2)+pow(u[i+1],2);
222  }
223  for (int i=0;i<bsize/2; i++) psd[i]=sqrt(psd[i]/(double)loops)*sqrt(1./Fs); // double side spectra 1/sqrt(Hz);
224  if(oneside) psd*=sqrt(2.); // one side spectra *sqrt(2);
225  _y=psd;nmax=bsize;
226  _x.resize(nmax);
227  for (int i=0;i<nmax/2;i++) _x.data[i]=i*Fs;
228  } else {
229  wavearray<double> _z(nmax);
230  double Fs=(double)td.rate()/(double)_y.size();
231  _y.FFTW(1);
232  for (int i=0;i<nmax;i+=2) _z.data[i/2]=sqrt(pow(_y.data[i],2)+pow(_y.data[i+1],2));
233  for (int i=0;i<nmax/2;i++) _y.data[i]=_z.data[i]*(1./Fs); // double side spectra 1/Hz
234  for (int i=0;i<nmax/2;i++) _x.data[i]=i*Fs;
235  if(oneside) for(int i=0;i<nmax/2;i++) _y.data[i]*=sqrt(2.); // one side spectra 1/Hz
236  }
237  nmax/=2;
238  ymin=1.e40;
239  ymax=0.;
240  for (int i=0;i<nmax;i++) {
241  if ((_x.data[i]>f1)&&(_x.data[i]<f2)) {
242  if (ymin > _y.data[i]) ymin=_y.data[i];
243  if (ymax < _y.data[i]) ymax=_y.data[i];
244  }
245  }
246  }
247 
248  G = new TGraph(nmax,_x.data,_y.data);
249  G->SetLineColor(col);
250  G->SetMarkerColor(col);
251 // G->Draw(opt);
252 // canvas->Update();
253  G->GetHistogram()->SetTitle("");
254  G->SetFillColor(kWhite);
255 
256  //G->GetXaxis()->SetNdivisions(70318);
257  G->GetXaxis()->SetTitleFont(42);
258  G->GetXaxis()->SetLabelFont(42);
259  G->GetXaxis()->SetLabelOffset(0.012);
260  G->GetXaxis()->SetTitleOffset(1.5);
261 
262  //G->GetYaxis()->SetNdivisions(508);
263  G->GetYaxis()->SetTitleFont(42);
264  G->GetYaxis()->SetLabelFont(42);
265  G->GetYaxis()->SetLabelOffset(0.01);
266  G->GetYaxis()->SetTitleOffset(1.4);
267 //
268  if(fft) {
269  G->GetHistogram()->GetXaxis()->SetRangeUser(f1,f2);
270  //G->GetHistogram()->GetYaxis()->SetRangeUser(ymin/2.,ymax*10.);
271  G->GetHistogram()->GetYaxis()->SetRangeUser(ymin/2.,ymax*1.1);
272  G->GetHistogram()->SetXTitle("freq, hz");
273  } else {
274  G->GetHistogram()->SetXTitle("time, s");
275  }
276  G->GetHistogram()->SetYTitle("magnitude");
277  G->Draw(opt);
278  graph.push_back(G);
279 // graph->Fit("fit","VMR");
280 
281  return;
282 }
283 
284 void watplot::plsmooth(WSeries<double> &w, int sfact, double t1, double t2, char* opt, int pal, int palId)
285 {
286 
287  t1 = t1==0. ? w.start() : t1;
288  t2 = t2==0. ? w.start()+w.size()/w.rate() : t2;
289 
290  int ni = 1<<w.pWavelet->m_Level;
291  int nb = int((t1-w.start())*w.rate()/ni);
292  int nj = int((t2-t1)*w.rate())/ni;
293  int ne = nb+nj;
294  double rate=w.rate();
295  double freq=rate/2/ni;
296  rate = rate/ni;
297 
298  double** iw = new double*[nj];
299  for(int j=0;j<nj;j++) iw[j] = new double[ni];
300 
301 // double wmax=0.;
303  for(int i=0;i<ni;i++) {
304  w.getLayer(wl,i);
305 // for(int j=nb;j<ne;j++) iw[j-nb][i] = fabs(wl.data[j]);
306 
307  for(int j=nb;j<ne;j++) {
308 //printf("DEB14 %d %d,%d,%d %d\n",j,j-nb,i,nj,wl.size());
309  iw[j-nb][i] = wl.data[j]*wl.data[j];
310 // if(wmax<iw[j-nb][i]) wmax=iw[j-nb][i];
311  }
312 
313  }
314 
315  int smfact=1;
316  for(int n=0;n<sfact;n++) smfact*=3;
317  double** sw = NULL;
318  int sni=ni;
319  int snj=nj;
320 
321  if(sfact==0) {
322  sw = new double*[snj];
323  for(int j=0;j<snj;j++) sw[j] = new double[sni];
324  for(int i=0;i<sni;i++) {
325  for(int j=0;j<snj;j++) {
326  sw[j][i] = iw[j][i];
327  }
328  }
329  }
330 
331  for(int n=1;n<=sfact;n++) {
332  if(sw!=NULL) {
333  for(int j=0;j<snj;j++) delete [] sw[j];
334  delete [] sw;
335  }
336  sni*=3;
337  snj*=3;
338  sw = new double*[snj];
339  for(int j=0;j<snj;j++) sw[j] = new double[sni];
340  for(int i=0;i<sni;i++) for(int j=0;j<snj;j++) sw[j][i] = 0.;
341  for(int i=1;i<sni-1;i++) {
342  for(int j=1;j<snj-1;j++) {
343  for(int j3=j-1;j3<=j+1;j3++) {
344  for(int i3=i-1;i3<=i+1;i3++) {
345 //printf("%d %d %d %d\n",i,j,i3/3,j3/3);
346  sw[j][i] += iw[(j3-j3%3)/3][(i3-i3%3)/3];
347  }
348  }
349  sw[j][i]/=9.;
350  }
351  }
352  for(int j=0;j<nj;j++) delete [] iw[j];
353  delete [] iw;
354  iw = new double*[snj];
355  for(int j=0;j<snj;j++) iw[j] = new double[sni];
356  for(int i=0;i<sni;i++) {
357  for(int j=0;j<snj;j++) {
358  iw[j][i] = sw[j][i];
359  }
360  }
361  }
362 
363  for(int j=0;j<snj;j++) delete [] iw[j];
364  delete [] iw;
365 
366  if(hist2D) { delete hist2D; hist2D=NULL; }
367  //hist2D=new TH2F("WTF","", snj, t1-w.start(), t2-w.start(), sni, 0., freq*ni);
368  hist2D=new TH2F("WTF","", snj, 0., t2-t1, sni, 0., freq*ni);
369  hist2D->SetXTitle("time, sec");
370  hist2D->SetYTitle("frequency, Hz");
371 
372  Int_t colors[30]={101,12,114,13,115,14,117,15,16,17,166,18,19,
373  167,0,0,167,19,18,166,17,16,15,117,14,115,13,114,12,101};
374  if(pal==0) gStyle->SetPalette(30,colors);
375  else if(pal<0) {SetPlotStyle(palId,pal);}
376  else {gStyle->SetPalette(1,0);gStyle->SetNumberContours(pal);}
377 /*
378  double swmax=0.;
379  for(int i=0;i<sni;i++) {
380  for(int j=0;j<snj;j++) {
381  if(swmax<sw[j][i]) swmax=sw[j][i];
382  }
383  }
384 */
385 
386  int snb=smfact*nb;
387  int sne=smfact*ne;
388  double srate=smfact*rate;
389  double sfreq=freq/smfact;
390  for(int i=0;i<sni;i++) {
391  for(int j=snb;j<sne;j++) {
392  double x = sw[j-snb][i];
393  //double x = sw[j-snb][i]*wmax/swmax;
394  //double x = sqrt(sw[j-snb][i]);
395  //double x = sqrt(sw[j-snb][i]*wmax/swmax);
396  //hist2D->Fill(double(j)/srate,(i+0.5)*sfreq,x);
397  hist2D->Fill(double(j-snb)/srate,(i+0.5)*sfreq,x);
398  }
399  }
400 
401  for(int j=0;j<snj;j++) delete [] sw[j];
402  delete [] sw;
403 
404  hist2D->SetStats(kFALSE);
405  hist2D->SetFillColor(kWhite);
406  hist2D->GetXaxis()->SetTitleFont(42);
407  hist2D->GetXaxis()->SetLabelFont(42);
408  hist2D->GetXaxis()->SetLabelOffset(0.012);
409  hist2D->GetXaxis()->SetTitleOffset(1.1);
410  hist2D->GetYaxis()->SetTitleFont(42);
411  hist2D->GetYaxis()->SetLabelFont(42);
412  hist2D->GetYaxis()->SetLabelOffset(0.01);
413  hist2D->GetZaxis()->SetLabelFont(42);
414  hist2D->SetTitleOffset(1.3,"Y");
415  if(opt) hist2D->Draw(opt);
416  else hist2D->Draw("COLZ");
417 
418  // change palette's width
419  canvas->Update();
420  TPaletteAxis *palette = (TPaletteAxis*)hist2D->GetListOfFunctions()->FindObject("palette");
421  palette->SetX1NDC(0.91);
422  palette->SetX2NDC(0.933);
423  palette->SetTitleOffset(0.92);
424  palette->GetAxis()->SetTickSize(0.01);
425  canvas->Modified();
426 
427  return;
428 }
429 
430 
431 double watplot::getmax(WSeries<double> &w, double t1, double t2)
432 {
433 //
434 // return max value of WSeries w object
435 //
436 // t1 : start of time interval in seconds
437 // t2 : end of time interval in seconds
438 //
439 
440  float x;
441 
442  t1 = t1==0. ? w.start() : t1;
443  t2 = t2==0. ? w.start()+w.size()/w.rate() : t2;
444 
445  int ni = 1<<w.pWavelet->m_Level;
446  int nb = int((t1-w.start())*w.rate()/ni);
447  int nj = int((t2-t1)*w.rate())/ni;
448  int ne = nb+nj;
449  double rate=w.rate();
450  rate = rate/ni;
451 
453  double xmax;
454 
455  for(int i=0;i<ni;i++) {
456  w.getLayer(wl,i);
457  for(int j=nb;j<ne;j++) {
458  x=fabs(wl.data[j]);
459  if(x>xmax) xmax=x;
460  }
461  }
462  return xmax;
463 }
464 
465 void watplot::plot(WSeries<double> &w, int mode, double t1, double t2, char* opt, int pal, int palId)
466 {
467 //
468 // Plot TF series
469 //
470 // w : WSeries<double> object
471 // mode : plot type
472 // if WSeries is WDM
473 // 0 : sqrt((E00+E90)/2)
474 // 1 : (E00+E90)/2
475 // 2 : sqrt((E00+E90)/2)
476 // 3 : amplitude:00
477 // 4 : energy:00
478 // 5 : |amplitude:00|
479 // 6 : amplitude:90
480 // 7 : energy:90
481 // 8 : |amplitude:90|
482 // if WSeries is wavelet
483 // 0 : amplitude
484 // 1 : energy
485 // 2 : |amplitude|
486 //
487 // t1 : start of time interval in seconds
488 // t2 : end of time interval in seconds
489 // opt : root draw options
490 // pal : palette
491 //
492 
493  t1 = t1==0. ? w.start() : t1;
494  t2 = t2==0. ? w.start()+w.size()/w.rate() : t2;
495 
496  if(t2 <= t1) {
497  printf("watplot::plot error: t2 must be greater then t1.\n");
498  exit(1);
499  }
500  if((t1 < w.start()) || (t1 > w.start()+w.size()/w.rate())) {
501  printf("watplot::plot error: t1 must be in this range [%0.12g,%0.12g]\n",
502  w.start(),w.start()+w.size()/w.rate());
503  exit(1);
504  }
505  if((t2 < w.start()) || (t2 > w.start()+w.size()/w.rate())) {
506  printf("watplot::plot error: t2 must be in this range [%0.12g,%0.12g]\n",
507  w.start(),w.start()+w.size()/w.rate());
508  exit(1);
509  }
510 
511  Int_t colors[30]={101,12,114,13,115,14,117,15,16,17,166,18,19,
512  167,0,0,167,19,18,166,17,16,15,117,14,115,13,114,12,101};
513  if(pal==0) gStyle->SetPalette(30,colors);
514  else if(pal<0) {SetPlotStyle(palId,pal);}
515  else {gStyle->SetPalette(1,0);gStyle->SetNumberContours(pal);}
516 
517  float x;
518  TString ztitle("");
519 
520  if(w.isWDM()) {
521 
522  double A00=1;
523  double A90=1;
524 
525  if(mode==0) mode=2;
526  if(mode==1) {ztitle="(E00+E90)/2";}
527  if(mode==2) {ztitle="sqrt((E00+E90)/2)";}
528 
529  if(mode==3) {ztitle="amplitude:00";A90=0;}
530  if(mode==4) {ztitle="energy:00";A00=sqrt(2.);A90=0;}
531  if(mode==5) {ztitle="|amplitude:00|";A90=0;}
532 
533  if(mode==6) {ztitle="amplitude:90";A00=0;}
534  if(mode==7) {ztitle="energy:90";A00=0;A90=sqrt(2.);}
535  if(mode==8) {ztitle="|amplitude:90|";A00=0;}
536 
537  mode*=-1;
538 
540  int M = w.getLevel();
541  double* map00 = wdm->pWWS;
542  double tsRate = w.rate();
543  int mF = int(w.size()/wdm->nSTS);
544  int nTC = w.size()/(M+1)/mF; // # of Time Coefficients
545  double* map90 = map00 + (mF-1)*(M+1)*nTC;
546 
547  //printf("nTC = %d, rate = %d M = %d\n", nTC, (int)w.rate(), M);
548 
549  // make Y bins:
550  double* yBins = new double[M+2];
551  double dF = tsRate/M/2.;
552  yBins[0] = 0;
553  yBins[1] = dF/2;
554  for(int i=2; i<=M; ++i) yBins[i] = yBins[1] + (i-1)*dF;
555  yBins[M+1] = tsRate/2.;
556 
557  const double scale = 1./w.wrate();
558  if(hist2D) { delete hist2D; hist2D=NULL; }
559  hist2D=new TH2F("WTF", "", 2*nTC, 0., nTC*scale, M+1, yBins);
560  hist2D->SetXTitle("time, sec");
561  hist2D->SetYTitle("frequency, Hz");
562  hist2D->GetXaxis()->SetRangeUser(t1-w.start(),t2-w.start());
563 
564  double v;
565  int it1 = (t1-w.start())/scale + 1;
566  int it2 = (t2-w.start())/scale + 1;
567  if(it2<=it1 || it2>nTC)it2 = nTC;
568 
569  map00+=it1*(M+1); map90+=it1*(M+1);
570  for(int i=it1; i<it2; ++i){
571  if(i){
572  v = procOpt(mode%3, A00*map00[0], A90*map90[0]); //first half-band
573  hist2D->SetBinContent( 2*i , 1, v);
574  hist2D->SetBinContent( 2*i+1 , 1, v);
575 
576  v = procOpt(mode%3, A00*map00[M], A90*map90[M]); //last half-band
577  hist2D->SetBinContent( 2*i , M+1, v);
578  hist2D->SetBinContent( 2*i+1 , M+1, v);
579  }
580  for(int j=1; j<M; ++j){
581  v = procOpt(mode%3, A00*map00[j], A90*map90[j]);
582  hist2D->SetBinContent( 2*i , j+1, v);
583  hist2D->SetBinContent( 2*i+1 , j+1, v);
584  }
585  map00+=M+1; map90+=M+1;
586  }
587 
588  delete [] yBins;
589 
590  } else {
591 
592  if(mode==0) ztitle="amplitude";
593  if(mode==1) ztitle="energy";
594  if(mode==2) ztitle="|amplitude|";
595 
596  int ni = 1<<w.pWavelet->m_Level;
597  int nb = int((t1-w.start())*w.rate()/ni);
598  int nj = int((t2-t1)*w.rate())/ni;
599  int ne = nb+nj;
600  double rate=w.rate();
601  double freq=rate/2/ni;
602  rate = rate/ni;
603 
604 // cout<<rate<<endl;
605 
606  if(hist2D) { delete hist2D; hist2D=NULL; }
607  hist2D=new TH2F("WTF","", nj, t1-w.start(), t2-w.start(), ni, 0., freq*ni);
608  hist2D->SetXTitle("time, sec");
609  hist2D->SetYTitle("frequency, Hz");
610 
612  double avr,rms;
613 
614  if(mode==0)
615  for(int i=0;i<ni;i++) {
616  w.getLayer(wl,i);
617  for(int j=nb;j<=ne;j++) {
618  x=wl.data[j];
619  hist2D->Fill((j+0.5)/rate,(i+0.5)*freq,x);
620  }
621  }
622 
623  if(mode==1)
624  for(int i=0;i<ni;i++) {
625  w.getLayer(wl,i);
626  for(int j=nb;j<=ne;j++) {
627  x=wl.data[j]*wl.data[j];
628  hist2D->Fill((j+0.5)/rate,(i+0.5)*freq,x);
629  }
630  }
631 
632  if(mode==2)
633  for(int i=0;i<ni;i++) {
634  w.getLayer(wl,i);
635  for(int j=nb;j<=ne;j++) {
636  x = fabs(wl.data[j]);
637  hist2D->Fill((j+0.5)/rate,(i+0.5)*freq,x);
638  }
639  }
640 
641  if(mode==4)
642  for(int i=0;i<ni;i++) {
643  w.getLayer(wl,i);
644  for(int j=nb;j<=ne;j++) {
645  x=(wl.data[j] > 0.) ? 1. : -1.;
646  hist2D->Fill((j+0.5)/rate,(i+0.5)*freq,x);
647  }
648  }
649 
650  if(mode==5)
651  for(int i=0;i<ni;i++) {
652  w.getLayer(wl,i);
653  for(int j=nb;j<=ne;j++) {
654  x=(wl.data[j] > 0.) ? 1. : -1.;
655  if(wl.data[j] == 0.) x=0.;
656  hist2D->Fill((j+0.5)/rate,(i+0.5)*freq,x);
657  }
658  }
659 
660  double sum = 0.;
661  int nsum = 0;
662 
663  if(mode==3){
664  for(int i=0;i<ni;i++) {
665  w.getLayer(wl,i);
666  wl.getStatistics(avr,rms);
667  for(int j=nb;j<ne-0;j++) {
668  x=(wl.data[j]-avr)/rms;
669  x*=x;
670  if(x>1){ sum += x; nsum++; }
671  hist2D->Fill(j/rate,(i+0.5)*freq,x);
672  }
673  }
674  }
675  }
676 
677  hist2D->SetStats(kFALSE);
678  hist2D->SetTitleFont(12);
679  hist2D->SetFillColor(kWhite);
680 
681  char title[256];
682  sprintf(title,"Scalogram (%s)",ztitle.Data());
683  hist2D->SetTitle(title);
684 
685  hist2D->GetXaxis()->SetNdivisions(506);
686  hist2D->GetXaxis()->SetLabelFont(42);
687  hist2D->GetXaxis()->SetLabelOffset(0.014);
688  hist2D->GetXaxis()->SetTitleOffset(1.4);
689  hist2D->GetYaxis()->SetTitleOffset(1.2);
690  hist2D->GetYaxis()->SetNdivisions(506);
691  hist2D->GetYaxis()->SetLabelFont(42);
692  hist2D->GetYaxis()->SetLabelOffset(0.01);
693  hist2D->GetZaxis()->SetLabelFont(42);
694  hist2D->GetZaxis()->SetNoExponent(false);
695  hist2D->GetZaxis()->SetNdivisions(506);
696 
697  hist2D->GetXaxis()->SetTitleFont(42);
698  hist2D->GetXaxis()->SetTitle("Time (sec)");
699  hist2D->GetXaxis()->CenterTitle(true);
700  hist2D->GetYaxis()->SetTitleFont(42);
701  hist2D->GetYaxis()->SetTitle("Frequency (Hz)");
702  hist2D->GetYaxis()->CenterTitle(true);
703 
704  hist2D->GetZaxis()->SetTitleOffset(0.6);
705  hist2D->GetZaxis()->SetTitleFont(42);
706  //hist2D->GetZaxis()->SetTitle(ztitle);
707  hist2D->GetZaxis()->CenterTitle(true);
708 
709  hist2D->GetXaxis()->SetLabelSize(0.03);
710  hist2D->GetYaxis()->SetLabelSize(0.03);
711  hist2D->GetZaxis()->SetLabelSize(0.03);
712 
713 
714  if(opt) hist2D->Draw(opt);
715  else hist2D->Draw("COLZ");
716 
717  // change palette's width
718  canvas->Update();
719  TPaletteAxis *palette = (TPaletteAxis*)hist2D->GetListOfFunctions()->FindObject("palette");
720  palette->SetX1NDC(0.91);
721  palette->SetX2NDC(0.933);
722  palette->SetTitleOffset(0.92);
723  palette->GetAxis()->SetTickSize(0.01);
724  canvas->Modified();
725 
726  return;
727 }
728 
729 void watplot::plot(skymap &sm, char* opt, int pal) {
730 //
731 // Draw skymap
732 //
733 // sm : skymap object
734 // opt : root draw options
735 // pal : palette
736 //
737 
738  int ni = sm.size(0); // number of theta layers
739  int nj = 0; // number of phi collumns
740 
741  for(int i=1; i<=ni; i++) if(nj<int(sm.size(i))) nj = sm.size(i);
742 
743  double t1 = sm.theta_1;
744  double t2 = sm.theta_2;
745  double p1 = sm.phi_1;
746  double p2 = sm.phi_2;
747  double dt = ni>1 ? (t2-t1)/(ni-1) : 0.;
748  double dp = nj>0 ? (p2-p1)/nj : 0.;
749 
750  if(hist2D) { delete hist2D; hist2D=NULL; }
751  hist2D=new TH2F("WTS","", nj,p1,p2, ni,90.-t2,90.-t1);
752  hist2D->SetXTitle("phi, deg.");
753  hist2D->SetYTitle("theta, deg.");
754 
755  Int_t colors[30]={101,12,114,13,115,14,117,15,16,17,166,18,19,
756  167,0,0,167,19,18,166,17,16,15,117,14,115,13,114,12,101};
757  if(pal==0) gStyle->SetPalette(30,colors);
758  else {gStyle->SetPalette(1,0);gStyle->SetNumberContours(pal);}
759 
761  double theta, phi;
762 
763  for(int i=0; i<ni; i++) {
764  theta = (t2+t1)/2.+(i-ni/2)*dt;
765  for(int j=0; j<nj; j++) {
766  phi = (p2+p1)/2.+(j-nj/2)*dp;
767  hist2D->Fill(phi,90.-theta,sm.get(theta,phi));
768  }
769  }
770 
771  hist2D->SetStats(kFALSE);
772  hist2D->SetFillColor(kWhite);
773 
774  //hist2D->GetXaxis()->SetNdivisions(70318);
775  hist2D->GetXaxis()->SetTitleFont(42);
776  hist2D->GetXaxis()->SetLabelFont(42);
777  hist2D->GetXaxis()->SetLabelOffset(0.012);
778  hist2D->GetXaxis()->SetTitleOffset(1.1);
779 
780  //hist2D->GetYaxis()->SetNdivisions(508);
781  hist2D->GetYaxis()->SetTitleFont(42);
782  hist2D->GetYaxis()->SetLabelFont(42);
783  hist2D->GetYaxis()->SetLabelOffset(0.01);
784 
785  hist2D->GetZaxis()->SetLabelFont(42);
786 //
787  if(opt) hist2D->Draw(opt);
788  else hist2D->Draw("COLZ");
789 
790  // change palette's width
791  canvas->Update();
792  TPaletteAxis *palette = (TPaletteAxis*)hist2D->GetListOfFunctions()->FindObject("palette");
793  palette->SetX1NDC(0.91);
794  palette->SetX2NDC(0.933);
795  palette->SetTitleOffset(0.92);
796  palette->GetAxis()->SetTickSize(0.01);
797  canvas->Modified();
798 
799  return;
800 }
801 
802 void watplot::plot(netcluster* pwc, int cid, int nifo, char type, int irate, char* opt, int pal, bool wp) {
803 //
804 // monster event display (only for 2G analysis)
805 // display pixels of all resolution levels
806 //
807 // pwc : pointer to netcluster object
808 // cid : cluster id
809 // nifo : number of detectors
810 // type : 'L' -> plot event likelihood pixels, 'N' -> plot event null pixels
811 // irate : select pixel rate
812 // 0 : all rates
813 // 1 : optimal rate
814 // x : rate x
815 // opt : root draw options
816 // pal : palette
817 // wp : true -> wavelet packet
818 //
819 
820  bool isPCs = std::isupper(type); // are Principal Components ?
821  type = std::toupper(type);
822 
823  if(type != 'L' && type != 'N') return;
824 
825  double RATE = pwc->rate; // original rate
826 
827  std::vector<int>* vint = &(pwc->cList[cid-1]); // pixel list
828 
829  int V = vint->size(); // cluster size
830  if(!V) return;
831 
832  // check if WDM (2G), else exit
833  netpixel* pix = pwc->getPixel(cid,0);
834  int mp= size_t(pwc->rate/pix->rate+0.1);
835  int mm= pix->layers; // number of wavelet layers
836  if(mm==mp) { // wavelet
837  cout << "watplot::plot - Error : monster event display is enabled only for WDM (2G)" << endl;
838  exit(1);
839  }
840 
841  bool optrate=false;
842  if(irate==1) { // extract optimal rate
843  optrate=true;
844  vector_int* pv = pwc->cRate.size() ? &(pwc->cRate[cid-1]) : NULL;
845  irate = pv!=NULL ? (*pv)[0] : 0;
846  }
847 
848  int minLayers=1000;
849  int maxLayers=0;
850  double minTime=1e20;
851  double maxTime=0.;
852  double minFreq=1e20;
853  double maxFreq=0.;
854  for(int j=0; j<V; j++) { // loop over the pixels
855  netpixel* pix = pwc->getPixel(cid,j);
856  if(!pix->core) continue;
857 
858  if((irate)&&(irate != int(pix->rate+0.5))) continue; // if irate!=0 skip rate!=irate
859 
860  if(pix->layers<minLayers) minLayers=pix->layers;
861  if(pix->layers>maxLayers) maxLayers=pix->layers;
862 
863  double dt = 1./pix->rate;
864  double time = int(pix->time/pix->layers)/double(pix->rate); // central bin time
865  time -= dt/2.; // begin bin time
866  if(time<minTime) minTime=time;
867  if(time+dt>maxTime) maxTime=time+dt;
868 
869  double freq = pix->frequency*pix->rate/2.;
870  if(freq<minFreq) minFreq=freq;
871  if(freq>maxFreq) maxFreq=freq;
872  }
873 
874  int minRate=RATE/(maxLayers-1);
875  int maxRate=RATE/(minLayers-1);
876 
877  double mindt = 1./maxRate;
878  double maxdt = 1./minRate;
879  double mindf = minRate/2.;
880  double maxdf = maxRate/2.;
881 
882  //cout << "minRate : " << minRate << "\t\t\t maxRate : " << maxRate << endl;
883  //cout << "minTime : " << minTime << "\t\t\t maxTime : " << maxTime << endl;
884  //cout << "minFreq : " << minFreq << "\t\t\t maxFreq : " << maxFreq << endl;
885  //cout << "mindt : " << mindt << "\t\t\t maxdt : " << maxdt << endl;
886  //cout << "mindf : " << mindf << "\t\t\t maxdf : " << maxdf << endl;
887 
888  double iminTime = minTime-maxdt;
889  double imaxTime = maxTime+maxdt;
890  int nTime = (imaxTime-iminTime)*maxRate;
891 
892  if(hist2D) { delete hist2D; hist2D=NULL; }
893  hist2D=new TH2F("WTF", "WTF", nTime, iminTime, imaxTime, 2*(maxLayers-1), 0, RATE/2);
894  hist2D->SetXTitle("time, sec");
895  hist2D->SetYTitle("frequency, Hz");
896 
897  Int_t colors[30]={101,12,114,13,115,14,117,15,16,17,166,18,19,
898  167,0,0,167,19,18,166,17,16,15,117,14,115,13,114,12,101};
899  if(pal==0) gStyle->SetPalette(30,colors);
900  else {gStyle->SetPalette(1,0);gStyle->SetNumberContours(pal);}
901 
902  hist2D->SetStats(kFALSE);
903  hist2D->SetTitleFont(12);
904  hist2D->SetFillColor(kWhite);
905 
906  hist2D->GetXaxis()->SetNdivisions(506);
907  hist2D->GetXaxis()->SetLabelFont(42);
908  hist2D->GetXaxis()->SetLabelOffset(0.014);
909  hist2D->GetXaxis()->SetTitleOffset(1.4);
910  hist2D->GetYaxis()->SetTitleOffset(1.2);
911  hist2D->GetYaxis()->SetNdivisions(506);
912  hist2D->GetYaxis()->SetLabelFont(42);
913  hist2D->GetYaxis()->SetLabelOffset(0.01);
914  hist2D->GetZaxis()->SetLabelFont(42);
915  hist2D->GetZaxis()->SetNoExponent(false);
916  hist2D->GetZaxis()->SetNdivisions(506);
917 
918  hist2D->GetXaxis()->SetTitleFont(42);
919  hist2D->GetXaxis()->SetTitle("Time (sec)");
920  hist2D->GetXaxis()->CenterTitle(true);
921  hist2D->GetYaxis()->SetTitleFont(42);
922  hist2D->GetYaxis()->SetTitle("Frequency (Hz)");
923  hist2D->GetYaxis()->CenterTitle(true);
924 
925  hist2D->GetZaxis()->SetTitleOffset(0.6);
926  hist2D->GetZaxis()->SetTitleFont(42);
927  hist2D->GetZaxis()->CenterTitle(true);
928 
929  hist2D->GetXaxis()->SetLabelSize(0.03);
930  hist2D->GetYaxis()->SetLabelSize(0.03);
931  hist2D->GetZaxis()->SetLabelSize(0.03);
932 
933  double dFreq = (maxFreq-minFreq)/10.>2*maxdf ? (maxFreq-minFreq)/10. : 2*maxdf ;
934  double mFreq = minFreq-dFreq<0 ? 0 : minFreq-dFreq;
935  double MFreq = maxFreq+dFreq>RATE/2 ? RATE/2 : maxFreq+dFreq;
936  hist2D->GetYaxis()->SetRangeUser(mFreq, MFreq);
937 
938  double dTime = (maxTime-minTime)/10.>2*maxdt ? (maxTime-minTime)/10. : 2*maxdt ;
939  double mTime = minTime-dTime<iminTime ? iminTime : minTime-dTime;
940  double MTime = maxTime+dTime>imaxTime ? imaxTime : maxTime+dTime;
941  hist2D->GetXaxis()->SetRangeUser(mTime,MTime);
942 
943  int npix=0;
944  double Null=0;
945  double Likelihood=0;
946  for(int n=0; n<V; n++) {
947  netpixel* pix = pwc->getPixel(cid,n);
948  if(!pix->core) continue;
949 
950  if((irate)&&(irate != int(pix->rate+0.5))) continue; // if irate!=0 skip rate!=irate
951 
952  double like=0;
953  double null=0;
954  if(wp) { // likelihoodWP
955  like = pix->likelihood>0. ? pix->likelihood : 0.;
956  null = pix->null>0. ? pix->null : 0.;
957  } else { // likelihood2G
958  double sSNR=0;
959  double wSNR=0;
960  for(int m=0; m<nifo; m++) {
961  sSNR += pow(pix->getdata('S',m),2); // snr whitened reconstructed signal 00
962  sSNR += pow(pix->getdata('P',m),2); // snr whitened reconstructed signal 90
963  wSNR += pow(pix->getdata('W',m),2); // snr whitened at the detector 00
964  wSNR += pow(pix->getdata('U',m),2); // snr whitened at the detector 90
965  }
966  if(!isPCs) {sSNR/=2;wSNR/=2;} // if not principal components we use (00+90)/2
967  like = sSNR;
968  null = wSNR-sSNR;
969  }
970 
971  int iRATE = int(pix->rate+0.5);
972  int M=maxRate/iRATE;
973  int K=2*(maxLayers-1)/(pix->layers-1);
974  double dt = 1./pix->rate;
975  double itime = int(pix->time/pix->layers)/double(pix->rate); // central bin time
976  itime -= dt/2.; // begin bin time
977  int i=(itime-iminTime)*maxRate;
978  int j=pix->frequency*K;
979  if(iRATE!=irate && irate!=0) continue;
980  Null+=null;
981  Likelihood+=like;
982  int L=0;int R=1;while (R < iRATE) {R*=2;L++;}
983  for(int m=0;m<M;m++) {
984  for(int k=0;k<K;k++) {
985  if(null<0) null=0;
986  double A = hist2D->GetBinContent(i+1+m,j+1+k-K/2);
987  if(type=='L') hist2D->SetBinContent(i+1+m,j+1+k-K/2,like+A);
988  else hist2D->SetBinContent(i+1+m,j+1+k-K/2,null+A);
989  }
990  }
991 
992  if(type=='L' && like>0) npix++;
993  if(type=='N' && null>0) npix++;
994  }
995 
996  char gtitle[256];
997  if(type=='L') sprintf(gtitle,"Likelihood %3.0f - dt(ms) [%.6g:%.6g] - df(hz) [%.6g:%.6g] - npix %d",
998  Likelihood,1000*mindt,1000*maxdt,mindf,maxdf,npix);
999  else sprintf(gtitle,"Null %3.0f - dt(ms) [%.6g:%.6g] - df(hz) [%.6g:%.6g] - npix %d",
1000  Null,1000*mindt,1000*maxdt,mindf,maxdf,npix);
1001 
1002  hist2D->SetTitle(gtitle);
1003 
1004  //cout << "Event Likelihood : " << Likelihood << endl;
1005  //cout << "Event Null : " << Null << endl;
1006 
1007  if(opt) hist2D->Draw(opt);
1008  else hist2D->Draw("COLZ");
1009 
1010  // change palette's width
1011  canvas->Update();
1012  TPaletteAxis *palette = (TPaletteAxis*)hist2D->GetListOfFunctions()->FindObject("palette");
1013  palette->SetX1NDC(0.91);
1014  palette->SetX2NDC(0.933);
1015  palette->SetTitleOffset(0.92);
1016  palette->GetAxis()->SetTickSize(0.01);
1017  canvas->Modified();
1018 
1019  return;
1020 }
1021 
1022 void watplot::plot(clusterdata* pcd, double inj_mchirp) {
1023 //
1024 // chirp event display (only for 2G analysis) : WARNING : Add check 2G !!!
1025 // display frequency vs time
1026 //
1027 // pcd : pointer to clusterdata structure
1028 // inj_mchirp : injected mchirp value
1029 //
1030 
1031  static TGraphErrors egraph;
1032 
1033  double x,y;
1034  double ex,ey;
1035  int np = pcd->chirp.GetN();
1036  egraph.Set(np);
1037  for(int i=0;i<np;i++) {
1038  pcd->chirp.GetPoint(i,x,y);
1039  ex=pcd->chirp.GetErrorX(i);
1040  ey=pcd->chirp.GetErrorY(i);
1041  egraph.SetPoint(i,x,y);
1042  egraph.SetPointError(i,ex,ey);
1043  }
1044  TF1* fit = &(pcd->fit);
1045 
1046  char title[256];
1047  if(inj_mchirp) {
1048  sprintf(title,"chirp mass : rec = %3.3f [%3.2f] inj = %3.3f , chi2 = %3.2f",
1049  pcd->mchirp,pcd->mchirperr,inj_mchirp,pcd->chi2chirp);
1050  } else {
1051  sprintf(title,"chirp mass : rec = %3.3f [%3.2f] , chi2 = %3.2f",
1052  pcd->mchirp,pcd->mchirperr,pcd->chi2chirp);
1053  }
1054  egraph.SetTitle(title);
1055  egraph.GetHistogram()->SetStats(kFALSE);
1056  egraph.GetHistogram()->SetTitleFont(12);
1057  egraph.SetFillColor(kWhite);
1058  egraph.SetLineColor(kBlack);
1059  egraph.GetXaxis()->SetNdivisions(506);
1060  egraph.GetXaxis()->SetLabelFont(42);
1061  egraph.GetXaxis()->SetLabelOffset(0.014);
1062  egraph.GetXaxis()->SetTitleOffset(1.4);
1063  egraph.GetYaxis()->SetTitleOffset(1.2);
1064  egraph.GetYaxis()->SetNdivisions(506);
1065  egraph.GetYaxis()->SetLabelFont(42);
1066  egraph.GetYaxis()->SetLabelOffset(0.01);
1067  egraph.GetXaxis()->SetTitleFont(42);
1068  egraph.GetXaxis()->SetTitle("Time (sec)");
1069  egraph.GetXaxis()->CenterTitle(true);
1070  egraph.GetYaxis()->SetTitleFont(42);
1071  egraph.GetYaxis()->SetTitle("Frequency^{-8/3}");
1072  egraph.GetYaxis()->CenterTitle(true);
1073  egraph.GetXaxis()->SetLabelSize(0.03);
1074  egraph.GetYaxis()->SetLabelSize(0.03);
1075 
1076  egraph.Draw("AP");
1077  fit->Draw("same");
1078 
1079  return;
1080 }
1081 
1082 void watplot::SetPlotStyle(int paletteId, int NCont) {
1083 //
1084 // set palette type
1085 //
1086 
1087  const Int_t NRGBs = 5;
1088 // const Int_t NCont = 255;
1089  NCont=abs(NCont);
1090 
1091  if (fabs(paletteId)==1) {
1092  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
1093  Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
1094  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
1095  Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
1096  if (paletteId<0) {
1097  TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
1098  } else {
1099  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
1100  }
1101  } else
1102  if (fabs(paletteId)==2) {
1103  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
1104  Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
1105  //Double_t red[NRGBs] = { 0.00, 0.00, 0.00, 1.00, 1.00 };
1106  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
1107  //Double_t green[NRGBs] = { 0.00, 1.00, 1.00, 1.00, 0.00 };
1108  Double_t blue[NRGBs] = { 1.00, 1.00, 0.00, 0.00, 0.00 };
1109  if (paletteId<0) {
1110  TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
1111  } else {
1112  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
1113  }
1114  } else
1115  if (fabs(paletteId)==3) {
1116  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
1117  Double_t red[NRGBs] = { 0.00, 0.09, 0.18, 0.09, 0.00 };
1118  Double_t green[NRGBs] = { 0.01, 0.02, 0.39, 0.68, 0.97 };
1119  Double_t blue[NRGBs] = { 0.17, 0.39, 0.62, 0.79, 0.97 };
1120  if (paletteId<0) {
1121  TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
1122  } else {
1123  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
1124  }
1125  } else
1126  if (fabs(paletteId)==4) {
1127  Double_t stops[NRGBs] = { 0.00, 0.50, 0.75, 0.875, 1.00 };
1128  Double_t red[NRGBs] = { 1.00, 1.00, 1.00, 1.00, 1.00 };
1129  Double_t green[NRGBs] = { 1.00, 0.75, 0.50, 0.25, 0.00 };
1130  Double_t blue[NRGBs] = { 0.00, 0.00, 0.00, 0.00, 0.00 };
1131  if (paletteId<0) {
1132  TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
1133  } else {
1134  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
1135  }
1136  } else
1137  if (fabs(paletteId)==5) { // Greyscale palette
1138  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
1139  Double_t red[NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 };
1140  Double_t green[NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 };
1141  Double_t blue[NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 };
1142  if (paletteId<0) {
1143  TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
1144  } else {
1145  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
1146  }
1147  }
1148  gStyle->SetNumberContours(NCont);
1149 
1150  return;
1151 }
1152 
1153 // wavescan (2D histogram)
1154 /*
1155 void watplot::scan(wavearray<double>& x, int M1, int M2, double=t1, double=t2, char* =c, int=nc)
1156 {
1157  int i,j,k;
1158  WDM<double>* wdm[10];
1159 
1160  if(M2-M1 > 10) M2 = M1+9;
1161  for(i=M1; i<=M2; i++) wdm[i-M1] = new WDM<double>(i,i);
1162 
1163 
1164  for(i=M1; i<=M2; i++) delete wdm[i-M1];
1165 
1166 }
1167 */
1168 
1169 
1170 void
1172 //
1173 // blackmanharris window
1174 //
1175 // w : in/out wavearray - must be initialized with w.size() & w.rate()
1176 //
1177 
1178  int size = (int)w.size();
1179 
1180  for (int i = 0; i < size; i++)
1181  {
1182  double f = 0.0;
1183  f = ((double) i) / ((double) (size - 1));
1184  w[i] = 0.35875 -
1185  0.48829 * cos(2.0 * TMath::Pi() * f) +
1186  0.14128 * cos(4.0 * TMath::Pi() * f) -
1187  0.01168 * cos(6.0 * TMath::Pi() * f);
1188  }
1189 
1190  double norm = 0;
1191  for (int i=0;i<size;i++) norm += pow(w[i],2);
1192  norm /= size;
1193  for (int i=0;i<size;i++) w[i] /= sqrt(norm);
1194 }
1195 
1196 void
1198 //
1199 // Save plot with name -> fname
1200 //
1201 
1202  if(canvas!=NULL && fname!="") {
1203  if(fname.Contains(".png")) {
1204  TString gname(fname);
1205  gname.ReplaceAll(".png",".gif");
1206  canvas->Print(gname);
1207  char cmd[1024];
1208  sprintf(cmd,"convert %s %s",gname.Data(),fname.Data());
1209  //cout << cmd << endl;
1210  gSystem->Exec(cmd);
1211  sprintf(cmd,"rm %s",gname.Data());
1212  //cout << cmd << endl;
1213  gSystem->Exec(cmd);
1214  } else {
1215  canvas->Print(fname);
1216  }
1217  }
1218 }
1219 
1220 void
1221 watplot::goptions(char* opt, int col, double t1, double t2, bool fft, float f1, float f2, bool psd, float t3, bool oneside) {
1222 //
1223 // Set graphical options
1224 //
1225 // opt : TGraph::Draw options
1226 // col : set TGraph line color
1227 // t1 : start of time interval in seconds
1228 // t2 : end of time interval in seconds
1229 // fft : true -> plot fft
1230 // f1 : set begin frequency (Hz)
1231 // f2 : set end frequency (Hz)
1232 // psd : true -> plot psd using blackmanharris window
1233 // t3 : is the chunk length (sec) used to produce the psd
1234 // oneside : true/false -> oneside/doubleside
1235 //
1236 
1237  this->ncol= 0;
1238  this->opt = opt;
1239  this->col = col;
1240  this->t1 = t1;
1241  this->t2 = t2;
1242  this->fft = fft;
1243  this->f1 = f1;
1244  this->f2 = f2;
1245  this->psd = psd;
1246  this->t3 = t3;
1247  this->oneside = oneside;
1248 
1249  for(size_t i=0; i<graph.size(); i++) {
1250  if(graph[i]) delete graph[i];
1251  }
1252  graph.clear();
1253 }
1254 
1255 void
1257 //
1258 // Set TGraph titles
1259 //
1260 // title : graph title
1261 // xtitle : x axis name
1262 // ytitle : y axis name
1263 //
1264 
1265  this->title=title;
1266  this->xtitle=xtitle;
1267  this->ytitle=ytitle;
1268  if(graph.size()>0) {
1269  if(title!="") graph[0]->SetTitle(title);
1270  if(xtitle!="") graph[0]->GetHistogram()->SetXTitle(xtitle);
1271  if(ytitle!="") graph[0]->GetHistogram()->SetYTitle(ytitle);
1272  canvas->Update();
1273  }
1274 }
1275 
1277 //
1278 // assignement operator
1279 //
1280 // graph >> x : fill x wavearray with watplot object graph with values
1281 //
1282 
1283  x=graph.data;
1284  return x;
1285 }
1286 
1288 //
1289 // assignement operator
1290 //
1291 // x >> graph : fill watplot object graph with x wavearray values
1292 //
1293 
1294  if(x.size()==0) {
1295  cout << "watplot::operator(>>) : input array with size=0" << endl;
1296  exit(1);
1297  }
1298 
1299  graph.opt.ToUpper();
1300  bool logx = false; if(graph.opt.Contains("LOGX")) {logx=true;graph.opt.ReplaceAll("LOGX","");}
1301  bool logy = false; if(graph.opt.Contains("LOGY")) {logy=true;graph.opt.ReplaceAll("LOGY","");}
1302 
1303  double t1 = graph.t1==0. ? x.start() : graph.t1;
1304  double t2 = graph.t2==0. ? x.start()+x.size()/x.rate() : graph.t2;
1305  double f1 = graph.f1==0. ? 0. : graph.f1;
1306  double f2 = graph.f2==0. ? x.rate()/2. : graph.f2;
1307 
1308  TString opt = graph.ncol ? "SAME" : graph.opt;
1309 
1310  graph.data=x;
1311  if(logx) graph.canvas->SetLogx();
1312  if(logy) graph.canvas->SetLogy();
1313  graph.plot(x, const_cast<char*>(opt.Data()), graph.col+graph.ncol, t1, t2, graph.fft, f1, f2, graph.psd, graph.t3, graph.oneside);
1314  if(graph.graph.size()>0) {
1315  if(graph.title!="") graph.graph[0]->SetTitle(graph.title);
1316  if(graph.xtitle!="") graph.graph[0]->GetHistogram()->SetXTitle(graph.xtitle);
1317  if(graph.ytitle!="") graph.graph[0]->GetHistogram()->SetYTitle(graph.ytitle);
1318  graph.canvas->Update();
1319  }
1320  graph.ncol++;
1321  return graph;
1322 }
1323 
1325 //
1326 // print operator
1327 //
1328 // graph >> fname : save watplot graph to fname file
1329 //
1330 
1331  if(graph.graph.size()==0) {
1332  cout << "watplot::operator (>>) - Error : No graph in the object" << endl;
1333  return fname;
1334  }
1335  TString name = fname;
1336  name.ReplaceAll(" ","");
1337 
1338  TObjArray* token = TString(name).Tokenize(TString("."));
1339  TObjString* ext_tok = (TObjString*)token->At(token->GetEntries()-1);
1340  TString ext = ext_tok->GetString();
1341  if(ext=="txt") {
1342 
1343  ofstream out;
1344  out.open(name.Data(),ios::out);
1345  if (!out.good()) {cout << "watplot::operator (>>) - Error : Opening File : " << name.Data() << endl;exit(1);}
1346 
1347  double x,y;
1348  int size=graph.graph[0]->GetN();
1349  for (int j=0;j<size;j++) {
1350  graph.graph[0]->GetPoint(j,x,y);
1351  if(x>graph.f1 && x<graph.f2) out << x << "\t" << y << endl;
1352  }
1353 
1354  out.close();
1355  } else {
1356  graph.print(name);
1357  }
1358  return fname;
1359 }
1360 
1362 //
1363 // print operator
1364 //
1365 // graph >> fname : save watplot graph to fname file
1366 //
1367 
1368  TString name = fname;
1369  graph >> name;
1370  return fname;
1371 }
std::vector< int > vector_int
Definition: cluster.hh:35
double f1
Definition: watplot.hh:207
std::vector< vector_int > cRate
Definition: netcluster.hh:398
void gtitle(TString title="", TString xtitle="", TString ytitle="")
Definition: watplot.cc:1256
char xtitle[1024]
TH1 * t1
Double_t green[nRGBs]
double M
Definition: DrawEBHH.C:13
bool oneside
Definition: watplot.hh:211
virtual void rate(double r)
Definition: wavearray.hh:141
int irate
#define np
TString opt
size_t frequency
Definition: netpixel.hh:111
float likelihood
Definition: netpixel.hh:114
WSeries< float > v[nIFO]
Definition: cwb_net.C:80
par [0] name
int n
Definition: cwb_net.C:28
TString("c")
int palette
Definition: DrawGnetwork2.C:17
double t3
Definition: watplot.hh:210
float mchirp
Definition: netcluster.hh:84
ofstream out
Definition: cwb_merge.C:214
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
TString title
Definition: watplot.hh:198
float theta
WDM< double > wdm(nLAYERS, nLAYERS, 6, 10)
wavearray< double > psd(33)
netpixel pix(nifo)
netcluster * pwc
Definition: cwb_job_obj.C:38
std::vector< TGraph * > graph
Definition: watplot.hh:194
return wmap canvas
Double_t blue[nRGBs]
double procOpt(int opt, double val00, double val90=0.)
Definition: watplot.hh:175
double theta_1
Definition: skymap.hh:321
Long_t size
Color_t colors[16]
size_t layers
Definition: netpixel.hh:112
int m
Definition: cwb_net.C:28
DataType_t * pWWS
Definition: WaveDWT.hh:141
std::vector< vector_int > cList
Definition: netcluster.hh:397
virtual void start(double s)
Definition: wavearray.hh:137
int j
Definition: cwb_net.C:28
i drho i
double rate
Definition: netcluster.hh:378
void plot(wavearray< double > &, char *=NULL, int=1, double=0., double=0., bool=false, float=0., float=0., bool=false, float=0., bool=false)
Definition: watplot.cc:150
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
Definition: ComputeSNR.C:75
bool psd
Definition: watplot.hh:209
bool core
Definition: netpixel.hh:120
double theta_2
Definition: skymap.hh:322
TH2F * hist2D
Definition: watplot.hh:193
size_t mode
watplot p2(const_cast< char *>("TFMap2"))
wavearray< double > w
Definition: Test1.C:27
void wrate(double r)
Definition: wseries.hh:120
virtual size_t size() const
Definition: wavearray.hh:145
TCanvas * canvas
Definition: watplot.hh:192
int getLevel()
Definition: wseries.hh:109
float phi
Double_t stops[nRGBs]
bool fft
Definition: watplot.hh:206
wavearray< double > freq
Definition: Regression_H1.C:79
double G
Definition: DrawEBHH.C:12
unsigned long nSTS
Definition: WaveDWT.hh:143
TF1 * f2
i() int(T_cor *100))
bool isWDM()
Definition: wseries.hh:208
double Pi
watplot p1(const_cast< char *>("TFMap1"))
double t2
Definition: watplot.hh:205
#define NCont
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)
char fname[1024]
double getmax(WSeries< double > &w, double t1, double t2)
Definition: watplot.cc:431
void goptions(char *opt=NULL, int col=1, double t1=0., double t2=0., bool fft=false, float f1=0., float f2=0., bool psd=false, float t3=0., bool oneside=false)
Definition: watplot.cc:1221
int k
TString xtitle
Definition: watplot.hh:199
static double A
Definition: geodesics.cc:26
TObjArray * token
float chi2chirp
Definition: netcluster.hh:86
void plsmooth(WSeries< double > &, int=0, double=0., double=0., char *=NULL, int=256, int=0)
Definition: watplot.cc:284
int m_Level
Definition: Wavelet.hh:115
size_t time
Definition: netpixel.hh:110
Definition: skymap.hh:63
#define NRGBs
int col
Definition: watplot.hh:203
double f2
Definition: watplot.hh:208
netpixel * getPixel(size_t n, size_t i)
Definition: netcluster.hh:413
int npix
wavearray< double > data
Definition: watplot.hh:197
double dt
virtual void FFTW(int=1)
Definition: wavearray.cc:896
#define RATE
float mchirperr
Definition: netcluster.hh:85
int nifo
char title[256]
Definition: SSeriesExample.C:1
TGraphErrors chirp
chirp fit parameters (don&#39;t remove ! fix crash when exit from CINT)
Definition: netcluster.hh:93
void print(TString fname)
Definition: watplot.cc:1197
TString ytitle
Definition: watplot.hh:200
double fabs(const Complex &x)
Definition: numpy.cc:55
char cmd[1024]
void blackmanharris(wavearray< double > &w)
Definition: watplot.cc:1171
void null()
Definition: watplot.hh:71
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
plot hist2D
void SetPlotStyle(int paletteId, int NCont=255)
Definition: watplot.cc:1082
double phi_2
Definition: skymap.hh:324
DataType_t * data
Definition: wavearray.hh:319
double getStatistics(double &mean, double &rms) const
Definition: wavearray.cc:2128
double * itime
double dFreq
Definition: DrawPSD.C:17
double phi_1
Definition: skymap.hh:323
double get(size_t i)
param: sky index
Definition: skymap.cc:699
WaveDWT< DataType_t > * pWavelet
Definition: wseries.hh:456
double t1
Definition: watplot.hh:204
float rate
Definition: netpixel.hh:113
size_t size()
Definition: skymap.hh:136
double getdata(char type='R', size_t n=0)
Definition: netpixel.hh:74
double shift[NIFO_MAX]
virtual void resize(unsigned int)
Definition: wavearray.cc:463
wavearray< double > y
Definition: Test10.C:31
wavearray< double > & operator>>(watplot &graph, wavearray< double > &x)
Definition: watplot.cc:1276
Double_t x0
float null
Definition: netpixel.hh:115
Double_t red[nRGBs]
float wSNR[PE_MAX_EVENTS][NIFO_MAX]
Definition: cwb_pereport.C:336
TString opt
Definition: watplot.hh:202
int ncol
Definition: watplot.hh:201
exit(0)
double avr