Logo coherent WaveBurst  
Config Reference Guide 		Logo
cwb_mkhtml_chunk.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 #include "ReadChunkList.C"
20 
21 #define CHUNK_FILE_LIST "Chunk_List.txt"
22 #define CHUNK_MAX_SIZE 100
23 
24 #define MAX_CED_ENTRIES 10
25 #define MIN_CED_RHO 4.5
26 //#define MIN_CED_RHO 0.0
27 
28 #define WWW_PUBLIC "https://ldas-jobs.ligo.caltech.edu/~waveburst/reports/"
29 #define WWW_LAG_MANUAL "https://gwburst.gitlab.io/documentation/latest/html/faq.html?highlight=lag#what-are-lags-and-how-to-use-them"
30 #define WWW_SLAG_MANUAL "https://gwburst.gitlab.io/documentation/latest/html/faq.html?highlight=lag#what-are-super-lags-and-how-to-use-them"
31 
32 
33 void GetChunkRange(TString run, TString search, int ichunk, double& xstart, double& xstop);
34 void ComputeSigmaFAP(double OBSERVATIONAL_TIME, double BACKGROUND_TIME, int TRIALS_FACTOR=1);
35 void ComputeSigmaFAP(double& fap, double& sigma, double ifar, double obs_time);
36 
37 void cwb_mkhtml_chunk(TString run, TString search, TString ibin, int ichunk, bool bbh, TString dir_bkg, TString dir_frg, TString odir, TString wlabel, TString net_file_name, int lag, int slag, int nIFO) {
38 
39  cout<<"cwb_mkhtml_chunk.C starts..."<<endl;
40  //cout << odir << endl;exit(0);
41 
42  netevent W(net_file_name,nIFO);
43 
44  gSystem->Exec("date");
45 
46  // Get detector names
47  TString net="";
48  vector<TString> ifo;
49  TList* ifoList = W.fChain->GetUserInfo();
50  detectorParams dParams[NIFO_MAX];
51  for (int n=0;n<ifoList->GetSize();n++) {
52  detector* pDetector = (detector*)ifoList->At(n);
53  dParams[n] = pDetector->getDetectorParams();
54  ifo.push_back(dParams[n].name);
55  net+=TString(dParams[n].name);
56  pDetector->print();
57  }
58  if(ifo.size()==0) {
59  cout << "ifo names are not contained in the input wave root file: " << net_file_name << endl;
60  gSystem->Exit(1);
61  }
62 
63  int rho_id;
64  if(search=="BBH" || search=="IMBHB") rho_id=1;
65  else rho_id=0;
66 
67  // get CWB_USER_URL
68  char cwb_user_url[1024] = WWW_PUBLIC;
69  if(gSystem->Getenv("CWB_USER_URL")!=NULL) {
70  strcpy(cwb_user_url,TString(gSystem->Getenv("CWB_USER_URL")).Data());
71  }
72 
73  char schunk[512];sprintf(schunk,"%02d",ichunk);
74  char search_title[512];
75  if(search=="BBH" || search=="IMBHB") {
76  // check if it is a combined search
77  TString tag = dir_bkg.Contains(".J_") ? "Combined " : "";
78  sprintf(search_title,"%s %s %s Bin%s %sSearch",run.Data(),net.Data(),search.Data(),ibin.Data(),tag.Data());
79  } else {
80  sprintf(search_title,"%s %s %s Bin%s Search",run.Data(),net.Data(),search.Data(),ibin.Data());
81  }
82 
83  double xstart=0;
84  double xstop=0;
85  GetChunkRange(run, search, ichunk, xstart, xstop);
86  cout << xstart << " " << xstop << endl;
87 
88  double interval = (xstop-xstart)/(24.*3600.);
89 
90  wat::Time beg_date(xstart);
91  wat::Time end_date(xstop);
92 
93  TString sbeg_date = beg_date.GetDateString();sbeg_date.Resize(19);
94  TString send_date = end_date.GetDateString();send_date.Resize(19);
95 
96  char period[1024];
97  sprintf(period,"GPS Interval [%d,%d]. UTC Interval %s - %s. Interval duration = %.2f days.",int(xstart),int(xstop),sbeg_date.Data(),send_date.Data(),interval);
98 
99  char box_title[1024];
100  if(lag==0 && slag==0)
101  sprintf(box_title,"Open Box Result");
102  else
103  sprintf(box_title,"Fake Open Box Result - ( <td><a href=\"%s\" target=\"_blank\">LAG</a></td> = %d - <td><a href=\"%s\" target=\"_blank\">SLAG</a></td> = %d )",WWW_LAG_MANUAL,lag,WWW_SLAG_MANUAL,slag);
104 // sprintf(box_title,"Fake Open Box Result - ( LAG = %d - SLAG = %d )",lag,slag);
105 
106  char bkg_rep[1024];
107  sprintf(bkg_rep,"<td><a href=\"%s/%s/%s/index.html\" target=\"_blank\">Background Report</a></td>",cwb_user_url,wlabel.Data(),dir_bkg.Data());
108 
109  // get livetime
110 
111  char fileliv[1024];
112  sprintf(fileliv,"../../%s/data/live.txt",dir_bkg.Data());
113  cout << fileliv << endl;
114  int countlag=0;
115  double OLIVETIME = GetLiveTime(fileliv,slag,lag,countlag); // get zero livetime
116  double LIVETIME = GetLiveTime(fileliv,-1,-1,countlag); // get total livetime
117  LIVETIME-=OLIVETIME; // non zero live time
118  countlag-=1; // subtract the zero lag
119  cout.precision(14);
120  cout << OLIVETIME << " " << LIVETIME << endl;
121  //ComputeSigmaFAP(OLIVETIME, LIVETIME);
122 
123  char livetime[1024];
124  sprintf(livetime,"Livetime - Background: %.2f years, Foreground: %.2f days",LIVETIME/(24.*3600.*365.),OLIVETIME/(24.*3600.));
125 
126  // create body.html file
127 
128  ofstream out;
129  char fileout[1024];
130  sprintf(fileout,"%s/body.html", odir.Data());
131  cout << fileout << endl;
132  out.open(fileout,ios::out);
133  if (!out.good()) {cout << "Error Opening File : " << fileout << endl;exit(1);}
134 
135  out << "<html>" << endl;
136 
137 // out << "<br>" << endl;
138  out << "<div align=\"center\"><font color=\"blue\"><h1>" << search_title << " : Chunk " << schunk << "</h1></font></div>" << endl;
139  out << "<div align=\"center\"><font color=\"red\"><h4>"<< box_title <<"</h4></font></div>" << endl;
140  out << "<div align=\"center\"><font color=\"black\"><h4>"<< period <<"</h4></font></div>" << endl;
141  out << "<div align=\"center\"><font color=\"black\"><h4>"<< livetime <<"</h4></font></div>" << endl;
142  out << "<div align=\"center\"><font color=\"black\"><h4>"<< bkg_rep <<"</h4></font></div>" << endl;
143 // out << "<br>" << endl;
144 
145  out << "<hr>" << endl;
146  out << "<br>" << endl;
147  out << "<br>" << endl;
148 
149  out << "<table>" << endl;
150  out << "<tr><td width=\"50%\"><div align=\"center\">" << endl;
151  out << "<font color=\"red\"><h2>FAR vs Rank</h2></font>" << endl;
152  out << "</div><div align=\"center\"><ul><br/>" << endl;
153  out << "<a class=\"image\" title=\"FARvsRank\">" << endl;
154  out << "<img src=\"FARvsRank.png\" width=\"470\"> </a>" << endl;
155  out << "</br></ul>" << endl;
156  out << "</div><br><br></td><td width=\"50%\"><div align=\"center\">" << endl;
157  out << "<font color=\"red\"><h2>Cumulative Number vs IFAR</h2></font>" << endl;
158  out << "</div><div align=\"center\"><ul><br/>" << endl;
159  out << "<a class=\"image\" title=\"CumulativeNumberVsIFAR_bbh_plot\">" << endl;
160  if(bbh) out << "<img src=\"CumulativeNumberVsIFAR_bbh_plot.png\" width=\"470\"> </a>" << endl;
161  else out << "<img src=\"CumulativeNumberVsIFAR_nobbh_plot.png\" width=\"470\"> </a>" << endl;
162  out << "</br></ul>" << endl;
163  out << "</div><br><br></td></tr>" << endl;
164  out << "</table>" << endl;
165 
166  out << "<table>" << endl;
167  out << "<tr><td width=\"50%\"><div align=\"center\">" << endl;
168  out << "<font color=\"red\"><h2>Background: Rank vs Time</h2></font>" << endl;
169  out << "</div><div align=\"center\"><ul><br/>" << endl;
170  out << "<a class=\"image\" title=\"rho_time\">" << endl;
171  out << "<img src=\"rho_time.gif\" width=\"470\"> </a>" << endl;
172  out << "</br></ul>" << endl;
173  out << "</div><br><br></td><td width=\"50%\"><div align=\"center\">" << endl;
174  out << "<font color=\"red\"><h2>Background: Rank vs Frequency</h2></font>" << endl;
175  out << "</div><div align=\"center\"><ul><br/>" << endl;
176  out << "<a class=\"image\" title=\"rho_frequency\">" << endl;
177  out << "<img src=\"rho_frequency.gif\" width=\"470\"> </a>" << endl;
178  out << "</br></ul>" << endl;
179  out << "</div><br><br></td></tr>" << endl;
180  out << "</table>" << endl;
181 
182 
183  char cut[1024];
184  sprintf(cut,"lag[%d]==%d && slag[%d]==%d",nIFO,lag,nIFO,slag);
185 
186  char sel[1024];
187  sprintf(sel,"rho[%d]:Entry$",rho_id);
188 
189  W.fChain->SetEstimate(W.fChain->GetEntries());
190  W.fChain->Draw(sel,cut,"goff");
191  Int_t nentries = (Int_t)W.fChain->GetSelectedRows();
192  Int_t *index = new Int_t[nentries];
193  double* entry = W.fChain->GetV2();
194  TMath::Sort(nentries,W.fChain->GetV1(),index,true);
195 
196  char os[1024];
197 
198  int pp_inetcc = 0;
199 
200  out << "<head>" << endl;
201  out << "<style type=\"text/css\">" << endl;
202  out << ".datagrid tr:hover td" << endl;
203  out << "{" << endl;
204  out << " background-color:#F1F1F2;" << endl;
205  out << "}" << endl;
206  out << "</style>" << endl;
207  out << "</head>" << endl;
208  out << "<b>" << endl;
209  out << "<hr>" << endl;
210  out << "<br>" << endl;
211  out << "<font color=\"red\" style=\"font-weight:bold;\"><center><p><h2>Foreground Loudest Event List</h2><p><center></font>" << endl;
212  if(bbh) out << "<a target=\"_blank\" href=\"CumulativeNumberVsIFAR_bbh_loudest.txt\">(Events sorted by IFAR: Cumulative FAP) </a>" << endl;
213  else out << "<a target=\"_blank\" href=\"CumulativeNumberVsIFAR_nobbh_loudest.txt\">(Events sorted by IFAR: Cumulative FAP) </a>" << endl;
214  out << "<br>" << endl;
215  out << "<br>" << endl;
216  out << "</html>" << endl;
217 
218  out << "<table border=0 cellpadding=2 class=\"datagrid\">" << endl;
219  out << "<tr align=\"center\">"<< endl;
220  out << "<td>CED</td>"<< endl;
221  out << "<td>rho</td>"<< endl;
222  out << "<td>IFAR(yr)</td>"<< endl;
223 // out << "<td>FAP</td>"<< endl;
224 // out << "<td>lag</td>"<< endl;
225 // out << "<td>slag</td>"<< endl;
226  out << "<td>SNRnet</td>"<< endl;
227 // out << "<td>Correlation["<<pp_inetcc<<"]</td>"<< endl;
228  out << "<td>Correlation</td>"<< endl;
229  out << "<td>Frequency(Hz)</td>"<< endl;
230  out << "<td>Bandwidth(Hz)</td>"<< endl;
231  out << "<td>Duration(sec)</td>"<< endl;
232 // out << "<td>run</td>"<< endl;
233  for (int nn=0;nn<nIFO;nn++) out << "<td>GPS " << ifo[nn] << "</td>"<< endl;
234  for (int nn=0;nn<nIFO;nn++) out << "<td>SNR " << ifo[nn] << "</td>"<< endl;
235  out << "<td>Combined</td>"<< endl;
236  out << "</tr>"<< endl;
237 
238  out << "<td></td>" << endl;
239 
240  int segEdge=10;
241 
242  float ifar;
243  W.fChain->SetBranchAddress("ifar",&ifar);
244  int combine=0;
245  if(dir_bkg.Contains(".J_")) W.fChain->SetBranchAddress("combine",&combine);
246 
247  int ecount=0;
248  for(int i=0; i<nentries; i++) {
249 
250  if(i>=MAX_CED_ENTRIES) continue;
251 
252  W.fChain->GetEntry(entry[index[i]]);
253 
254  if(W.rho[rho_id]<=MIN_CED_RHO) continue;
255 
256  double fap;
257  double sigma;
258  ComputeSigmaFAP(fap, sigma, ifar, OLIVETIME);
259 
260  cout.precision(2);
261  cout << "EVENT " << i << "\trho = " << W.rho[rho_id] << "\tnetcc = " << W.netcc[0] << "\tifar (years) = " << ifar/(24.*3600.*365.)
262  << "\tobs (days) = " << OLIVETIME/(24.*3600.) << "\tfap = " << fap << "\tsigma = " << sigma << "\tcombine = " << combine << endl;
263  ecount++;
264  //cout << ecount << " ";cout.flush();
265  out << "<tr align=\"center\">"<< endl;
266  sprintf(os,"run==%i",W.run);
267 // TString s_ced=GetFileLabel(W.fChain,W.run,W.lag[nIFO],W.slag[nIFO],segEdge,wlabel);
268 
269  char label[1024];
270  int start = TMath::Nint(W.start[0]-W.left[0])+segEdge;
271  int stop = TMath::Nint(W.stop[0]+W.right[0])-TMath::Nint(W.start[0]-W.left[0])-2*segEdge;
272  sprintf(label,"%i_%i_%s_slag%i_lag%i_%i_job%i", start, stop, wlabel.Data(), (int)W.slag[nIFO], (int)W.lag[nIFO], 1, W.run);
273  TString s_ced=label;
274 
275  char namedir[1024];
276  sprintf(namedir,"");
277  for (int nn=0; nn<nIFO; nn++) sprintf(namedir,"%s%s",namedir,ifo[nn].Data());
278  for (int nn=0; nn<nIFO; nn++) {
279  sprintf(namedir,"%s_%.3f",namedir,W.start[nn]);
280  }
281  sprintf(os,"<td><a href=\"%s/%s/%s/ced/ced_%s/%s\" target=\"_blank\">%i</a></td>",cwb_user_url,wlabel.Data(),dir_frg.Data(),s_ced.Data(),namedir,ecount);
282 
283 // cout << os << endl;
284  out << os << endl;
285  sprintf(os,"<td>%.2f</td>",W.rho[rho_id]);
286  out << os << endl;
287  sprintf(os,"<td>%.2g</td>",ifar/(24.*3600.*365.));
288  out << os << endl;
289 // sprintf(os,"<td>%.2g</td>",fap);
290 // out << os << endl;
291 // sprintf(os,"<td>%i</td>",(int)W.lag[nIFO]);
292 // out << os << endl;
293 // sprintf(os,"<td>%i</td>",(int)W.slag[nIFO]);
294 // out << os << endl;
295  sprintf(os,"<td>%3.1f</td>",sqrt(W.likelihood));
296  out << os << endl;
297  sprintf(os,"<td>%.2f</td>",W.netcc[pp_inetcc]);
298  out << os << endl;
299  sprintf(os,"<td>%i</td>",(int)W.frequency[0]);
300  out << os << endl;
301  sprintf(os,"<td>%i</td>",(int)W.bandwidth[1]);
302  out << os << endl;
303  sprintf(os,"<td>%.3f</td>",W.duration[1]);
304  out << os << endl;
305  for (int nn=0; nn<nIFO; nn++) {
306  sprintf(os,"<td>%.2f</td>",W.time[nn]);
307  out << os << endl;
308  }
309  for (int nn=0; nn<nIFO; nn++) {
310  sprintf(os,"<td>%3.1f</td>",sqrt(W.sSNR[nn]));
311  out << os << endl;
312  }
313  sprintf(os,"<td>%d</td>",combine);
314  out << os << endl;
315 
316  out << "</tr>" << endl;
317 
318  } // End event loop
319 
320  out << "</table>" << endl;
321  out << "<p>" << endl;
322  out << endl;
323 
324  exit(0);
325 }
326 
327 void GetChunkRange(TString run, TString search, int ichunk, double& xstart, double& xstop) {
328 
329  // get CWB_CONFIG
330  char cwb_config_env[1024] = "";
331  if(gSystem->Getenv("CWB_CONFIG")!=NULL) {
332  strcpy(cwb_config_env,TString(gSystem->Getenv("CWB_CONFIG")).Data());
333  }
334 
335  char chunk_file_list[1024];
336  sprintf(chunk_file_list,"%s/%s/CHUNKS/%s/%s",cwb_config_env,run.Data(),search.Data(),CHUNK_FILE_LIST);
337  cout << chunk_file_list << endl;
338 
339  int chunk[CHUNK_MAX_SIZE];
340  double start[CHUNK_MAX_SIZE];
341  double stop[CHUNK_MAX_SIZE];
342 
343  int nChunks = ReadChunkList(chunk_file_list,chunk,start,stop);
344 
345  for(int i=0;i<nChunks;i++) {
346  if(chunk[i]==ichunk) {
347  xstart = start[i];
348  xstop = stop[i];
349  }
350  }
351 }
352 
353 void ComputeSigmaFAP(double OBSERVATIONAL_TIME, double BACKGROUND_TIME, int TRIALS_FACTOR) {
354 
355  double N = OBSERVATIONAL_TIME*(1./BACKGROUND_TIME);
356  double FAP = 1-exp(-N*TRIALS_FACTOR);
357  double Gsigma = sqrt(2)*TMath::ErfcInverse(FAP);
358 
359  double FAP_from_Gsigma = TMath::Erfc(Gsigma*1./sqrt(2)); // xcheck
360 
361  cout << endl;
362  cout << "-----------------------------------------------" << endl;
363  cout << "OBSERVATIONAL_TIME : " << OBSERVATIONAL_TIME/(24.*3600.) << " days" << endl;
364  cout << "BACKGROUND_TIME : " << BACKGROUND_TIME/(24.*3600.*365.) << " years" << endl;
365  cout << "-----------------------------------------------" << endl;
366  cout << "FAP : " << FAP << endl;
367  cout << "Gaussian sigma : " << Gsigma << endl;
368  cout << "-----------------------------------------------" << endl;
369  cout << endl;
370 
371 }
372 
373 void ComputeSigmaFAP(double& fap, double& sigma, double ifar, double obs_time) {
374 
375  double N = obs_time/ifar;
376  double FAP = 1-exp(-N);
377  double Gsigma = sqrt(2)*TMath::ErfcInverse(FAP);
378 
379  double FAP_from_Gsigma = TMath::Erfc(Gsigma*1./sqrt(2))/2; // xcheck
380 
381 /*
382  cout << endl;
383  cout << "-----------------------------------------------" << endl;
384  cout << "OBSERVATIONAL_TIME : " << obs_time/(24.*3600.) << " days" << endl;
385  cout << "BACKGROUND_TIME : " << ifar/(24.*3600.*365.) << " years" << endl;
386  cout << "-----------------------------------------------" << endl;
387  cout << "FAP : " << FAP << endl;
388  cout << "Gaussian sigma : " << Gsigma << endl;
389  cout << "-----------------------------------------------" << endl;
390  cout << endl;
391 */
392 
393  fap = FAP;
394  sigma = Gsigma;
395 }
396 
GetLiveTime
void GetLiveTime(int nwave, TChain **live, int *lag, int *slag, int *bin, TString *run, double *liveZero, double *liveTot)
Definition: Make_PP_IFAR.C:1351
segEdge
segEdge
Definition: user_parameters.C:35
nIFO
nIFO
Definition: user_parameters.C:9
strcpy
strcpy(analysis,"2G")
pp_inetcc
pp_inetcc
Definition: user_pparameters.C:17
ReadChunkList
int ReadChunkList(TString ifile, int *chunk=NULL, double *start=NULL, double *stop=NULL)
Definition: ReadChunkList.C:22
n
shift breaksw case n
Definition: cwb_clchunk.csh:75
GetChunkRange
void GetChunkRange(TString run, TString search, int ichunk, double &xstart, double &xstop)
Definition: cwb_mkhtml_chunk.C:327
MIN_CED_RHO
#define MIN_CED_RHO
Definition: cwb_mkhtml_chunk.C:25
WWW_PUBLIC
#define WWW_PUBLIC
Definition: cwb_mkhtml_chunk.C:28
WWW_SLAG_MANUAL
#define WWW_SLAG_MANUAL
Definition: cwb_mkhtml_chunk.C:30
ComputeSigmaFAP
void ComputeSigmaFAP(double OBSERVATIONAL_TIME, double BACKGROUND_TIME, int TRIALS_FACTOR=1)
Definition: cwb_mkhtml_chunk.C:353
CHUNK_FILE_LIST
#define CHUNK_FILE_LIST
Definition: cwb_mkhtml_chunk.C:21
WWW_LAG_MANUAL
#define WWW_LAG_MANUAL
Definition: cwb_mkhtml_chunk.C:29
MAX_CED_ENTRIES
#define MAX_CED_ENTRIES
Definition: cwb_mkhtml_chunk.C:24
CHUNK_MAX_SIZE
#define CHUNK_MAX_SIZE
Definition: cwb_mkhtml_chunk.C:22
cwb_mkhtml_chunk
void cwb_mkhtml_chunk(TString run, TString search, TString ibin, int ichunk, bool bbh, TString dir_bkg, TString dir_frg, TString odir, TString wlabel, TString net_file_name, int lag, int slag, int nIFO)
Definition: cwb_mkhtml_chunk.C:37
exit
exit(0)
sprintf
sprintf(tag,"wave_%s", data_label)
tag
char tag[256]
Definition: cwb_mklinks.C:46
name
par[0] name
Definition: CWB_Plugin_EBBH_Config.C:14
cWB_conf.search
string search
Definition: cWB_conf.py:63
cWB_conf.run
string run
Definition: cWB_conf.py:6
cWB_conf.label
string label
Definition: cWB_conf.py:18