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