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