Logo coherent WaveBurst  
Library Reference Guide
Logo
DrawNRI2.C
Go to the documentation of this file.
1 //
2 // Draw Network Response Index skymap
3 // Author : Gabriele Vedovato
4 
5 #define L1_ENABLED
6 #define H1_ENABLED
7 #define V1_ENABLED
8 #define H2_ENABLED
9 //#define G1_ENABLED
10 //#define T1_ENABLED
11 //#define A1_ENABLED
12 //#define A2_ENABLED
13 //#define J1_ENABLED
14 
15 #define H1_SIGMA 1
16 #define L1_SIGMA 1
17 #define G1_SIGMA 1
18 #define V1_SIGMA 1
19 #define T1_SIGMA 1
20 //#define H2_SIGMA 2
21 #define H2_SIGMA 1
22 #define A1_SIGMA 1
23 #define A2_SIGMA 1
24 #define J1_SIGMA 1
25 
26 
27 //#define NTRIES 1000000
28 #define NTRIES 100000
29 
30 //#define USE_NOISE
31 
32 //#define WRITE_PLOT
33 #define PLOT_POSTFIX "_WM"
34 
35 #define SNR 10
36 
37 #define RESOLUTION 1
38 //#define RESOLUTION 2
39 //#define RESOLUTION 4
40 
41 
42 //#define WRONG_DIRECTION
43 
44 
45 //#define COORDINATES "cWB"
46 #define COORDINATES "Geographic"
47 
48 #define PROJECTION ""
49 //#define PROJECTION "hammer"
50 //#define PROJECTION "sinusoidal"
51 
52 #define DISPLAY_WORLD_MAP
53 
54 //#define DISPLAY_PERC_UNDER2
55 
56 //#define NET_GAMMA 0.02
57 #define NET_GAMMA 0.2
58 //#define NET_GAMMA 0.0
59 
60 
61 static inline int net9(double*, double, double*, double, double);
62 int GetSkyMapSensitivity(double*& x, double*& y, double*& z, double*& t,
63  TString& title, TString& ofileName, int& ndet);
64 
65 
66 void DrawNRI2() {
67 
68  int ndet=0;
69  double *x,*y,*z,*t;
70  TString title;
72 
73  int size = GetSkyMapSensitivity(x,y,z,t,title,ofileName,ndet);
74  if (ndet==0) {cout << "No detector !!!" << endl;exit(0);}
75 
76  TH2D* Fx = new TH2D("Fx","Fx", 360*RESOLUTION, 0, 360, 180*RESOLUTION, 0, 180);
77  TH2D* h2 = new TH2D("angle","angle", 360*RESOLUTION, 0, 360, 180*RESOLUTION, 0, 180);
78  h2->SetStats(kFALSE);
79 
80  h2->GetXaxis()->SetNdivisions(70318);
81  h2->GetXaxis()->SetLabelFont(42);
82  h2->GetXaxis()->SetLabelOffset(0.012);
83  h2->GetXaxis()->SetTitleOffset(1.1);
84  h2->GetXaxis()->SetTitleFont(72);
85  h2->GetYaxis()->SetNdivisions(409);
86  h2->GetYaxis()->SetLabelFont(42);
87  h2->GetYaxis()->SetLabelOffset(0.01);
88  h2->GetZaxis()->SetLabelFont(42);
89  h2->GetXaxis()->SetTitleFont(42);
90  h2->GetXaxis()->SetTitle("Phi");
91  h2->GetXaxis()->CenterTitle(true);
92  h2->GetYaxis()->SetTitleFont(42);
93  h2->GetYaxis()->SetTitle("Theta");
94  h2->GetYaxis()->CenterTitle(true);
95  h2->GetZaxis()->SetLabelOffset(0.00001);
96  h2->GetZaxis()->SetNoExponent(false);
97 
98  h2->SetTitle(title);
99 #ifdef DISPLAY_PERC_UNDER2
100  int ncnt_under2[360*RESOLUTION][180*RESOLUTION];
101  for(int i=0;i<360*RESOLUTION;i++) for(int j=0;j<180*RESOLUTION;j++) ncnt_under2[i][j]=0;
102 #endif
103  int ncnt[360*RESOLUTION][180*RESOLUTION];
104  for(int i=0;i<360*RESOLUTION;i++) for(int j=0;j<180*RESOLUTION;j++) ncnt[i][j]=0;
105  for (int i=0;i<size;i++) {
106  int ii=int(x[i]); if(ii>359) ii=359;
107  int jj=int(y[i]); if(jj>179) jj=179;
108  ii*=RESOLUTION;
109  jj*=RESOLUTION;
110  ncnt[ii][jj]++;
111 #ifdef DISPLAY_PERC_UNDER2
112  if(z[i]<2) ncnt_under2[ii][jj]++;
113  Fx->SetBinContent(ii+1,jj+1,t[i]);
114 #else
115  double binc = h2->GetBinContent(ii+1,jj+1);
116  h2->SetBinContent(ii+1,jj+1,binc+z[i]);
117 #endif
118  }
119  for(int i=0;i<360*RESOLUTION;i++) for(int j=0;j<180*RESOLUTION;j++) {
120 #ifdef DISPLAY_PERC_UNDER2
121  if(ncnt[i][j]>0) h2->SetBinContent(i+1,j+1,100.*(double(ncnt_under2[i][j])/double(ncnt[i][j])));
122 #else
123  double binc = h2->GetBinContent(i+1,j+1);
124  if(ncnt[i][j]>0) binc/=ncnt[i][j];
125  h2->SetBinContent(i+1,j+1,binc);
126 #endif
127  }
128  //h2->Draw("colfz");
129 
130  size=180*RESOLUTION*360*RESOLUTION;
131  double* xx = new double[size];
132  double* yy = new double[size];
133  double* zz = new double[size];
134 
135  int cnt=0;
136  for(int i=0;i<360*RESOLUTION;i++) for(int j=0;j<180*RESOLUTION;j++) {
137 #ifdef DISPLAY_PERC_UNDER2
138  double perc_under2 = h2->GetBinContent(i+1,j+1);
139  double fx = Fx->GetBinContent(i+1,j+1);
140 #else
141  double nri = h2->GetBinContent(i+1,j+1);
142 #endif
143  xx[cnt]=i;
144  yy[cnt]=j;
145 #ifdef DISPLAY_PERC_UNDER2
146  zz[cnt]=perc_under2;
147  //if(perc_under2>90) zz[cnt]=perc_under2;
148  //if(perc_under2>0) zz[cnt]=fx;
149  //if(perc_under2>0) zz[cnt]=fx*(100.-perc_under2);
150  //zz[cnt]=fx;
151 #else
152  zz[cnt]=nri;
153 #endif
154  cnt++;
155  }
156 
157  gskymap* gSM = new gskymap(int(7));
158  gSM->SetOptions(PROJECTION,COORDINATES,RESOLUTION);
159 
160 #ifdef DISPLAY_WORLD_MAP
161  gSM->SetWorldMap();
162 #endif
163  TH2D* hh2 = (TH2D*)gSM->GetHistogram();
164 // hh2->GetZaxis()->SetRangeUser(0,1);
165  gSM->SetTitle(title);
166  if(COORDINATES=="Geographic") {
167  for (int i=0;i<size;i++) {
168  CwbToGeographic(xx[i],yy[i],xx[i],yy[i]);
169  }
170  }
171  gSM->FillData(size, xx, yy, zz);
172 // gSM->SetPalette(1);
173  gSM->Draw(0);
174 
175 
176 #ifdef WRITE_PLOT
177  cout << "Print File : " << ofileName.Data() << endl;
178  canvas->Print(ofileName);
179 #endif
180 
181 }
182 
183 #define NDETECTORS 9
184 
185 int
186 GetSkyMapSensitivity(double*& x, double*& y, double*& z, double*& t,
187  TString& title, TString& ofileName, int& ndet) {
188 
189  bool detector_selected[NDETECTORS];
190 #ifdef H1_ENABLED
191  detector_selected[0]=true; // LHO1
192 #else
193  detector_selected[0]=false; // LHO1
194 #endif
195 #ifdef L1_ENABLED
196  detector_selected[1]=true; // LLO
197 #else
198  detector_selected[1]=false; // LLO
199 #endif
200 #ifdef G1_ENABLED
201  detector_selected[2]=true; // GEO
202 #else
203  detector_selected[2]=false; // GEO
204 #endif
205 #ifdef V1_ENABLED
206  detector_selected[3]=true; // VIRGO
207 #else
208  detector_selected[3]=false; // VIRGO
209 #endif
210 #ifdef T1_ENABLED
211  detector_selected[4]=true; // TAMA
212 #else
213  detector_selected[4]=false; // TAMA
214 #endif
215 #ifdef H2_ENABLED
216  detector_selected[5]=true; // LHO2
217 #else
218  detector_selected[5]=false; // LHO2
219 #endif
220 #ifdef A1_ENABLED
221  detector_selected[6]=true; // AIGO
222 #else
223  detector_selected[6]=false; // AIGO
224 #endif
225 #ifdef A2_ENABLED
226  detector_selected[7]=true; // AIGO2
227 #else
228  detector_selected[7]=false; // AIGO2
229 #endif
230 #ifdef J1_ENABLED
231  detector_selected[8]=true; // LCGT
232 #else
233  detector_selected[8]=false; // LCGT
234 #endif
235 
236  TString detector_name[NDETECTORS];
237  detector_name[0]="H1"; // LHO1
238  detector_name[1]="L1"; // LLO
239  detector_name[2]="G1"; // GEO
240  detector_name[3]="V1"; // VIRGO
241  detector_name[4]="T1"; // TAMA
242  detector_name[5]="H2"; // LHO2
243  detector_name[6]="A1"; // AIGO
244  detector_name[7]="A2"; // AIGO2
245  detector_name[8]="J1"; // LCGT
246 
247  double detector_sensitivity[NDETECTORS];
248  detector_sensitivity[0]=H1_SIGMA; // LHO1
249  detector_sensitivity[1]=L1_SIGMA; // LLO
250  detector_sensitivity[2]=G1_SIGMA; // GEO
251  detector_sensitivity[3]=V1_SIGMA; // VIRGO
252  detector_sensitivity[4]=T1_SIGMA; // TAMA
253  detector_sensitivity[5]=H2_SIGMA; // LHO2
254  detector_sensitivity[6]=A1_SIGMA; // AIGO
255  detector_sensitivity[7]=A2_SIGMA; // AIGO2
256  detector_sensitivity[8]=J1_SIGMA; // LCGT
257 
258  TString sdetectors;
259  if (detector_selected[0]) sdetectors+="LHO1 ";
260  if (detector_selected[1]) sdetectors+="LLO ";
261  if (detector_selected[2]) sdetectors+="GEO ";
262  if (detector_selected[3]) sdetectors+="VIRGO ";
263  if (detector_selected[4]) sdetectors+="TAMA ";
264  if (detector_selected[5]) sdetectors+="LHO2 ";
265  if (detector_selected[6]) sdetectors+="AIGO ";
266  if (detector_selected[7]) sdetectors+="AIGO2 ";
267  if (detector_selected[8]) sdetectors+="LCGT ";
268 
269  for (int k=0;k<NDETECTORS;k++) {
270  if (detector_selected[k]) {
271  ndet++;
272  }
273  }
274 
275  char hacc_title[256];
276 #ifdef USE_NOISE
277  sprintf(hacc_title,"%s Network Response Index (gamma=%2.2f)",sdetectors.Data(),NET_GAMMA);
278 #else
279  sprintf(hacc_title,"%s Network Response Index (gamma=%2.2f)",sdetectors.Data(),NET_GAMMA);
280 #endif
281  title = hacc_title;
282 
283 
284  char ofile_title[256];
285 #ifdef USE_NOISE
286  sprintf(ofile_title,"%sNRI%s",sdetectors.Data(),PLOT_POSTFIX);
287 #else
288  sprintf(ofile_title,"%sNRI%s",sdetectors.Data(),PLOT_POSTFIX);
289 #endif
290  if(PROJECTION=="hammer") sprintf(ofile_title,"%s_hammer.png",ofile_title);
291  else sprintf(ofile_title,"%s.png",ofile_title);
292 
293  ofileName = ofile_title;
294  ofileName.ReplaceAll(" ","_");
295 
297 
298  for(int n=0;n<9;n++) D[n] = new detector((char*)detector_name[n].Data());
299 
300  skymap sm(0.4,0,180,0,360);
301  int L = sm.size();
302 
303  double pi = TMath::Pi();
306  double X[NDETECTORS];
307  int size = NTRIES;
308  x = new double[size];
309  y = new double[size];
310  z = new double[size];
311  t = new double[size];
312  int cnt=0;
313  gRandom->SetSeed(0);
314  for (int i=0;i<NTRIES;i++) {
315  if(i%100000==0) cout << i << endl;
316  int idx = int(gRandom->Uniform(0,L));
317 
318  //double phi=sm.getPhi(idx);
319  //double theta=sm.getTheta(idx);
320  //double psi=0;
321 
322  double phi=gRandom->Uniform(0,360);
323  double theta=gRandom->Uniform(0,180);
324  double psi=gRandom->Uniform(0,180);
325 
326  double Hp=gRandom->Uniform(-1,1);
327  double Hx=gRandom->Uniform(-1,1);
328 
329  for (int k=0;k<NDETECTORS;k++) {
330  if (detector_selected[k]) {
331  F[k] = D[k]->antenna(theta,phi,psi);
332  X[k] = Hp*F[k].real()+Hx*F[k].imag();
333  }
334  }
335 
336  double gp=0;
337  double gx=0;
338  double gI=0;
339  for (int k=0;k<NDETECTORS;k++) {
340  if (detector_selected[k]) {
341  gp+=F[k].real()*F[k].real();
342  gx+=F[k].imag()*F[k].imag();
343  gI+=F[k].real()*F[k].imag();
344  }
345  }
346  double gR = (gp-gx)/2.;
347  double gr = (gp+gx)/2.;
348  double gc = sqrt(gR*gR+gI*gI); // norm of complex antenna pattern
349 
350  double Xp=0;
351  double Xx=0;
352  for (int k=0;k<NDETECTORS;k++) {
353  if (detector_selected[k]) {
354  Xp+=X[k]*F[k].real();
355  Xx+=X[k]*F[k].imag();
356  }
357  }
358  double uc = (Xp*gx - Xx*gI); // u cos of rotation to PCF
359  double us = (Xx*gp - Xp*gI); // u sin of rotation to PCF
360  double vc = (gp*uc + gI*us); // (u*f+)/|u|^2 - 'cos' for v
361  double vs = (gx*us + gI*uc); // (u*fx)/|u|^2 - 'sin' for v
362  double u[NDETECTORS];
363  double v[NDETECTORS];
364  double um=0;
365  double vm=0;
366  for (int k=0;k<NDETECTORS;k++) {
367  u[k]=0;
368  v[k]=0;
369  if (detector_selected[k]) {
370  u[k]=F[k].real()*uc+F[k].imag()*us;
371  um+=u[k]*u[k];
372  v[k]=-F[k].imag()*vc+F[k].real()*vs;
373  vm+=v[k]*v[k];
374  }
375  }
376 
377 // Compute Network Detector Index
378  double gamma=NET_GAMMA;
379  double GAMMA = 1.-gamma*gamma; // network regulator
380  double ndi=net9(u,um,v,vm,GAMMA);
381 
382  x[cnt]=phi;
383  y[cnt]=theta;
384  z[cnt]=ndi;
385  cnt++;
386  }
387  return size;
388 }
389 
390 inline int net9(double* u, double um, double* v, double vm, double g) {
391  double q = (1.-g)*um;
392 
393  return int(u[0]*u[0]>q) - int((u[0]*u[0]/um+v[0]*v[0]/vm)>g) +
394  int(u[1]*u[1]>q) - int((u[1]*u[1]/um+v[1]*v[1]/vm)>g) +
395  int(u[2]*u[2]>q) - int((u[2]*u[2]/um+v[2]*v[2]/vm)>g) +
396  int(u[3]*u[3]>q) - int((u[3]*u[3]/um+v[3]*v[3]/vm)>g) +
397  int(u[4]*u[4]>q) - int((u[4]*u[4]/um+v[4]*v[4]/vm)>g) +
398  int(u[5]*u[5]>q) - int((u[5]*u[5]/um+v[5]*v[5]/vm)>g) +
399  int(u[6]*u[6]>q) - int((u[6]*u[6]/um+v[6]*v[6]/vm)>g) +
400  int(u[7]*u[7]>q) - int((u[7]*u[7]/um+v[7]*v[7]/vm)>g) +
401  int(u[8]*u[8]>q) - int((u[8]*u[8]/um+v[8]*v[8]/vm)>g);
402 }
403 
wavearray< double > t(hp.size())
TH2D * GetHistogram()
Definition: gskymap.hh:138
gskymap * gSM
static double g(double e)
Definition: GNGen.cc:116
#define PROJECTION
Definition: DrawNRI2.C:48
#define V1_SIGMA
Definition: DrawNRI2.C:18
WSeries< float > v[nIFO]
Definition: cwb_net.C:80
int n
Definition: cwb_net.C:28
wavearray< double > z
Definition: Test10.C:32
double imag() const
Definition: wavecomplex.hh:70
TString("c")
#define L1_SIGMA
Definition: DrawNRI2.C:16
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
float theta
return wmap canvas
void Draw(int dpaletteId=1, Option_t *option="colfz")
Definition: gskymap.cc:460
Long_t size
int j
Definition: cwb_net.C:28
i drho i
int GetSkyMapSensitivity(double *&x, double *&y, double *&z, double *&t, TString &title, TString &ofileName, int &ndet)
Definition: DrawNRI2.C:186
#define T1_SIGMA
Definition: DrawNRI2.C:19
double pi
Definition: TestChirp.C:18
#define H2_SIGMA
Definition: DrawNRI2.C:21
#define NET_GAMMA
Definition: DrawNRI2.C:57
#define NDETECTORS
Definition: DrawNRI2.C:183
void FillData(int size, double *phi, double *theta, double *binc)
Definition: gskymap.cc:394
float phi
wavecomplex antenna(double, double, double=0.)
param: source theta,phi, polarization angle psi in degrees
Definition: detector.cc:490
float psi
wavearray< double > xx
Definition: TestFrame1.C:11
TGraph * gr
double real() const
Definition: wavecomplex.hh:69
i() int(T_cor *100))
#define PLOT_POSTFIX
Definition: DrawNRI2.C:33
double Pi
gwavearray< double > * gx
void DrawNRI2()
Definition: DrawNRI2.C:66
#define H1_SIGMA
Definition: DrawNRI2.C:15
double D[50000]
#define A1_SIGMA
Definition: DrawNRI2.C:22
#define G1_SIGMA
Definition: DrawNRI2.C:17
int k
#define NTRIES
Definition: DrawNRI2.C:28
double F
Definition: skymap.hh:63
wavearray< double > yy
Definition: TestFrame5.C:12
TString ofileName
Definition: MergeTrees.C:37
void SetTitle(TString title)
Definition: gskymap.hh:152
char title[256]
Definition: SSeriesExample.C:1
void SetWorldMap(bool drawWorldMap=true)
Definition: gskymap.hh:154
#define A2_SIGMA
Definition: DrawNRI2.C:23
#define RESOLUTION
Definition: DrawNRI2.C:37
void CwbToGeographic(double ilongitude, double ilatitude, double &olongitude, double &olatitude)
Definition: skycoord.hh:414
int cnt
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
#define J1_SIGMA
Definition: DrawNRI2.C:24
size_t size()
Definition: skymap.hh:136
wavearray< double > y
Definition: Test10.C:31
#define COORDINATES
Definition: DrawNRI2.C:46
void SetOptions(TString projection="hammer", TString coordinate="Geographic", double resolution=1, bool goff=false)
Definition: gskymap.cc:84
exit(0)
#define GAMMA(TYPE)
Definition: xroot.hh:5
static int net9(double *, double, double *, double, double)
Definition: DrawNRI2.C:390