4 cout<<
"cwb_report_prod_2.C starts..."<<endl;
21 gStyle->SetTitleOffset(1.4,
"X");
23 gStyle->SetTitleOffset(2.0,
"Y");
24 gStyle->SetLabelOffset(0.010,
"X");
25 gStyle->SetLabelOffset(0.010,
"Y");
26 gStyle->SetLabelFont(42,
"X");
27 gStyle->SetLabelFont(42,
"Y");
28 gStyle->SetTitleFont(42,
"X");
29 gStyle->SetTitleFont(42,
"Y");
30 gStyle->SetTitleSize(0.05,
"X");
31 gStyle->SetTitleSize(0.04,
"Y");
32 gStyle->SetLabelSize(0.06,
"X");
33 gStyle->SetLabelSize(0.06,
"Y");
35 gStyle->SetTitleH(0.060);
36 gStyle->SetTitleW(0.95);
37 gStyle->SetTitleY(0.99);
38 gStyle->SetTitleFont(12,
"D");
39 gStyle->SetTitleColor(kBlue,
"D");
40 gStyle->SetTextFont(12);
41 gStyle->SetTitleFillColor(kWhite);
42 gStyle->SetLineColor(kWhite);
43 gStyle->SetNumberContours(256);
44 gStyle->SetCanvasColor(kWhite);
45 gStyle->SetStatBorderSize(1);
46 gStyle->SetOptStat(kFALSE);
49 gStyle->SetFrameBorderMode(0);
69 TCanvas *
c1 =
new TCanvas(
"c",
"C",0,0,600,600);
76 c1->SetRightMargin(0.1517039);
77 c1->SetTopMargin(0.0772727);
78 c1->SetBottomMargin(0.103939);
79 gStyle->SetPalette(1,0);
81 TChain
wave(
"waveburst");
82 TChain
live(
"liveTime");
100 cout <<
"cwb_report_prod_2.C Error : Leaf " << veto_name
101 <<
" is not present in tree" << endl;
109 TIter
next(
wave.GetListOfBranches());
110 while ((branch=(TBranch*)
next())) {
111 if (
TString(
"slag").CompareTo(branch->GetName())==0) slagFound=
true;
115 cout <<
"cwb_report_prod_2.C : Error - wave tree do not contains slag leaf : "
133 if(bifo) XDQF[nXDQF++]=
VDQF[
j];
142 nameVETO[
j]=veto_name;
144 W.fChain->SetBranchAddress(nameVETO[
j].Data(),&VETO[
j]);
147 gSystem->Exec(
"date");
175 wave.SetEstimate(nsize);
176 wave.Draw(TString::Format(
"rho[%d]",
pp_irho).Data(),
"",
"goff");
184 double far_drho = double(far_rho_max-far_rho_min)/double(far_rho_bin);
208 cout <<
"Start CWB::Toolbox::getLiveTime : be patient, it takes a while ..." << endl;
213 for(j=0; j<Wlag[
nIFO].
size(); j++) {
214 lagslag[std::make_pair(
int(Wslag[nIFO].data[j]),
int(Wlag[nIFO].data[j]))] =
j;
217 Rlag =
Tlag; Rlag = 0.;
218 Elag =
Tlag; Elag = 0.;
219 Olag =
Tlag; Olag = 0.;
224 cout <<
"cwb_report_prod_2.C Error : Error opening " << fname << endl;
232 TChain live1(
"liveTime");
238 fprintf(flive,
"slag=%i\t lag=%i\t ",0,0);
239 for (
int jj=0; jj<
nIFO; jj++)
fprintf(flive,
"%s:slag=%5.2f\tlag=%5.2f\t",IFO[jj],0.,0.);
240 fprintf(flive,
"live=%9.2f\n",liveZero);
241 printf(
"live time: zero lags = %10.1f \n",liveZero);
244 fprintf(flive,
"slag=%i\t lag=%i\t ",(
int)Wslag[nIFO].data[j],(
int)Wlag[nIFO].data[j]);
245 for (
int jj=0; jj<
nIFO; jj++) {
246 fprintf(flive,
"%s:slag=%5.2f\tlag=%5.2f\t",IFO[jj],Wslag[jj].data[j],Wlag[jj].data[j]);
250 fprintf(flive,
"nonzero lags live=%10.1f \n",liveTot);
252 printf(
"total live time: non-zero lags = %10.1f \n",liveTot);
253 cout <<
"End CWB::Toolbox::getLiveTime : ";cout.flush();gSystem->Exec(
"/bin/date");
259 for(
int i=0;i<
nlag;i++)
if(Wlag[nIFO][i]<min_lag) min_lag=Wlag[
nIFO][
i];
260 for(
int i=0;i<nlag;i++) if(Wlag[nIFO][i]>
max_lag) max_lag=Wlag[nIFO][i];
261 int*
iy_lag =
new int[max_lag+1];
262 for(
int i=0;i<=
max_lag;i++) iy_lag[i]=-1;
263 for(
int i=0;i<
nlag;i++) iy_lag[
int(Wlag[nIFO][i])-
min_lag]=
i;
264 double** y_lag =
new double*[
nlag];
265 double** y_lag_veto =
new double*[
nlag];
266 for(
int i=0;i<
nlag;i++) {
272 for (
int i=0;i<
nlag;i++)
for (
int j=0; j<
pp_rho_bin; j++) {y_lag[
i][
j] = 0.;y_lag_veto[
i][
j] = 0.;}
276 TH1F*
xCor =
new TH1F(
"xCor",
"",100,0.4,1.);
277 xCor->SetTitleOffset(1.3,
"Y");
279 xCor->GetXaxis()->SetTitle(
"network correlation");
280 xCor->GetYaxis()->SetTitle(
"events");
282 TH1F*
dens =
new TH1F(
"dens",
"",100,0.,3);
283 dens->SetTitleOffset(1.3,
"Y");
285 dens->GetXaxis()->SetTitle(
"-log10(cluster density)");
286 dens->GetYaxis()->SetTitle(
"events");
288 TH1F*
hpf =
new TH1F(
"hpf",
"",100,0.,1.);
289 hpf->SetTitleOffset(1.3,
"Y");
291 hpf->GetXaxis()->SetTitle(
"penalty factor");
292 hpf->GetYaxis()->SetTitle(
"events");
294 TH1F*
hvED =
new TH1F(
"hvED",
"",200,0,2.);
295 hvED->SetTitleOffset(1.3,
"Y");
297 hvED->GetXaxis()->SetTitle(
"hvED consistency");
298 hvED->GetYaxis()->SetTitle(
"events");
300 TH1F*
ecor =
new TH1F(
"ecor",
"",200,0,10.);
301 ecor->SetTitleOffset(1.3,
"Y");
303 ecor->GetXaxis()->SetTitle(
"correlated amplitude");
304 ecor->GetYaxis()->SetTitle(
"events");
306 TH1F*
ECOR =
new TH1F(
"ECOR",
"",200,0,10.);
307 ECOR->SetTitleOffset(1.3,
"Y");
309 ECOR->GetXaxis()->SetTitle(
"#rho(ECOR))");
310 ECOR->GetYaxis()->SetTitle(
"events");
312 TH1F*
Like =
new TH1F(
"Like",
"",100,1,4.);
313 Like->SetTitleOffset(1.3,
"Y");
315 Like->GetXaxis()->SetTitle(
"log10(likelihood)");
316 Like->GetYaxis()->SetTitle(
"events");
318 TH1F*
Mchirp =
new TH1F(
"Mchirp",
"",100,0.,40.);
319 Mchirp->SetTitleOffset(1.3,
"Y");
320 Mchirp->SetTitle(
"After pp cuts");
321 Mchirp->GetXaxis()->SetTitle(
"reconstructed chirp mass(M_{#odot})");
322 Mchirp->GetYaxis()->SetTitle(
"events");
324 TH1F*
McErr =
new TH1F(
"McErr",
"",100,0.,1.);
325 McErr->SetTitleOffset(1.3,
"Y");
326 McErr->SetTitle(
"After pp cuts");
327 McErr->GetXaxis()->SetTitle(
"pixels used in reconstruction/total");
328 McErr->GetYaxis()->SetTitle(
"events");
330 TH1F*
cutF =
new TH1F(
"cutF",
"",24,60,2048.);
331 cutF->SetTitleOffset(1.2,
"Y");
332 cutF->GetXaxis()->SetTitle(
"frequency, Hz");
334 cutF->SetStats(kFALSE);
335 cutF->SetFillColor(2);
341 rhocc->SetTitle(
"0 < cc < 1");
342 rhocc->SetTitleOffset(1.3,
"Y");
343 rhocc->GetXaxis()->SetTitle(
"network correlation");
344 rhocc->GetYaxis()->SetTitle(
"#rho");
345 rhocc->SetStats(kFALSE);
346 rhocc->SetMarkerStyle(20);
347 rhocc->SetMarkerSize(0.5);
348 rhocc->SetMarkerColor(1);
350 TH2F*
ECOR_cc =
new TH2F(
"ECOR_cc",
"",100,0.,1.,500,0,15.);
351 ECOR_cc->SetTitleOffset(1.3,
"Y");
352 ECOR_cc->GetXaxis()->SetTitle(
"network correlation");
353 ECOR_cc->GetYaxis()->SetTitle(
"#rho(ECOR)");
354 ECOR_cc->SetStats(kFALSE);
355 ECOR_cc->SetMarkerStyle(20);
356 ECOR_cc->SetMarkerSize(0.5);
357 ECOR_cc->SetMarkerColor(1);
359 TH2F*
ecor_cc =
new TH2F(
"ecor_cc",
"",100,0.,1.,500,0,15.);
360 ecor_cc->SetTitleOffset(1.3,
"Y");
361 ecor_cc->GetXaxis()->SetTitle(
"network correlation");
362 ecor_cc->GetYaxis()->SetTitle(
"#rho(ecor)");
363 ecor_cc->SetStats(kFALSE);
364 ecor_cc->SetMarkerStyle(20);
365 ecor_cc->SetMarkerSize(0.5);
366 ecor_cc->SetMarkerColor(1);
369 rho_subnet->SetTitle(
"0 < subnet < 1");
370 rho_subnet->SetTitleOffset(1.3,
"Y");
371 rho_subnet->GetXaxis()->SetTitle(
"subnet");
372 rho_subnet->GetYaxis()->SetTitle(
"#rho");
373 rho_subnet->SetMarkerStyle(20);
374 rho_subnet->SetMarkerColor(1);
375 rho_subnet->SetMarkerSize(0.5);
376 rho_subnet->SetStats(kFALSE);
379 rho_vED->SetTitle(
"0 < vED < 1");
380 rho_vED->SetTitleOffset(1.3,
"Y");
381 rho_vED->GetXaxis()->SetTitle(
"vED");
382 rho_vED->GetYaxis()->SetTitle(
"#rho");
383 rho_vED->SetMarkerStyle(20);
384 rho_vED->SetMarkerColor(1);
385 rho_vED->SetMarkerSize(0.5);
386 rho_vED->SetStats(kFALSE);
391 rho_pf->SetTitle(
"0 < penalty < 1");
392 rho_pf->GetXaxis()->SetTitle(
"penalty");
395 rho_pf->SetTitle(
"chi2");
396 rho_pf->GetXaxis()->SetTitle(
"log10(chi2)");
398 rho_pf->SetTitleOffset(1.3,
"Y");
399 rho_pf->GetYaxis()->SetTitle(
"#rho");
400 rho_pf->SetMarkerStyle(20);
401 rho_pf->SetMarkerColor(1);
402 rho_pf->SetMarkerSize(0.5);
403 rho_pf->SetStats(kFALSE);
405 TH2F*
pf_cc =
new TH2F(
"pf_cc",
"",100,0.,1.,100,0,1.);
406 pf_cc->SetTitleOffset(1.3,
"Y");
407 pf_cc->GetXaxis()->SetTitle(
"network correlation");
408 pf_cc->GetYaxis()->SetTitle(
"penalty");
409 pf_cc->SetMarkerStyle(20);
410 pf_cc->SetMarkerColor(1);
411 pf_cc->SetMarkerSize(0.8);
412 pf_cc->SetStats(kFALSE);
414 TH2F*
pf_vED =
new TH2F(
"pf_vED",
"",100,0.,1.,100,0,1.);
415 pf_vED->SetTitleOffset(1.3,
"Y");
416 pf_vED->GetXaxis()->SetTitle(
"vED consistency");
417 pf_vED->GetYaxis()->SetTitle(
"penalty");
418 pf_vED->SetMarkerStyle(20);
419 pf_vED->SetMarkerColor(1);
420 pf_vED->SetMarkerSize(0.8);
421 pf_vED->SetStats(kFALSE);
423 TH2F*
pf_ed =
new TH2F(
"pf_ed",
"",100,0.,1.,100,0.,1.);
424 pf_ed->SetTitleOffset(1.3,
"Y");
425 pf_ed->GetXaxis()->SetTitle(
"energy disbalance");
426 pf_ed->GetYaxis()->SetTitle(
"penalty");
427 pf_ed->SetMarkerStyle(20);
428 pf_ed->SetMarkerColor(1);
429 pf_ed->SetMarkerSize(0.8);
430 pf_ed->SetStats(kFALSE);
432 TH2F*
rho_L =
new TH2F(
"rho_L",
"",500,0.,5.,500,0,5.);
433 rho_L->SetTitleOffset(1.3,
"Y");
434 rho_L->GetYaxis()->SetTitle(
"log10(average SNR)");
435 rho_L->GetXaxis()->SetTitle(
"log10(2*#rho^2)");
436 rho_L->SetStats(kFALSE);
437 rho_L->SetMarkerStyle(20);
438 rho_L->SetMarkerSize(0.5);
439 rho_L->SetMarkerColor(1);
442 rho_mchirp->SetTitle(
"rho vs Mchirp(after pp cuts)");
443 rho_mchirp->SetTitleOffset(1.3,
"Y");
444 rho_mchirp->GetXaxis()->SetTitle(
"Mchirp");
445 rho_mchirp->GetYaxis()->SetTitle(
"#rho");
446 rho_mchirp->SetMarkerStyle(20);
447 rho_mchirp->SetMarkerColor(1);
448 rho_mchirp->SetMarkerSize(0.5);
449 rho_mchirp->SetStats(kFALSE);
455 for (
int qq=0;qq<sm_ra_dec.size(); qq++) sm_ra_dec.set(qq,0);
473 int hrho_col[2] = {1,3};
474 for(
int j=0; j<2; j++) {
477 hrho[
j]->SetTitleOffset(1.3,
"Y");
479 hrho[
j]->GetXaxis()->SetTitle(_title);
480 hrho[
j]->GetYaxis()->SetTitle(
"events");
481 hrho[
j]->GetYaxis()->SetTitleOffset(2.0);
482 hrho[
j]->GetYaxis()->CenterTitle(
true);
483 hrho[
j]->SetLineColor(hrho_col[j]);
484 hrho[
j]->SetFillColor(hrho_col[j]);
485 hrho[
j]->SetStats(kFALSE);
486 hrho[
j]->SetMinimum(0.5);
489 hrho[0]->SetTitle(
"full events (black), after pp cuts (green)");
493 for(
int j=0; j<nIFO+1; j++){
495 hF[
j] =
new TH1F(ch1,ch1,nBins,T_bgn,T_end);
496 hF[
j]->GetXaxis()->SetTitle(fname);
497 hF[
j]->GetXaxis()->SetTitleSize(0.1);
498 hF[
j]->GetXaxis()->SetLabelSize(0.1);
499 hF[
j]->GetXaxis()->SetTitleOffset(0.9);
501 hF2[
j] =
new TH1F(ch1,ch1,24,32,2048.);
502 hF2[
j]->GetXaxis()->SetTitle(
"frequency, (80 Hz per bin)");
503 hF2[
j]->GetXaxis()->SetTitleSize(0.1);
504 hF2[
j]->GetXaxis()->SetLabelSize(0.1);
505 hF2[
j]->GetXaxis()->SetTitleOffset(0.9);
507 sprintf(_title,
"%s : fraction (%%)",IFO[j-1]);
508 hF[
j]->GetYaxis()->SetTitle(_title);
509 hF2[
j]->GetYaxis()->SetTitle(_title);
511 hF[
j]->GetYaxis()->SetTitleSize(0.12);
512 hF[
j]->GetYaxis()->SetLabelSize(0.1);
513 hF[
j]->GetYaxis()->SetTitleOffset(0.4);
514 hF[
j]->SetStats(kFALSE);
515 hF[
j]->SetLineColor(1);
516 hF[
j]->SetFillColor(1);
517 hF[
j]->SetTitle(
"full events");
518 hF2[
j]->GetYaxis()->SetTitleSize(0.12);
519 hF2[
j]->GetYaxis()->SetLabelSize(0.1);
520 hF2[
j]->GetYaxis()->SetTitleOffset(0.4);
521 hF2[
j]->SetStats(kFALSE);
522 hF2[
j]->SetLineColor(1);
523 hF2[
j]->SetFillColor(1);
524 hF2[
j]->SetTitle(
"full events");
529 for(
int j=0; j<3; j++){
531 hR[
j] =
new TH1F(fname,
"",nBins,T_bgn,T_end);
532 hR[
j]->SetTitleOffset(1.3,
"Y");
534 hR[
j]->GetXaxis()->SetTitle(fname);
535 hR[
j]->GetYaxis()->SetTitle(
"rate, Hz");
536 hR[
j]->GetYaxis()->SetTitleOffset(2.0);
537 hR[
j]->SetTitle(
"full events (black), after pp cuts (green)");
538 hR[
j]->SetStats(kFALSE);
539 hR[
j]->SetLineColor(hR_col[j]);
540 hR[
j]->SetFillColor(hR_col[j]);
549 fprintf(ftrig,
"# -/+ - not passed/passed final selection cuts\n");
550 fprintf(ftrig,
"# 1 - effective correlated amplitude rho \n");
551 fprintf(ftrig,
"# 2 - correlation coefficient 0/1 (1G/2G)\n");
552 fprintf(ftrig,
"# 3 - correlation coefficient 2\n");
553 fprintf(ftrig,
"# 4 - correlation coefficient 3\n");
554 fprintf(ftrig,
"# 5 - correlated amplitude \n");
555 fprintf(ftrig,
"# 6 - time shift \n");
556 fprintf(ftrig,
"# 7 - time super shift \n");
557 fprintf(ftrig,
"# 8 - likelihood \n");
558 fprintf(ftrig,
"# 9 - penalty factor \n");
559 fprintf(ftrig,
"# 10 - energy disbalance \n");
560 fprintf(ftrig,
"# 11 - central frequency \n");
561 fprintf(ftrig,
"# 12 - bandwidth \n");
562 fprintf(ftrig,
"# 13 - duration \n");
563 fprintf(ftrig,
"# 14 - number of pixels \n");
564 fprintf(ftrig,
"# 15 - frequency resolution \n");
565 fprintf(ftrig,
"# 16 - cwb run number \n");
567 for(
int j=0; j<
nIFO; j++)
fprintf(ftrig,
"# %2d - time for %s detector \n",++i,IFO[j]);
568 for(
int j=0; j<
nIFO; j++)
fprintf(ftrig,
"# %2d - sSNR for %s detector \n",++i,IFO[j]);
569 for(
int j=0; j<
nIFO; j++)
fprintf(ftrig,
"# %2d - hrss for %s detector \n",++i,IFO[j]);
570 fprintf(ftrig,
"# %2d - PHI \n",++i);
571 fprintf(ftrig,
"# %2d - THETA \n",++i);
572 fprintf(ftrig,
"# %2d - PSI \n",++i);
575 cout<<
"total events: "<<ntrg<<endl;
576 printf(
"data start GPS: %15.5f \n",sTARt);
577 printf(
"data stop GPS: %15.5f \n",sTOp);
591 cout <<
"cwb_report_prod_2.C : Start event loop ..." << endl;
592 for(
int i=0; i<
ntrg; i++) {
595 if(i%10000==0) cout <<
"cwb_report_prod_2.C : loop 1 - processed: " << i <<
"/" << ntrg << endl;
601 for(
int j=0;j<
nXDQF;j++) {
603 if((XDQF[j].
cat==
CWB_CAT2)&&(!VETO[j])) {countVETO[
j]++;bveto=
true;}
604 if((XDQF[j].
cat==
CWB_CAT3)&&(VETO[j])) {countVETO[
j]++;saveVETO[
j]=0;save_veto=0;}
605 if((XDQF[j].
cat==
CWB_HVETO)&&(VETO[j])) {countVETO[
j]++;saveVETO[
j]=0;save_veto=0;}
609 double WLag = W.lag[
nIFO];
610 double WSLag = W.slag[
nIFO];
613 if((WLag==0)&&(WSLag==0))
continue;
617 if((WLag==0)&&(WSLag==0))
continue;
621 if((WLag==0)&&(WSLag==0))
continue;
631 rUn =
int(W.run+0.1);
632 hR[0]->Fill(
float(W.gps[0]),
Trun.
data[rUn]);
635 hR[1]->Fill(
float(W.gps[0]));
639 for(
int j=0;j<
nFCUT;j++) {
640 if((W.frequency[0]>
lowFCUT[j])&&(W.frequency[0]<=
highFCUT[j])) bcut=
true;
644 vED = W.neted[0]/W.ecor;
649 acor = sqrt(W.ecor/nIFO);
652 if (
T_scc>0) {
if(W.netcc[2] <
T_scc) save =
false;}
653 if (
T_pen>0) {
if(W.penalty <
T_pen) save =
false;}
660 if (
T_efrac>0) {
if(W.chirp[5]*W.chirp[3] + log10(W.size[1])/10<
T_efrac) save =
false;}
669 pf_cc->Fill(xcor,W.penalty);
670 if(W.frequency[0]>0 && W.duration[0]>0) dens->Fill(-log10(W.size[0]*0.5/W.duration[0]/W.frequency[0]));
671 pf_vED->Fill(vED,W.penalty);
672 pf_ed->Fill(vED,W.penalty);
673 hvED->Fill(
fabs(vED));
674 hpf->Fill(W.penalty);
678 for(
int j=0;j<
nXDQF;j++) {
682 for(
int j=0;j<
nXDQF;j++) {
687 int wslag=0; wslag=W.slag[
nIFO];
688 fprintf(ftrig,
"%5.4f %4.3f %4.3f %4.3f %5.2f %7.3f %7d %5.1e %4.3f %4.3f ",
689 rho,xcor,xcor2,xcor3,acor,W.lag[nIFO],wslag,W.likelihood,W.penalty,vED);
690 fprintf(ftrig,
"%4.0f %4.0f %4.3f %3d %3.1d %5d ",
691 W.frequency[0],W.bandwidth[0],
692 W.duration[0],W.size[0],W.rate[0],
int(W.run));
694 for(j=0; j<
nIFO; j++)
fprintf(ftrig,
"%14.4f ",W.time[j]);
695 for(j=0; j<
nIFO; j++)
fprintf(ftrig,
"%5.1e ",W.sSNR[j]);
696 for(j=0; j<
nIFO; j++)
fprintf(ftrig,
"%5.1e ",W.hrss[j]);
697 if(TMath::IsNaN(W.psi[0])) W.psi[0]=-1000;
698 fprintf(ftrig,
"%4.2f %4.2f %4.2f ",W.phi[0], W.theta[0], W.psi[0]);
703 int indx = lagslag[std::make_pair(
int(W.slag[nIFO]),
int(W.lag[nIFO]))];
704 Rlag.
data[indx] += 1.;
715 if(rho >
Xcut.
data[j] && save_veto==1) {y_lag_veto[iy_lag[
int(W.lag[nIFO])-
min_lag]][
j] += 1.;}
728 for(
int i=0;i<
nlag;i++) {
730 delete [] y_lag_veto[
i];
734 for(
int i=0; i<
ntrg; i++) {
737 if(i%10000==0) cout <<
"cwb_report_prod_2.C : loop 2 - processed: " << i <<
"/" << ntrg << endl;
743 sm_ra_dec.add(sm_ra_dec.getSkyIndex(-(W.theta[2]-90.),W.phi[2]),1);
747 vED = W.neted[0]/W.ecor;
752 for(j=0; j<
nIFO; j++) {
753 if(W.snr[j]>mSNR) {detId=j+1; mSNR=W.snr[
j];}
758 rho_L->Fill(log10(2*rho*rho),log10(W.likelihood/nIFO));
760 rhocc->Fill(xcor,rho);
761 ecor_cc->Fill(xcor,W.rho[0]);
762 ECOR_cc->Fill(xcor,W.rho[1]);
763 rho_subnet->Fill(W.netcc[2],rho);
764 rho_vED->Fill(vED,rho);
765 float chi2 = W.penalty>0 ? log10(W.penalty) : 0;
766 if(
TString(analysis)==
"1G") rho_pf->Fill(W.penalty,rho);
767 else rho_pf->Fill(chi2,rho);
769 float Wgps=float(W.gps[0]);
770 float Wfreq=float(W.frequency[0]);
773 hF[detId]->Fill(Wgps,100.);
775 hF2[detId]->Fill(Wfreq,100.);
778 ecor->Fill(sqrt(W.ecor/nIFO));
779 ECOR->Fill(sqrt(W.ECOR/nIFO));
780 Like->Fill(log10(W.likelihood));
782 Mchirp->Fill(W.chirp[1]);
783 McErr->Fill(W.chirp[5]);
784 rho_mchirp->Fill(W.chirp[1],rho);
786 if (!
SAVE[i])
continue;
795 for(
int j=0;j<
nXDQF;j++) {
796 cout <<
"Events vetoed by " << nameVETO[
j].Data() <<
" - " << countVETO[
j] << endl;
803 if(Wlag[nIFO].data[j]>0.) {
804 fprintf(filelags,
"%6.3f %6.3f %d %d\n",
805 Wslag[nIFO].data[j],Wlag[nIFO].data[j],(
int)Tlag.
data[j],(
int)Rlag.
data[j]);
811 if(pp_rho_vs_rate_no_multiplicity) cout <<
"start rate vs mutiplicity ..." << endl;
813 if(pp_rho_vs_rate_no_multiplicity) {
835 if(!out_far.good()) {cout <<
"Error Opening File " << fname << endl;
exit(1);}
839 if(!out_efar.good()) {cout <<
"Error Opening File " << fname << endl;
exit(1);}
845 if(y>=1) out_far << x <<
"\t\t" << y << endl;
846 if(y<1) out_far << x <<
"\t\t" << y << endl;
848 if(y>=1) out_efar << x <<
"\t\t" << y <<
"\t\t" << ex <<
"\t" << ey << endl;
849 if(y<1) out_efar << x <<
"\t\t" << y <<
"\t" << ex <<
"\t" << ey << endl;
884 for (
int k=1; k<
nlag; k++) {
886 if(Tlag[k]>0) rate=y_lag[
k][
j]/Tlag[
k];
893 sigma= sqrt((sum-elag*mean*mean)/elag);
898 cout << fname << endl;
908 for (
int k=1; k<nlag-1; k++) {
910 if(Tlag[k]>0) rate=y_lag_veto[
k][
j]/Tlag[
k];
917 sigma= sqrt((sum-elag*mean*mean)/elag);
922 cout << fname << endl;
928 pf_vED->Draw(
"colz");
930 c1->Update(); c1->SaveAs(fname);
935 c1->Update(); c1->SaveAs(fname);
940 c1->Update(); c1->SaveAs(fname);
945 c1->Update(); c1->SaveAs(fname);
948 c1 =
new TCanvas(
"c",
"C",0,0,1400,400*nIFO);
949 c1->SetBorderMode(0);
951 c1->SetBorderSize(2);
954 c1->SetBottomMargin(0.143939);
955 c1->SetRightMargin(0.1517039);
956 c1->SetTopMargin(0.0772727);
960 for(i=1; i<nIFO+1; i++) {
963 hF[
i]->Divide(hF[0]); hF[
i]->SetMaximum(100.);
964 gStyle->SetOptStat(kFALSE);
965 gPad->SetBottomMargin(0.2);
967 hF[
i]->GetYaxis()->SetLabelSize(0.06);
968 hF[
i]->GetXaxis()->SetTimeDisplay(1);
976 for(i=1; i<nIFO+1; i++) {
979 hF2[
i]->Divide(hF2[0]); hF2[
i]->SetMaximum(100.);
980 gStyle->SetOptStat(kFALSE);
981 gPad->SetBottomMargin(0.2);
983 hF2[
i]->GetYaxis()->SetLabelSize(0.06);
996 c1 =
new TCanvas(
"c",
"C",0,0,1400,1050);
997 c1->SetBorderMode(0);
999 c1->SetBorderSize(2);
1002 c1->SetBottomMargin(0.143939);
1004 c1->SetRightMargin(0.10);
1005 c1->SetLeftMargin(0.15);
1006 c1->SetTopMargin(0.0772727);
1010 hR[1]->Divide(hR[0]);
1012 hR[1]->SetMinimum(1.
e-5);
1013 gStyle->SetOptStat(kFALSE);
1015 hR[1]->GetXaxis()->SetTimeDisplay(1);
1016 hR[1]->GetXaxis()->SetLabelSize(0.04);
1017 for(i=2; i<3; i++) { hR[
i]->Divide(hR[0]); hR[
i]->Draw(
"same"); }
1019 c1->Update(); c1->SaveAs(fname);
1022 c1->SetLogy(kFALSE);
1026 gr2->SetTitle(_title);
1027 gr2->SetLineColor(1);
1028 gr2->SetMarkerStyle(20);
1029 gr2->SetMarkerColor(1);
1030 gr2->SetMarkerSize(1);
1032 gr2->GetHistogram()->SetXTitle(
"lag");
1033 gr2->GetHistogram()->SetYTitle(
"rate, Hz");
1034 gr2->GetHistogram()->GetXaxis()->SetNdivisions(508);
1037 c1->Update(); c1->SaveAs(fname);
1040 c1->SetLogy(kFALSE);
1044 gr2->SetTitle(
"live vs lag");
1045 gr2->SetLineColor(1);
1046 gr2->SetMarkerStyle(20);
1047 gr2->SetMarkerColor(1);
1048 gr2->SetMarkerSize(1);
1050 gr2->GetHistogram()->GetXaxis()->SetNdivisions(508);
1051 gr2->GetHistogram()->SetXTitle(
"lag");
1052 gr2->GetHistogram()->SetYTitle(
"live, days");
1056 c1->Update(); c1->SaveAs(fname);
1060 if(
pp_rho_log) c1->SetLogx(kTRUE);
else c1->SetLogx(kFALSE);
1063 gr1->SetTitle(
"after pp-cuts");
1064 gr1->SetLineColor(1);
1065 gr1->SetMarkerStyle(20);
1066 gr1->SetMarkerColor(1);
1067 gr1->SetMarkerSize(1);
1068 gr1->SetMinimum(0.5/liveTot);
1071 gr1->GetHistogram()->SetXTitle(
"#rho");
1072 gr1->GetHistogram()->SetYTitle(
"rate, Hz");
1075 c1->Update(); c1->SaveAs(fname);
1080 if(
pp_rho_log) c1->SetLogx(kTRUE);
else c1->SetLogx(kFALSE);
1081 TGraphErrors* gr1_veto;
1083 gr1_veto->SetTitle(
"after pp-cuts & vetoes");
1084 gr1_veto->SetLineColor(1);
1085 gr1_veto->SetMarkerStyle(20);
1086 gr1_veto->SetMarkerColor(1);
1087 gr1_veto->SetMarkerSize(1);
1088 gr1_veto->SetMinimum(0.5/liveTot);
1090 gr1_veto->Draw(
"AP");
1091 gr1_veto->GetHistogram()->SetXTitle(
"#rho");
1092 gr1_veto->GetHistogram()->SetYTitle(
"rate, Hz");
1095 c1->Update(); c1->SaveAs(fname);
1098 if(
pp_rho_log) c1->SetLogx(kTRUE);
else c1->SetLogx(kFALSE);
1099 gr1->SetTitle(
"after pp-cuts (black), after pp-cuts & vetoes (red) ");
1101 gr1_veto->SetLineColor(2);
1102 gr1_veto->SetMarkerColor(2);
1103 gr1_veto->Draw(
"P,same");
1105 c1->Update(); c1->SaveAs(fname);
1109 if(pp_rho_vs_rate_no_multiplicity) {
1112 if(
pp_rho_log) c1->SetLogx(kTRUE);
else c1->SetLogx(kFALSE);
1113 TGraphErrors* gr1_multi;
1115 gr1_multi->SetLineColor(1);
1116 gr1_multi->SetMarkerStyle(20);
1117 gr1_multi->SetMarkerColor(1);
1118 gr1_multi->SetMarkerSize(1);
1119 gr1_multi->SetMinimum(0.5/liveTot);
1121 gr1_multi->Draw(
"AP");
1122 gr1_multi->GetHistogram()->SetXTitle(
"#rho");
1123 gr1_multi->GetHistogram()->SetYTitle(
"rate, Hz");
1126 c1->Update(); c1->SaveAs(fname);
1129 if(
pp_rho_log) c1->SetLogx(kTRUE);
else c1->SetLogx(kFALSE);
1131 gr1->SetTitle(
"after pp-cuts (black), no-multiplicity (green) ");
1132 gr1_multi->SetLineColor(kGreen);
1133 gr1_multi->SetMarkerColor(kGreen);
1134 gr1_multi->Draw(
"P,same");
1135 sprintf(fname,
"%s/rate_threshold_multi_compare.eps",
netdir);
1136 c1->Update(); c1->SaveAs(fname);
1140 c1->SetLogy(kFALSE);
1143 c1->Update(); c1->SaveAs(fname);
1149 c1->Update(); c1->SaveAs(fname);
1155 c1->Update(); c1->SaveAs(fname);
1157 gStyle->SetOptStat(kTRUE);
1162 c1->Update(); c1->SaveAs(fname);
1165 if(
pp_rho_log) c1->SetLogx(kTRUE);
else c1->SetLogx(kFALSE);
1168 c1->Update();c1->SaveAs(fname);
1172 hrho[1]->Draw(
"same");
1174 c1->Update(); c1->SaveAs(fname);
1177 gROOT->LoadMacro(wat_dir+
"/tools/cwb/postproduction/burst/h12ascii.C");
1187 c1->Update(); c1->SaveAs(fname);
1192 c1->Update(); c1->SaveAs(fname);
1197 c1->Update(); c1->SaveAs(fname);
1200 c1->SetLogy(kFALSE);
1203 c1->Update(); c1->SaveAs(fname);
1206 c1->SetLogx(kFALSE);
1207 if(
pp_rho_log) c1->SetLogy(kTRUE);
else c1->SetLogy(kFALSE);
1208 rho_subnet->Draw(
"colz");
1210 c1->Update(); c1->SaveAs(fname);
1213 c1->SetLogx(kFALSE);
1214 if(
pp_rho_log) c1->SetLogy(kTRUE);
else c1->SetLogy(kFALSE);
1215 rho_vED->Draw(
"colz");
1217 c1->Update(); c1->SaveAs(fname);
1220 c1->SetLogx(kFALSE);
1221 if(
pp_rho_log) c1->SetLogy(kTRUE);
else c1->SetLogy(kFALSE);
1222 rho_pf->Draw(
"colz");
1224 c1->Update(); c1->SaveAs(fname);
1227 c1->SetLogx(kFALSE);
1228 if(
pp_rho_log) c1->SetLogy(kTRUE);
else c1->SetLogy(kFALSE);
1229 rhocc->Draw(
"colz");
1231 c1->Update(); c1->SaveAs(fname);
1236 c1->Update(); c1->SaveAs(fname);
1241 c1->Update(); c1->SaveAs(fname);
1245 if(
pp_rho_log) c1->SetLogy(kTRUE);
else c1->SetLogy(kFALSE);
1246 rho_mchirp->Draw(
"colz");
1248 c1->Update(); c1->SaveAs(fname);
1249 c1->SetLogx(kFALSE);
1255 gsm_theta_phi->
SetTitle(
"Latitude - Longitude");
1257 gsm_theta_phi->
Draw();
1259 gsm_theta_phi->
GetCanvas()->SaveAs(fname);
1260 delete gsm_theta_phi;
1263 gsm_ra_dec->
SetTitle(
"right ascension - declination");
1267 gsm_ra_dec->
GetHistogram()->GetXaxis()->SetTitle(
"right ascension, degree");
1268 gsm_ra_dec->
GetHistogram()->GetYaxis()->SetTitle(
"declination, degree");
wavearray< double > Yfar(far_rho_bin)
wavearray< double > Ycut_veto(pp_rho_bin)
virtual size_t size() const
void Export(TString fname="")
printf("total live time: non-zero lags = %10.1f \n", liveTot)
wavearray< double > Yerr_veto(pp_rho_bin)
wavearray< double > Trun(500000)
wavearray< double > Yerr_multi(pp_rho_bin)
wavearray< double > yCUT(pp_rho_bin)
wavearray< double > Xfar(far_rho_bin)
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<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t 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
wavearray< int > SAVE(ntrg)
cout<< "Start CWB::Toolbox::getLiveTime : be patient, it takes a while ..."<< endl;liveTot=CWB::Toolbox::getLiveTime(nIFO, live, Trun, Wlag, Wslag, Tlag, Tdlag, cwb_lag_number, cwb_slag_number);std::map< std::pair< int, int >, int > lagslag
wavearray< int > Wsel(ntrg)
void Import(TString umacro="")
void Draw(int dpaletteId=0, Option_t *option="colfz")
wavearray< double > yERR(pp_rho_bin)
wavearray< double > Xerr(pp_rho_bin)
size_t getSkyIndex(double th, double ph)
param: theta param: phi
wavearray< double > Yerr(pp_rho_bin)
ofstream out_lf(fname, ios::out)
wavearray< double > Ycut(pp_rho_bin)
wavearray< double > Xcut(pp_rho_bin)
TIter next(wave.GetListOfBranches())
void SetGalacticDisk(double gpsGalacticDisk=0.0)
wavearray< int > LAG_PASS(ntrg)
wavearray< double > Tdlag
void SetTitle(TString title)
wavearray< double > Ycut_multi(pp_rho_bin)
void SetWorldMap(bool drawWorldMap=true)
fprintf(flive,"nonzero lags live=%10.1f \n", liveTot)
double fabs(const Complex &x)
strcpy(RunLabel, RUN_LABEL)
void h12ascii(TH1 *h, char *fname)
sprintf(fname,"%s/live.txt", netdir)
void set(size_t i, double a)
param: sky index param: value to set
wavearray< double > Wlag[NIFO_MAX+1]
detectorParams detParms[4]
void add(size_t i, double a)
param: sky index param: value to add
virtual void resize(unsigned int)
void SetSingleDetectorMode()
wavearray< double > Wslag[NIFO_MAX+1]