Logo coherent WaveBurst  
Library Reference Guide 		Logo
DrawGWGC.C
Go to the documentation of this file.
1 //
2 // Draw GWGC Distribution
3 // Author : Gabriele Vedovato
4 
5 //#define RESOLUTION 2
6 #define RESOLUTION 4
7 
8 #define GWGCCatalog "$CWB_GWAT/data/GWGCCatalog_Rev1d8.txt"
9 #define DISTANCE_THR 50 // Mpc
10 
11 // celestial coordinates
12 #define SGRA_NAME "SgrA*"
13 #define SGRA_DEC -29.11667
14 #define SGRA_RA 266.41667
15 
16 #define NGC0224_NAME "M31"
17 #define NGC0224_DEC 41.2687
18 #define NGC0224_RA 10.6846
19 
20 #define NGC0292_NAME "SMC"
21 #define NGC0292_DEC -72.8002
22 #define NGC0292_RA 13.1583
23 
24 #define ESO056_115_NAME "LMC"
25 #define ESO056_115_DEC -69.7561
26 #define ESO056_115_RA 80.8941
27 
28 #define COORDINATES "Geographic"
29 #define PROJECTION ""
30 //#define PROJECTION "hammer"
31 
32 //#define WRITE_PLOT
33 
34 #define NMAX 2136800
35 
36 #define SKYMASK_FILE_NAME "SkyMaskCC_GWGC50MPC_Rev1d7_R0d4.txt"
37 //#define SKYMASK_FILE_NAME "SkyMaskCC_BRST50MPC_S6_Rev1d7_R0d4.txt"
38 
39 //#define SKYMASK_FILE_NAME "SkyMaskCC_BRST50MPC_S6_Rev1d7_R0d4_PROB.txt"
40 #define SKYMASK_RESOLUTION 0.4
41 
42 //#define SKYMASK_FILE_NAME "SkyMaskCC_BRST50MPC_S6_Rev1d7_R0d4_PHPAD.txt"
43 //#define SKYMASK_RESOLUTION 0.4
44 
45 //#define SKYMASK_FILE_NAME "SkyMaskCC_BRST50MPC_S6_Rev1d7_R0d2_PHPAD.txt"
46 //#define SKYMASK_RESOLUTION 0.2
47 
48 //#define SKYMASK_FILE_NAME "SkyMaskCC_BRST50MPC_S6_Rev1d7_R0d1_PHPAD.txt"
49 //#define SKYMASK_RESOLUTION 0.1
50 
51 #define REJECTED_SKY_PIXEL_PERCENTAGE 1.0
52 //#define ADD_PHI_PADDING
53 //#define OUTPUT_PROBABILITY_SKYMAP
54 
55 int ReadGWGCCatalog(double*& x, double*& y, double*& z, int& nGWGC);
56 int ReadBRST50MPC_S6(double*& x, double*& y, double*& z, int& nGWGC);
57 int ReadBRST50MPC_S6_ToSkymask(skymap& sm);
58 
59 void DrawGWGC() {
60 
61  double *x,*y,*z;
62 
63  x = new double[NMAX];
64  y = new double[NMAX];
65  z = new double[NMAX];
66 
67  int nGWGC=0;
68  int size = ReadGWGCCatalog(x,y,z,nGWGC);
69  //int size = ReadBRST50MPC_S6(x,y,z,nGWGC);
70 
71 #ifdef SKYMASK_FILE_NAME
72 
73  skymap sm(SKYMASK_RESOLUTION,0,180,0,360);
74  int L = sm.size();
75  //nGWGC=ReadGWGCCatalogToSkymap(sm);
76  nGWGC=ReadBRST50MPC_S6_ToSkymask(sm);
77 
78 #ifdef ADD_PHI_PADDING
79  skymap sm2(sm);
80  for (int l=0;l<L;l++) {
81  if(sm.get(l)) {
82  double th = sm.getTheta(l);
83  double ph = sm.getPhi(l);
84  double sp = sm.getPhiStep(l);
85 
86  int ll = sm.getSkyIndex(th,ph-sp); // neighbour sky index
87  int lr = sm.getSkyIndex(th,ph+sp); // neighbour sky index
88 
89  sm2.set(ll,sm.get(ll)+1);
90  sm2.set(lr,sm.get(lr)+1);
91  }
92  }
93  for (int l=0;l<L;l++) sm.set(l,sm2.get(l));
94 #endif
95 
96  ofstream out;
97  out.open(SKYMASK_FILE_NAME, ios::out);
98  if (!out.good()) {cout << "Error Opening File : " << SKYMASK_FILE_NAME << endl;exit(1);}
99  cout << "Write File : " << SKYMASK_FILE_NAME << endl;
100 #ifdef OUTPUT_PROBABILITY_SKYMAP
101  double norm=0;
102  for (int l=0;l<L;l++) norm+=sm.get(l);
103  for (int l=0;l<L;l++) out << l << " " << sm.get(l)/norm << endl;
104  //for (int l=0;l<L;l++) out << l << " " << sm.get(l) << endl;
105  out.close();
106 #else
107  double perc=0;
108  for (int l=0;l<L;l++) {
109  int mask = sm.get(l)>0 ? 1 : 0;
110  out << l << " " << mask << endl;
111  perc+=mask;
112  }
113  out.close();
114  perc/=L;
115  cout << "SkyMaskCC perc : " << perc*100 << endl;
116 #endif
117 
118 #endif
119 
120  char title[256];
121  sprintf(title,"Statistics of sources in the GWGC catalog out to a %d Mpc measured distance (# %d)",DISTANCE_THR,nGWGC);
122  //sprintf(title,"Statistics of sources in the BRST50MPC_S6 MDC SET out to a %d Mpc measured distance (# %d)",DISTANCE_THR,nGWGC);
123 
124  gskymap* gSM = new gskymap(SKYMASK_RESOLUTION,0,180,0,360);
125  gSM->SetOptions(PROJECTION,COORDINATES,RESOLUTION/2);
126  gSM->SetTitle(title);
127  gSM->SetGridxColor(kWhite);
128  gSM->SetGridyColor(kWhite);
129  gSM->SetGalacticDisk(0);
130  gSM->SetGalacticDiskColor(kYellow);
131  gSM->FillData(size, x, y, z);
132  //gSM->FillData(sm);
133  gSM->SetZaxisTitle("Mpc");
134  gSM->Draw(-2);
135 
136  gSM->DrawMarker(SGRA_RA-180, SGRA_DEC, 4, 2.0, kYellow);
137  gSM->DrawText(SGRA_RA-5-180, SGRA_DEC+5, SGRA_NAME, 0.04, kYellow);
138 /*
139  gSM->DrawMarker(NGC0292_RA-180, NGC0292_DEC, 4, 2.0, kYellow);
140  gSM->DrawText(NGC0292_RA-10-180, NGC0292_DEC+10, NGC0292_NAME, 0.04, kYellow);
141 
142  gSM->DrawMarker(ESO056_115_RA-180, ESO056_115_DEC, 4, 2.0, kYellow);
143  gSM->DrawText(ESO056_115_RA-0-180, ESO056_115_DEC+10, ESO056_115_NAME, 0.04, kYellow);
144 
145  gSM->DrawMarker(NGC0224_RA-180, NGC0224_DEC, 4, 2.0, kYellow);
146  gSM->DrawText(NGC0224_RA+5-180, NGC0224_DEC-10, NGC0224_NAME, 0.04, kYellow);
147 */
148 #ifdef WRITE_PLOT
149  char ofileName[256]="GWGCCatalog_Rev1d8_plot.png";
150 
151 #ifdef SKYMASK_FILE_NAME
152  TString ofileLabel = gSystem->ExpandPathName(SKYMASK_FILE_NAME);
153  ofileLabel.ReplaceAll(".txt","");
154  ofileLabel.ReplaceAll(".","");
155  ofileLabel.ReplaceAll("/","");
156  sprintf(ofileName,"%s_Distance%dMpc%s.PNG",ofileLabel.Data(),DISTANCE_THR,PROJECTION);
157  cout << "ofileName : " << ofileName << endl;
158 #endif
159 
160  cout << "Write : " << ofileName << endl;
161  gSM->Print(ofileName);
162 #endif
163 
164  delete [] x;
165  delete [] y;
166  delete [] z;
167 }
168 
169 int
170 ReadGWGCCatalogToSkymap(skymap& sm) {
171 
172  ifstream in;
173  in.open(gSystem->ExpandPathName(GWGCCatalog),ios::in);
174  if (!in.good()) {cout << "Error Opening File : " << gSystem->ExpandPathName(GWGCCatalog) << endl;exit(1);}
175 
176  char iline[1024];
177  in.getline(iline,1024); // ski first line (header)
178  nGWGC=0;
179  while (1) {
180 
181  in.getline(iline,1024);
182  if (!in.good()) break;
183  //cout << iline << endl;
184  TObjArray* tok = TString(iline).Tokenize(TString('|'));
185 
186  TObjString* tname = (TObjString*)tok->At(1);
187  TObjString* tra = (TObjString*)tok->At(2);
188  TObjString* tdec = (TObjString*)tok->At(3);
189  TObjString* tdist = (TObjString*)tok->At(14);
190 
191  TString name = tname->GetString();
192  double ra = tra->GetString().Atof();
193  double dec = tdec->GetString().Atof();
194  double dist = tdist->GetString().Atof();
195  if (dist<DISTANCE_THR) {
196 
197  double th = -(dec-90);
198  double ph = ra*360./24.;
199  int ind = sm.getSkyIndex(th,ph);
200  sm.set(ind,sm.get(ind)+1);
201 
202  nGWGC++;
203  }
204  }
205  in.close();
206 
207  return nGWGC;
208 }
209 
210 
211 int
212 ReadGWGCCatalog(double*& x, double*& y, double*& z, int& nGWGC) {
213 
214  ifstream in;
215  in.open(gSystem->ExpandPathName(GWGCCatalog),ios::in);
216  if (!in.good()) {cout << "Error Opening File : " << gSystem->ExpandPathName(GWGCCatalog) << endl;exit(1);}
217 
218  char iline[1024];
219  in.getline(iline,1024); // ski first line (header)
220  nGWGC=0;
221  while (1) {
222 
223  in.getline(iline,1024);
224  if (!in.good()) break;
225  //cout << iline << endl;
226  TObjArray* tok = TString(iline).Tokenize(TString('|'));
227 
228  TObjString* tname = (TObjString*)tok->At(1);
229  TObjString* tra = (TObjString*)tok->At(2);
230  TObjString* tdec = (TObjString*)tok->At(3);
231  TObjString* tdist = (TObjString*)tok->At(14);
232 
233  TString name = tname->GetString();
234  double ra = tra->GetString().Atof();
235  double dec = tdec->GetString().Atof();
236  double dist = tdist->GetString().Atof();
237  if (dist<DISTANCE_THR) {
238 
239  //Geographic coordinates
240  double ph = ra*360./24.+180.;
241  double th = dec;
242 //if(ph<0) cout << "ph " << ph << endl;
243 //if(th<0) cout << "th " << th << endl;
244  //double ph = ra*360./24.;
245  //double th = -(dec-90);
246  //double th = -(-dec-90);
247 
248  x[nGWGC]=ph;
249  y[nGWGC]=th;
250  z[nGWGC]=dist;
251  nGWGC++;
252 if (name.CompareTo("NGC0224")==0) { // M31
253 //if (name.CompareTo("NGC0292")==0) { // SMC
254 //if (name.CompareTo("ESO056-115")==0) { // LMC
255 //if (fabs(dist-0.05)<0.01) {
256  cout << name.Data() << " ra : " << ra << " dec : " << dec << " " << dist << endl;
257  cout << "PH : " << ph << " TH : " << th << " " << dist << endl;
258 }
259  }
260  }
261  in.close();
262  cout << "nGWGC : " << nGWGC << endl;
263 
264  int size=nGWGC;
265  for (int i=0;i<360*RESOLUTION;i++) {
266  for (int j=0;j<180*RESOLUTION;j++) {
267  double ph = i/(double)RESOLUTION;
268  double th = j/(double)RESOLUTION;
269  x[size]=ph;
270  y[size]=th;
271  z[size]=DISTANCE_THR;
272  size++;
273  }
274  }
275 
276  Int_t *index = new Int_t[size];
277  TMath::Sort(size,z,index,true);
278  double T[NMAX];
279  for (int i=0;i<size;i++) T[i]=x[index[i]];
280  for (int i=0;i<size;i++) x[i]=T[i];
281  for (int i=0;i<size;i++) T[i]=y[index[i]];
282  for (int i=0;i<size;i++) y[i]=T[i];
283  for (int i=0;i<size;i++) T[i]=z[index[i]];
284  for (int i=0;i<size;i++) z[i]=T[i];
285  delete [] index;
286 
287  cout << "size : " << size << endl;
288  return size;
289 }
290 
291 
292 int
293 ReadBRST50MPC_S6_ToSkymask(skymap& sm) {
294 
295  double Pi = TMath::Pi();
296 
297  TFile *ifile = TFile::Open(SKYMASK_FILE_NAME);
298  TTree* itree = (TTree *) gROOT->FindObject("MDC");
299  int isize = itree->GetEntries();
300  //itree->Draw("15*2.5e-21/SimHrss:External_x:External_phi:External_psi","","goff");
301  itree->Draw("15*2.5e-21/SimHrss:External_x:External_phi:EarthCtrGPS","","goff");
302  double* idistance = itree->GetV1();
303  double* itheta = itree->GetV2();
304  double* iphi = itree->GetV3();
305  //double* ipsi = itree->GetV4();
306  double* igps = itree->GetV4();
307 
308  int nGWGC=0;
309  for(int n=0;n<isize;n++) {
310 
311  double dist = idistance[n];
312 
313  if (dist<DISTANCE_THR) {
314 
315  double th = 0;
316  th = acos(itheta[n]);
317  th*= 180/Pi;
318 
319  double ph = 0;
320  ph = iphi[n] > 0 ? iphi[n] : 2*Pi+iphi[n];
321  ph*= 180/Pi;
322  ph = sm.phi2RA(ph, igps[n]);
323 
324  int ind = sm.getSkyIndex(th,ph);
325  sm.set(ind,sm.get(ind)+1);
326 
327  nGWGC++;
328  }
329  }
330 
331  return nGWGC;
332 }
333 
334 
335 int
336 ReadBRST50MPC_S6(double*& x, double*& y, double*& z, int& nGWGC) {
337 
338  double Pi = TMath::Pi();
339 
340  x = new double[NMAX];
341  y = new double[NMAX];
342  z = new double[NMAX];
343 
344  TFile *ifile = TFile::Open(SKYMASK_FILE_NAME);
345  TTree* itree = (TTree *) gROOT->FindObject("MDC");
346  int isize = itree->GetEntries();
347  //itree->Draw("15*2.5e-21/SimHrss:External_x:External_phi:External_psi","","goff");
348  itree->Draw("15*2.5e-21/SimHrss:External_x:External_phi:EarthCtrGPS","","goff");
349  double* idistance = itree->GetV1();
350  double* itheta = itree->GetV2();
351  double* iphi = itree->GetV3();
352  //double* ipsi = itree->GetV4();
353  double* igps = itree->GetV4();
354 
355  skymap sm;
356 
357  nGWGC=0;
358  for(int n=0;n<isize;n++) {
359 
360  double dist = idistance[n];
361 
362  if (dist<DISTANCE_THR) {
363 
364  double th = 0;
365  th = acos(itheta[n]);
366  th*= 180/Pi;
367  th = -(th-90);
368 
369  double ph = 0;
370  ph = iphi[n] > 0 ? iphi[n] : 2*Pi+iphi[n];
371  ph*= 180/Pi;
372  //ph = sm.RA2phi(ph, igps[n]);
373  ph = sm.phi2RA(ph, igps[n]);
374  ph-=180;
375 /*
376  double psi = 0;
377  psi = ipsi[n];
378  psi*= 180/Pi;
379 */
380  x[nGWGC]=ph;
381  y[nGWGC]=th;
382  z[nGWGC]=dist;
383  nGWGC++;
384  }
385  }
386  cout << "nGWGC : " << nGWGC << endl;
387 
388  int size=nGWGC;
389  for (int i=0;i<360*RESOLUTION;i++) {
390  for (int j=0;j<180*RESOLUTION;j++) {
391  double ph = i/(double)RESOLUTION;
392  double th = j/(double)RESOLUTION;
393  x[size]=ph;
394  y[size]=th;
395  z[size]=DISTANCE_THR;
396  size++;
397  }
398  }
399 
400  Int_t *index = new Int_t[size];
401  TMath::Sort(size,z,index,true);
402  double T[NMAX];
403  for (int i=0;i<size;i++) T[i]=x[index[i]];
404  for (int i=0;i<size;i++) x[i]=T[i];
405  for (int i=0;i<size;i++) T[i]=y[index[i]];
406  for (int i=0;i<size;i++) y[i]=T[i];
407  for (int i=0;i<size;i++) T[i]=z[index[i]];
408  for (int i=0;i<size;i++) z[i]=T[i];
409  delete [] index;
410 
411  cout << "size : " << size << endl;
412  return size;
413 }
414 
gSM
gskymap * gSM
Definition: cwb_draw_antpat.C:67
tdist
TObjString * tdist
Definition: ConvertGWGC.C:42
dist
double dist
Definition: ConvertGWGC.C:48
sm2
skymap * sm2
Definition: TestTObjectSkymap.C:27
SKYMASK_RESOLUTION
#define SKYMASK_RESOLUTION
Definition: DrawGWGC.C:40
name
par [0] name
Definition: CWB_Plugin_HEN_SIM_Config.C:93
SGRA_NAME
#define SGRA_NAME
Definition: DrawGWGC.C:12
n
int n
Definition: cwb_net.C:28
z
wavearray< double > z
Definition: Test10.C:32
TString
TString("c")
Definition: cwb_report_skymap.C:129
out
ofstream out
Definition: cwb_merge.C:214
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<< "\ 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
Definition: BaseFunctionWDM.C:51
gskymap::DrawMarker
void DrawMarker(double phi, double theta, int marker, Size_t msize=1, Color_t tcolor=1)
Definition: gskymap.cc:742
gskymap::DrawText
void DrawText(double phi, double theta, TString text, double tsize=0.04, Color_t tcolor=1)
Definition: gskymap.cc:799
ph
TH2F * ph
Definition: ReadHistFromJobFile.C:17
skymap::getTheta
double getTheta(size_t i)
Definition: skymap.hh:224
gskymap::Draw
void Draw(int dpaletteId=1, Option_t *option="colfz")
Definition: gskymap.cc:460
size
Long_t size
Definition: cwb_condor_check.C:66
skymap::getPhiStep
double getPhiStep(size_t i)
Definition: skymap.hh:182
j
int j
Definition: cwb_net.C:28
i
i drho i
Definition: cwb_epparameters.C:88
ReadGWGCCatalogToSkymap
int ReadGWGCCatalogToSkymap(skymap &sm)
Definition: DrawGWGC.C:170
gskymap::SetGridxColor
void SetGridxColor(Color_t colorGridx=kBlack)
Definition: gskymap.hh:142
isize
int isize
Definition: cwb_condor_create_ced.C:68
ReadGWGCCatalog
int ReadGWGCCatalog(double *&x, double *&y, double *&z, int &nGWGC)
Definition: DrawGWGC.C:212
DrawGWGC
void DrawGWGC()
Definition: DrawGWGC.C:59
ReadBRST50MPC_S6
int ReadBRST50MPC_S6(double *&x, double *&y, double *&z, int &nGWGC)
Definition: DrawGWGC.C:336
GWGCCatalog
#define GWGCCatalog
Definition: DrawGWGC.C:8
skymap::getSkyIndex
size_t getSkyIndex(double th, double ph)
param: theta param: phi
Definition: skymap.cc:720
gskymap::FillData
void FillData(int size, double *phi, double *theta, double *binc)
Definition: gskymap.cc:394
ra
double ra
Definition: ConvertGWGC.C:46
NMAX
#define NMAX
Definition: DrawGWGC.C:34
gskymap::SetGalacticDiskColor
void SetGalacticDiskColor(Color_t colorGalacticDisk=kBlack)
Definition: gskymap.hh:165
RESOLUTION
#define RESOLUTION
Definition: DrawGWGC.C:6
Pi
double Pi
Definition: DrawPhaseShift.C:12
l
int l
Definition: cwb_report_cbc.C:775
gskymap::SetGalacticDisk
void SetGalacticDisk(double gpsGalacticDisk=0.0)
Definition: gskymap.hh:163
DISTANCE_THR
#define DISTANCE_THR
Definition: DrawGWGC.C:9
skymap
Definition: skymap.hh:63
gskymap::SetGridyColor
void SetGridyColor(Color_t colorGridy=kBlack)
Definition: gskymap.hh:146
ReadBRST50MPC_S6_ToSkymask
int ReadBRST50MPC_S6_ToSkymask(skymap &sm)
Definition: DrawGWGC.C:293
ifile
TFile * ifile
Definition: cwb_dump_config.C:41
ofileName
TString ofileName
Definition: MergeTrees.C:37
skymap::phi2RA
double phi2RA(double ph, double gps)
Definition: skymap.hh:212
ind
int * ind
Definition: cwb_fix_live_slag_missed.C:107
gskymap::SetTitle
void SetTitle(TString title)
Definition: gskymap.hh:152
SGRA_RA
#define SGRA_RA
Definition: DrawGWGC.C:14
COORDINATES
#define COORDINATES
Definition: DrawGWGC.C:28
skymap::getPhi
double getPhi(size_t i)
Definition: skymap.hh:164
title
char title[256]
Definition: SSeriesExample.C:1
T
double T
Definition: testWDM_4.C:11
in
ifstream in
Definition: CWB_Plugin_HEN_BKG_Config.C:26
index
wavearray< int > index
Definition: TestEulerCharacteristic.C:5
SGRA_DEC
#define SGRA_DEC
Definition: DrawGWGC.C:13
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
mask
double mask
Definition: cwb1G_parameters.C:140
tra
TObjString * tra
Definition: ConvertGWGC.C:40
itree
TTree * itree
Definition: cwb_dump_events.C:45
skymap::set
void set(size_t i, double a)
param: sky index param: value to set
Definition: skymap.hh:122
tdec
TObjString * tdec
Definition: ConvertGWGC.C:41
skymap::get
double get(size_t i)
param: sky index
Definition: skymap.cc:699
skymap::size
size_t size()
Definition: skymap.hh:136
PROJECTION
#define PROJECTION
Definition: DrawGWGC.C:29
gskymap::Print
void Print(TString pname)
Definition: gskymap.cc:1122
y
wavearray< double > y
Definition: Test10.C:31
gskymap::SetOptions
void SetOptions(TString projection="hammer", TString coordinate="Geographic", double resolution=1, bool goff=false)
Definition: gskymap.cc:84
gskymap::SetZaxisTitle
void SetZaxisTitle(TString zAxisTitle)
Definition: gskymap.hh:150
SKYMASK_FILE_NAME
#define SKYMASK_FILE_NAME
Definition: DrawGWGC.C:36
exit
exit(0)