Logo coherent WaveBurst  
Library Reference Guide 		Logo
CWB_Plugin_Skymap.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 #define XIFO 4
20 
21 #pragma GCC system_header
22 
23 #include "cwb.hh"
24 #include "config.hh"
25 #include "network.hh"
26 #include "wavearray.hh"
27 #include "TString.h"
28 #include "TObjArray.h"
29 #include "TObjString.h"
30 #include "TRandom.h"
31 #include "TComplex.h"
32 #include "TMath.h"
33 #include "mdc.hh"
34 #include <vector>
35 
36 
37 void
38 CWB_Plugin(TFile* jfile, CWB::config* cfg, network* net, WSeries<double>* x, TString ifo, int type) {
39 //!EVENT_REPORT
40 // This plugin save skymap to output root file
41 
42  cout << endl;
43  cout << "-----> CWB_Plugin_Skymap.C" << endl;
44  cout << "ifo " << ifo.Data() << endl;
45  cout << "type " << type << endl;
46  cout << endl;
47 
48  if(type==CWB_PLUGIN_OLIKELIHOOD) {
49 
50  int i,j,k,n,m,l;
51  int nIFO = net->ifoListSize();
52  int K = net->nLag;
53  int M = net->mdc__ID.size();
54  int ID;
55  char search = net->tYPe;
56 
57  wavearray<double> id;
58 
59  bool batch = gROOT->IsBatch();
60  gROOT->SetBatch(true);
61 
62  watplot SMS(const_cast<char*>("sms"),200,20,800,400);
63 
64  char ifostr[20] = "";
65  char strtime[1024] = "";
66  char fname[1024];
67 
68  int minTimeDet=nIFO;
69  double minTime=1e40;
70  double eventTime[NIFO];
71  double lagTime[NIFO];
72  int ifoid[NIFO],sortid[NIFO];
73  double factor = 1;
74 
75 
76  //Fill in all skymaps
77  double old_cc = net->netCC;
78  double old_rho = net->netRHO;
79  net->netCC = -1;
80  net->netRHO = 0;
81 
82  // search output root file in the system list
83  TFile* froot = NULL;
84  TList *files = (TList*)gROOT->GetListOfFiles();
85  if (files) {
86  TIter next(files);
87  TSystemFile *file;
88  TString fname;
89  bool check=false;
90  while ((file=(TSystemFile*)next())) {
91  fname = file->GetName();
92  // set output root file as the current file
93  if(fname.Contains("wave_")) {froot=(TFile*)file;froot->cd();}
94  }
95  if(!froot) {
96  cout << "CWB_Plugin_skymap.C : Error - output root file not found" << endl;
97  exit(1);
98  }
99  } else {
100  cout << "CWB_Plugin_skymap.C : Error - output root file not found" << endl;
101  exit(1);
102  }
103 
104  TDirectory* cdskymap = (TDirectory*)froot->Get("skymap");
105  if(cdskymap==NULL) cdskymap = froot->mkdir("skymap");
106 
107  netevent* EVT = new netevent(nIFO);
108 
109  int rate = int(2*net->getifo(0)->TFmap.resolution(0)+0.5);
110 
111  for(k=0; k<K; k++) { // loop over the lags
112 
113  id = net->getwc(k)->get(const_cast<char*>("ID"), 0, 'L', rate);
114 
115  for(j=0; j<(int)id.size(); j++) { // loop over cluster index
116 
117  ID = size_t(id.data[j]+0.5); // cluster id
118 
119  EVT->output(NULL,net,factor,ID,k); // get reconstructed parameters
120  if(EVT->rho[1] < cfg->cedRHO) continue; // skip events under rho threshold
121 
122  int masterDet=0;
123  int lagMin=2147483647;
124  for(n=0; n<nIFO; n++) if(EVT->lag[n]<lagMin) {lagMin=int(EVT->lag[n]);masterDet=n;}
125 
126  net->likelihood(search, net->acor, ID, k); // exec likelihood search
127 
128  // get event network time
129  double gps_start = EVT->time[masterDet]-EVT->duration[masterDet];
130  double gps_stop = EVT->time[masterDet]+EVT->duration[masterDet];
131 
132  // create a unique label
133  for(n=0; n<nIFO; n++) sprintf(strtime, "%s_%.3f",strtime,EVT->start[n]);
134 
135  //likelihood skymaps
136  TMarker inj(EVT->phi[1],90.-EVT->theta[1], 29); // injected pos (white star)
137  TMarker rec(EVT->phi[0],90.-EVT->theta[0], 29); // recostructed pos (black star)
138  TMarker det(EVT->phi[3],90.-EVT->theta[3], 20); // detected pos (black dot )
139  inj.SetMarkerSize(2.0); inj.SetMarkerColor(kWhite);
140  rec.SetMarkerSize(2.0); rec.SetMarkerColor(kBlack);
141  det.SetMarkerSize(1.0); det.SetMarkerColor(kBlack);
142 
143  char evtdir[64];
144  sprintf(evtdir, "run_%d_id_%d",net->nRun,ID);
145  TDirectory* cdevtdir = (TDirectory*)cdskymap->Get(evtdir);
146  if(cdevtdir==NULL) cdevtdir = cdskymap->mkdir(evtdir);
147  cdevtdir->cd();
148 
149  SMS.canvas->cd();
150  sprintf(fname, "sensitivity_plus");
151  SMS.plot(net->nSensitivity, 0);
152  rec.Draw(); det.Draw(); if(M) inj.Draw();
153  SMS.canvas->Write(fname);
154  SMS.clear();
155  sprintf(fname, "sensitivity_cross");
156  SMS.plot(net->nAlignment, 0);
157  rec.Draw(); det.Draw(); if(M) inj.Draw();
158  SMS.canvas->Write(fname);
159  SMS.clear();
160  sprintf(fname, "skystat");
161  SMS.plot(net->nSkyStat, 0);
162  rec.Draw(); det.Draw(); if(M) inj.Draw();
163  SMS.canvas->Write(fname);
164  SMS.clear();
165  sprintf(fname, "likelihood");
166  SMS.plot(net->nLikelihood, 0);
167  rec.Draw(); det.Draw(); if(M) inj.Draw();
168  SMS.canvas->Write(fname);
169  SMS.clear();
170  sprintf(fname, "null_energy");
171  SMS.plot(net->nNullEnergy, 0);
172  rec.Draw(); det.Draw(); if(M) inj.Draw();
173  SMS.canvas->Write(fname);
174  SMS.clear();
175  sprintf(fname, "corr_energy");
176  SMS.plot(net->nCorrEnergy, 0);
177  rec.Draw(); det.Draw(); if(M) inj.Draw();
178  SMS.canvas->Write(fname);
179  SMS.clear();
180  sprintf(fname, "penalty");
181  SMS.plot(net->nPenalty, 0);
182  rec.Draw(); det.Draw(); if(M) inj.Draw();
183  SMS.canvas->Write(fname);
184  SMS.clear();
185  sprintf(fname, "disbalance");
186  SMS.plot(net->nDisbalance, 0);
187  rec.Draw(); det.Draw(); if(M) inj.Draw();
188  SMS.canvas->Write(fname);
189  SMS.clear();
190  sprintf(fname, "correlation");
191  SMS.plot(net->nCorrelation, 0);
192  rec.Draw(); det.Draw(); if(M) inj.Draw();
193  SMS.canvas->Write(fname);
194  SMS.clear();
195  sprintf(fname, "netindex");
196  SMS.plot(net->nNetIndex, 0);
197  rec.Draw(); det.Draw(); if(M) inj.Draw();
198  SMS.canvas->Write(fname);
199  SMS.clear();
200  sprintf(fname, "ellipticity");
201  SMS.plot(net->nEllipticity, 0);
202  rec.Draw(); det.Draw(); if(M) inj.Draw();
203  SMS.canvas->Write(fname);
204  SMS.clear();
205  sprintf(fname, "polarisation");
206  SMS.plot(net->nPolarisation, 0);
207  rec.Draw(); det.Draw(); if(M) inj.Draw();
208  SMS.canvas->Write(fname);
209  SMS.clear();
210  sprintf(fname, "probability");
211  SMS.plot(net->nProbability, 0);
212  rec.Draw(); det.Draw(); if(M) inj.Draw();
213  SMS.canvas->Write(fname);
214  SMS.clear();
215 
216  } // End loop on found events
217  } // End loop on lags
218 
219  // restore NET thresholds
220  net->netCC = old_cc;
221  net->netRHO = old_rho;
222  }
223 
224  return;
225 }
226 
network::getifo
detector * getifo(size_t n)
param: detector index
Definition: network.hh:436
cfg
CWB::config * cfg
Definition: TestCWB_Plugin.C:23
network::netCC
double netCC
Definition: network.hh:596
NIFO
#define NIFO
Definition: wat.hh:74
network::nLag
size_t nLag
Definition: network.hh:573
netevent::rho
Float_t * rho
biased null statistics
Definition: netevent.hh:112
M
double M
Definition: DrawEBHH.C:13
netevent::start
Double_t * start
cluster duration = stopW-startW
Definition: netevent.hh:99
netevent::duration
Float_t * duration
max cluster time relative to segment start
Definition: netevent.hh:98
CWB::config::cedRHO
double cedRHO
Definition: config.hh:298
factor
float factor
Definition: cwb_report_cbc.C:777
CWB_Plugin
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *net, WSeries< double > *x, TString ifo, int type)
COHERENCE.
Definition: CWB_Plugin_Skymap.C:38
n
int n
Definition: cwb_net.C:28
network::nSkyStat
skymap nSkyStat
Definition: network.hh:628
TString
TString("c")
Definition: cwb_report_skymap.C:129
network::nRun
size_t nRun
Definition: network.hh:572
ID
int ID
Definition: TestMDC.C:70
network::nPenalty
skymap nPenalty
Definition: network.hh:624
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
netcluster::get
wavearray< double > get(char *name, size_t index=0, char atype='R', int type=1, bool=true)
param: string with parameter name param: index in the amplitude array, which define detector param: c...
Definition: netcluster.cc:2207
cWB_conf.search
string search
Definition: cWB_conf.py:110
ifostr
char ifostr[64]
Definition: cwb_report_skymap.C:142
SMS
#define SMS
Definition: FrDisplay.cc:126
network::nPolarisation
skymap nPolarisation
Definition: network.hh:630
size
Long_t size
Definition: cwb_condor_check.C:66
m
int m
Definition: cwb_net.C:28
network::mdc__ID
std::vector< size_t > mdc__ID
Definition: network.hh:615
j
int j
Definition: cwb_net.C:28
i
i drho i
Definition: cwb_epparameters.C:88
watplot::plot
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
network::nCorrEnergy
skymap nCorrEnergy
Definition: network.hh:625
ifo
char ifo[NIFO_MAX][8]
Definition: cwb1G_parameters.C:17
net
network ** net
NOISE_MDC_SIMULATION.
Definition: CWB_Plugin_HEN_BKG_Config.C:4
network::ifoListSize
size_t ifoListSize()
Definition: network.hh:431
network::nLikelihood
skymap nLikelihood
Definition: network.hh:622
watplot::clear
void clear()
Definition: watplot.hh:58
nIFO
#define nIFO
Definition: DrawSensitivitiesS5.C:11
watplot::canvas
TCanvas * canvas
Definition: watplot.hh:192
jfile
jfile
Definition: cwb_job_obj.C:43
int
i() int(T_cor *100))
l
int l
Definition: cwb_report_cbc.C:775
next
TIter next(twave->GetListOfBranches())
netevent::lag
Float_t * lag
time between consecutive events
Definition: netevent.hh:84
fname
char fname[1024]
Definition: cwb_report_prod_1.C:128
k
int k
Definition: cwb_report_prod_2.C:168
network::acor
double acor
Definition: network.hh:585
netevent::output
void output(TTree *=NULL, network *=NULL, double=0., size_t=0, int=-1)
Definition: netevent.cc:1193
network::likelihood
long likelihood(char='E', double=sqrt(2.), int=0, size_t=0, int=-1, bool=false)
Definition: network.cc:4446
froot
TFile * froot
Definition: TestReadWriteDetectorObject.C:11
detector::TFmap
WSeries< double > TFmap
Definition: detector.hh:354
netevent::phi
Float_t * phi
sqrt(h+*h+ + hx*hx)
Definition: netevent.hh:87
network::nAlignment
skymap nAlignment
Definition: network.hh:620
EVT
netevent EVT(itree, nifo)
netevent::theta
Float_t * theta
[0]-reconstructed, [1]-injected phi angle, [2]-RA
Definition: netevent.hh:88
netevent::time
Double_t * time
beam pattern coefficients for hx
Definition: netevent.hh:94
network::nDisbalance
skymap nDisbalance
Definition: network.hh:627
WSeries::resolution
double resolution(int=0)
Definition: wseries.hh:155
network::getwc
netcluster * getwc(size_t n)
param: delay index
Definition: network.hh:439
network::nNetIndex
skymap nNetIndex
Definition: network.hh:626
sprintf
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
network::nSensitivity
skymap nSensitivity
list of wdm tranformations
Definition: network.hh:619
K
int K
Definition: CWB_Plugin_Sim4_BRST_Config.C:87
id
Long_t id
Definition: cwb_condor_check.C:66
network::nNullEnergy
skymap nNullEnergy
Definition: network.hh:623
network::nCorrelation
skymap nCorrelation
Definition: network.hh:621
check
int check
Definition: CWB_configPlugin_Test.C:5
network::tYPe
char tYPe
Definition: network.hh:588
network::netRHO
double netRHO
Definition: network.hh:597
network::nProbability
skymap nProbability
Definition: network.hh:631
exit
exit(0)
network::nEllipticity
skymap nEllipticity
Definition: network.hh:629