Logo coherent WaveBurst  
Library Reference Guide
Logo
DrawSkyDistributionPRC.C
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 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 //
20 // Draw injected/reconstructed locations of the detected events
21 // Note : this macro is used to generate the PRC report (cwb_report merge_label prc)
22 // Author : Gabriele Vedovato
23 
24 #include "TROOT.h"
25 #include "TSystem.h"
26 #include "TFile.h"
27 #include "TTree.h"
28 #include <fstream>
29 #include <iostream>
30 #include "TGraph.h"
31 #include "TMultiGraph.h"
32 #include "TCanvas.h"
33 #include "TH1F.h"
34 #include "TMath.h"
35 #include "TH1F.h"
36 #include "TH2F.h"
37 #include <TComplex.h>
38 #include <TStyle.h>
39 #include <TRandom.h>
40 #include "TVector3.h"
41 #include "TRotation.h"
42 #include "Math/Vector3Dfwd.h"
43 #include "Math/Rotation3D.h"
44 #include "constants.hh"
45 
46 #define RESOLUTION 2
47 //#define RESOLUTION 4
48 
49 //#define COORDINATES "cWB"
50 #define COORDINATES "Geographic"
51 
52 #define PROJECTION ""
53 //#define PROJECTION "hammer"
54 //#define PROJECTION "sinusoidal"
55 
56 //#define WRITE_PLOT
57 
58 #define DRAW_EVENTS
59 
60 //#define DRAW_ANGULAR_DISTANCE
61 #define ANGULAR_DISTANCE 90
62 
63 //#define DRAW_PP
64 
65 #define DISPLAY_WORLD_MAP
66 #define WORLD_MAP_DIR "$CWB_GWAT/data/"
67 
68 using namespace CWB;
69 
70 void
72  detectorParams* dP, float T_win, int pp_inetcc, float T_cor, int pp_irho, float T_cut,
73  float T_vED, float T_pen, float T_ifar, bool binj=true, TString polarization="TENSOR", bool save_antpat=false) {
74 
75  // binj=true -> select injected directions
76  // binj=false -> select reconstructed directions
77 
78  if (!gROOT->GetClass("Polar3DVector")) gSystem->Load("libMathCore");
79 
80  using namespace ROOT::Math;
81 
82  double r2d = 180./TMath::Pi();
83  double d2r = TMath::Pi()/180.;
84 
85  double speedLight = watconstants::SpeedOfLightInVacuo();
86 
87  int nIFO = ifo.size();
88 
90  for(int n=0;n<nIFO;n++) {
91  if(ifo[n]!="") pD[n] = new detector((char*)ifo[n].Data()); // built in detector
92  else pD[n] = new detector(dP[n]); // user detector
93  }
94  for(int n=0;n<nIFO;n++) cout << n << " " << pD[n]->Name << endl;
95 
96  gnetwork* gNET = new gnetwork;
97  for(int i=0; i<nIFO; i++) gNET->add(pD[i]);
98 
99  if(polarization=="SCALAR") {
100  // setting for SGW
101  for(int n=0;n<(int)gNET->ifoListSize();n++) {
102  detector *d = gNET->getifo(n);
104  }
105  }
106 
107  // setup skymap options
108  gskymap* gSM = gNET->GetGskymap();
110 
111  TH2D* h2 = (TH2D*)gSM->GetHistogram();
112  h2->GetXaxis()->SetTitleSize(0.05);
113  h2->GetXaxis()->SetLabelSize(0.05);
114  h2->GetYaxis()->SetTitleSize(0.05);
115  h2->GetYaxis()->SetLabelSize(0.05);
116  h2->GetYaxis()->SetLabelFont(42);
117  h2->GetYaxis()->SetLabelFont(42);
118  h2->GetXaxis()->SetTitleFont(42);
119  h2->GetYaxis()->SetTitleFont(42);
120  h2->GetZaxis()->SetRangeUser(0,1.0);
121 
122  if(save_antpat) { // draw antenna patterns
123 
124  char ofileName[1024];
125  TString world_map = gSystem->ExpandPathName(WORLD_MAP_DIR);
126 
127  gSM->SetWorldMapPath(world_map.Data());
128  gSM->SetWorldMap();
129  gNET->DrawAntennaPattern(2); // |fx|/|f+|
130  gNET->DrawSitesShortLabel(kBlack);
131  gNET->DrawSites(kBlack,2.0);
132  gNET->DrawSitesArms(1000000,kWhite,3.0);
133  sprintf(ofileName,"%s/antpat_alignment.png",odir.Data());
134  cout << "Write : " << ofileName << endl;
135  gSM->Print(ofileName);
136 
137  gSM->SetWorldMapPath(world_map.Data());
138  gSM->SetWorldMap();
139  gNET->DrawAntennaPattern(3); // sqrt(|F+|^2+|Fx|^2)
140  gNET->DrawSitesShortLabel(kBlack);
141  gNET->DrawSites(kBlack,2.0);
142  gNET->DrawSitesArms(1000000,kWhite,3.0);
143  sprintf(ofileName,"%s/antpat_sensitivity.png",odir.Data());
144  cout << "Write : " << ofileName << endl;
145  gSM->Print(ofileName);
146  }
147 
148  gSM->SetWorldMap(false);
149 #ifdef DRAW_PP
150  gNET->DrawAntennaPattern(3,0,true,3);
151  cout << "order : " << gSM->getOrder() << endl;
152  cout << "size : " << gSM->size() << endl;
153  int apsize = gSM->size();
154  double* ap = new double[apsize];
155  for(int i=0;i<apsize;i++) ap[i] = gSM->get(i);
156  //for(int i=0;i<apsize;i++) cout << i << " " << gSM->get(i) << endl;
157  Int_t *iap = new Int_t[apsize];
158  TMath::Sort(apsize,ap,iap,true);
159  //for(int i=0;i<apsize;i++) cout << i << " " << ap[iap[i]] << endl;
160 #else
161  gNET->DrawAntennaPattern(3,0,true);
162 #endif
163 
164  // compute direction of D1 D2 axis -> vector v12
165  XYZVector D1(pD[0]->Rv[0],pD[0]->Rv[1],pD[0]->Rv[2]);
166  XYZVector D2(pD[1]->Rv[0],pD[1]->Rv[1],pD[1]->Rv[2]);
167  D1=D1/speedLight;
168  D2=D2/speedLight;
169  XYZVector D12 = D1-D2;
170  TVector3 vD12(D12.X(),D12.Y(),D12.Z());
171 
172  double th12 = r2d*vD12.Theta();
173  double ph12 = r2d*vD12.Phi();
174  cout << "coordinates D12 " << ph12 << " " << th12 << endl;
175  Polar3DVector v12(1, d2r*th12, d2r*ph12);
176 
177  // open reconstructed event file
178 #ifdef DRAW_EVENTS
179 
180  char wave_file_name[1024];
181  sprintf(wave_file_name,"merge/wave_%s.%s.root",data_label.Data(),merge_label.Data());
182 
183  TFile* ifile = TFile::Open(wave_file_name);
184  if(ifile==NULL) {cout<<"Error opening file : "<<wave_file_name<<endl;exit(1);}
185  TTree* itree = (TTree *) gROOT->FindObject("waveburst");
186  if(itree==NULL) {cout<<"Error opening tree : "<<"waveburst"<<endl;exit(1);}
187 
188  char cut[1024];
189  char tmp[1024];
190  sprintf(cut,"abs(time[0]-time[%d])<%f && netcc[%d]>%f && rho[%d]>%f",
191  nIFO,T_win,pp_inetcc,T_cor,pp_irho,T_cut);
192  if(T_vED>0) {strcpy(tmp,cut);sprintf(cut,"%s && neted[0]/ecor<%f",tmp,T_vED);}
193  if(T_pen>0) {strcpy(tmp,cut);sprintf(cut,"%s && penalty>%f",tmp,T_pen);}
194  if(T_ifar>0) {strcpy(tmp,cut);sprintf(cut,"%s && ifar>(24.*3600.*365.)*%f",tmp,T_ifar);}
195 
196  itree->Draw("theta[0]:phi[0]:theta[1]:phi[1]",cut,"goff");
197  int isize=itree->GetSelectedRows();
198  cout << "isize : " << isize << endl;
199 
200  double* th0 = itree->GetV1();
201  double* ph0 = itree->GetV2();
202  double* th1 = itree->GetV3();
203  double* ph1 = itree->GetV4();
204 
205 #ifdef DRAW_ANGULAR_DISTANCE
206  TH1F* h1 = new TH1F("hist","hist",100,0,180);
207 #endif
208 
209  // create gskymap to store events
210  gskymap* gEVT = new gskymap((int)3);
212  *gEVT = 0;
213 
214  // draw events
215  double markerSize = isize>10000 ? 0.3 : 0.4;
216  double ph,th;
217  for (int i=0;i<isize;i++) {
218 
219  // compute distance (dOmega) between injected and reconstructed directions
220  Polar3DVector v0(1, d2r*th0[i], d2r*ph0[i]);
221  Polar3DVector v1(1, d2r*th1[i], d2r*ph1[i]);
222  //cout << th0[i] << " " << ph0[i] << " " << th1[i] << " " << ph1[i] << endl;
223  double dot = v0.Dot(v1);
224  double dOmega = r2d*TMath::ACos(dot);
225 
226  // compute distance (dOmega12) between injected and D1 D2 axis
227  double dot12 = v12.Dot(v1);
228  double dOmega12 = r2d*TMath::ACos(dot12);
229 
230  // discart events outside the ring (width=5degrees) at distance ANGULAR_DISTANCE
231 // if(fabs(dOmega12-ANGULAR_DISTANCE)>5) continue;
232 
233 #ifdef DRAW_ANGULAR_DISTANCE
234  h1->Fill(dOmega);
235 #endif
236 
237  if(binj) {ph=ph1[i]; th=th1[i];}
238  else {ph=ph0[i]; th=th0[i];}
239 
240  if(COORDINATES=="cWB") {
241  gSM->DrawMarker(ph, th, 20, markerSize, kBlack); // cWB
242  }
243  if(COORDINATES=="Geographic") {
244  double phi,theta;
245  CwbToGeographic(ph,th,phi,theta);
246  gSM->DrawMarker(phi,theta, 20, markerSize, kBlack); // Geographic
247 
248  int index = gEVT->getSkyIndex(th,ph);
249  gEVT->set(index,gEVT->get(index)+1);
250  }
251  }
252 #endif
253 
254  // draw circle at distance -ANGULAR_DISTANCE from D1 D2 axis
255  double phi=ph12;
256  double theta=th12-ANGULAR_DISTANCE;
257  if(COORDINATES=="Geographic") {
258  CwbToGeographic(phi,theta,phi,theta);
259  }
260  gNET->DrawCircles(phi,theta,(Color_t)kWhite,1,1,true);
261 
262  gSM->GetCanvas()->Update();
263 
264  // draw angular distance
265 #ifdef DRAW_ANGULAR_DISTANCE
266  gStyle->SetLineColor(kBlack);
267  h1->Draw("HIST");
268 #endif
269 
270  // write plot
271  char ofname[1024];
272  if(binj) {
273  sprintf(ofname,"%s/inj_detected_antpat.png",odir.Data());
274  } else {
275  sprintf(ofname,"%s/rec_detected_antpat.png",odir.Data());
276  }
277  cout << "Write : " << ofname << endl;
278  gSM->Print(ofname);
279  //gSystem->Exit(0);
280 
281 #ifdef DRAW_PP
282  TCanvas* canvas = new TCanvas;
283  canvas->SetGridx();
284  canvas->SetGridy();
285 
286  double* x = new double[apsize];
287  double* y = new double[apsize];
288  for(int i=0;i<apsize;i++) {
289  //x[i]=ap[iap[i]];
290  x[i]=pow(ap[iap[i]],4);
291  y[i]=gEVT->get(iap[i]);
292 //x[i]=pow(ap[i],4);
293 //y[i]=gEVT->get(i);
294  }
295  for(int i=1;i<apsize;i++) {
296  x[i]+=x[i-1];
297  y[i]+=y[i-1];
298  }
299  for(int i=0;i<apsize;i++) {
300  x[i]/=x[apsize-1];
301  y[i]/=y[apsize-1];
302  }
303  TGraph* gr = new TGraph(apsize,x,y);
304  gr->Draw("ALP");
305  gr->SetLineColor(kRed);
306  gr->SetMarkerColor(kRed);
307 
308  TH1F* hist = gr->GetHistogram();
309  hist->SetStats(kFALSE);
310  char title[256];
311  sprintf(title,"%s - pp plot",TITLE);
312  hist->SetTitle(title);
313  //hist->GetXaxis()->SetTitle("(F+^{2}+Fx^{2}) Probability");
314  hist->GetXaxis()->SetTitle("(F+^{2}+Fx^{2})^{2} Probability");
315  hist->GetYaxis()->SetTitle("Frequentist Probability");
316  hist->GetXaxis()->SetRangeUser(0,1);
317  hist->GetYaxis()->SetRangeUser(0,1);
318 
319 // for(int i=0;i<apsize;i++) cout << i << " " << ap[iap[i]] << " " << gEVT->get(iap[i]) << endl;
320 // gEVT->Draw();
321 #endif
322 
323 }
324 
detector * getifo(size_t n)
param: detector index
Definition: network.hh:436
TH2D * GetHistogram()
Definition: gskymap.hh:138
gskymap * gSM
double T_ifar
double T_pen
gnetwork * gNET
size_t add(detector *)
param: detector structure return number of detectors in the network
Definition: network.cc:2559
void DrawSitesArms(double mlength=600000., Color_t lcolor=kBlack, Size_t lwidth=1.0, Style_t lstyle=1)
Definition: gnetwork.cc:314
Definition: ced.hh:42
void DrawAntennaPattern(int polarization=-1, int dpaletteId=0, bool btitle=true, int order=6)
Definition: gnetwork.cc:674
void DrawCircles(double phi, double theta, double gps, Color_t lcolor=kBlack, Width_t lwidth=1, Style_t lstyle=1, bool labels=false)
Definition: gnetwork.cc:1308
int n
Definition: cwb_net.C:28
TString("c")
#define TITLE
Definition: TestSiteJ1.C:44
void set(size_t i, double a)
Definition: gskymap.hh:128
double T_cor
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
char odir[1024]
float theta
void DrawSkyDistributionPRC(TString data_label, TString odir, TString merge_label, vector< TString > ifo, detectorParams *dP, float T_win, int pp_inetcc, float T_cor, int pp_irho, float T_cut, float T_vED, float T_pen, float T_ifar, bool binj=true, TString polarization="TENSOR", bool save_antpat=false)
void DrawMarker(double phi, double theta, int marker, Size_t msize=1, Color_t tcolor=1)
Definition: gskymap.cc:742
TH2F * ph
return wmap canvas
#define ANGULAR_DISTANCE
i pp_inetcc
int polarization
void setPolarization(POLARIZATION polarization=TENSOR)
Definition: detector.hh:306
i drho i
int isize
char ifo[NIFO_MAX][8]
size_t ifoListSize()
Definition: network.hh:431
#define PROJECTION
#define nIFO
size_t getSkyIndex(double th, double ph)
param: theta param: phi
Definition: skymap.cc:720
TString world_map
Definition: DrawGNET.C:16
char data_label[512]
Definition: test_config1.C:160
#define WORLD_MAP_DIR
void DrawSitesShortLabel(Color_t tcolor=kBlack, Size_t tsize=0.052, Font_t tfont=32)
Definition: gnetwork.cc:426
float phi
TGraph * gr
detector D1
i() int(T_cor *100))
double Pi
int getOrder()
Definition: skymap.hh:314
const int NIFO_MAX
Definition: wat.hh:22
double * tmp
Definition: testWDM_5.C:31
void DrawSites(Color_t mcolor=kBlack, Size_t msize=2.0, Style_t mstyle=20)
Definition: gnetwork.cc:257
char cut[512]
TCanvas * GetCanvas()
Definition: gskymap.hh:137
void SetWorldMapPath(TString worldMapPath)
Definition: gskymap.hh:156
TString merge_label
#define RESOLUTION
TFile * ifile
TString ofileName
Definition: MergeTrees.C:37
gskymap * GetGskymap()
Definition: gnetwork.hh:44
char title[256]
Definition: SSeriesExample.C:1
void SetWorldMap(bool drawWorldMap=true)
Definition: gskymap.hh:154
wavearray< int > index
double T_win
void CwbToGeographic(double ilongitude, double ilatitude, double &olongitude, double &olatitude)
Definition: skycoord.hh:414
strcpy(RunLabel, RUN_LABEL)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
double T_cut
TTree * itree
double get(size_t i)
param: sky index
Definition: skymap.cc:699
double T_vED
size_t size()
Definition: skymap.hh:136
void Print(TString pname)
Definition: gskymap.cc:1122
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
i drho pp_irho
#define COORDINATES
detector ** pD
exit(0)
double SpeedOfLightInVacuo()
Definition: constants.hh:114