35 gStyle->SetTitleOffset(1.4,
"X");
36 gStyle->SetTitleOffset(1.2,
"Y");
37 gStyle->SetLabelOffset(0.014,
"X");
38 gStyle->SetLabelOffset(0.010,
"Y");
39 gStyle->SetLabelFont(42,
"X");
40 gStyle->SetLabelFont(42,
"Y");
41 gStyle->SetTitleFont(42,
"X");
42 gStyle->SetTitleFont(42,
"Y");
43 gStyle->SetLabelSize(0.03,
"X");
44 gStyle->SetLabelSize(0.03,
"Y");
46 gStyle->SetTitleH(0.050);
47 gStyle->SetTitleW(0.95);
48 gStyle->SetTitleY(0.98);
49 gStyle->SetTitleFont(12,
"D");
50 gStyle->SetTitleColor(kBlue,
"D");
51 gStyle->SetTextFont(12);
52 gStyle->SetTitleFillColor(kWhite);
54 gStyle->SetNumberContours(256);
55 gStyle->SetCanvasColor(kWhite);
56 gStyle->SetStatBorderSize(1);
59 gStyle->SetStatX(0.878);
60 gStyle->SetStatY(0.918);
61 gStyle->SetStatW(0.200);
62 gStyle->SetStatH(0.160);
64 gStyle->SetStatColor(0);
67 gStyle->SetFrameBorderMode(0);
70 cout<<
"cwb_report_sim.C starts..."<<endl;
73 Color_t
colors[
NMDC_MAX] = { 6, 3, 2, 8,43, 7, 8, 4, 4, 2,43, 1, 3, 1, 6, 7,
74 6, 3, 2, 8,43, 7, 8, 4, 4, 2,43, 1, 3, 1, 6, 7,
75 6, 3, 2, 8,43, 7, 8, 4, 4, 2,43, 1, 3, 1, 6, 7,
76 6, 3, 2, 8,43, 7, 8, 4, 4, 2,43, 1, 3, 1, 6, 7 };
77 Style_t markers[
NMDC_MAX]= {20,21,22,23,29,29,21,22,20,21,22,23,29,20,21,22,
78 21,22,23,29,29,21,22,20,21,22,23,29,20,21,22,20,
79 21,22,23,29,29,21,22,20,21,22,23,29,20,21,22,20,
80 21,22,23,29,29,21,22,20,21,22,23,29,20,21,22,20 };
115 imdc_name,imdc_fcentral,imdc_fbandwidth);
118 cout <<
"cwb_report_sim.C : Error - no injection types or bad format" << endl;
119 cout <<
"mdc injection file list : " <<
mdc_inj_file << endl;
120 cout <<
"format must be : mdc_set mdc_type mdc_name mdc_fcentral mdc_fbandwidth" << endl;
121 cout <<
"process terminated !!!" << endl;
126 imdc_flow[
i] = imdc_fcentral[
i]-imdc_fbandwidth[
i];
129 imdc_fhigh[
i] = imdc_fcentral[
i]+imdc_fbandwidth[
i];
135 for(
int j=0;
j<
ninj;
j++) imdc_index[imdc_type[
j]]=
j;
137 cout <<
j <<
" " << imdc_index[
j] << endl;
138 if(imdc_index[
j]>=ninj) {
140 cout <<
"cwb_report_sim.C : Error - mdc type must be < num max type inj = " << ninj << endl;
141 cout <<
"check mdc_type injection file list : " <<
mdc_inj_file << endl;
142 cout <<
"format : mdc_set mdc_type mdc_name mdc_fcentral mdc_fbandwidth" << endl;
143 cout <<
"process terminated !!!" << endl;
152 for(
int j=0;
j<
nset;
j++)
if(imdc_set[
i]==imdc_set_name[
j]) bnew=
false;
153 if(bnew) imdc_set_name[nset++]=imdc_set[
i];
155 cout <<
"nset : " << nset << endl;
158 for(
int j=0;j<
ninj;j++)
if(imdc_set[j]==imdc_set_name[
i]) imdc_iset[
j]=
i;
162 for(
int j=0;j<
nset;j++)
if(imdc_set[
i]==imdc_set_name[j]) imdc_tset[imdc_type[
i]]=
j;
166 cout <<
i <<
" " << imdc_name[
i] <<
" " << imdc_set[
i] <<
" " << imdc_iset[
i]
167 <<
" " << imdc_tset[
i] <<
" " << imdc_type[
i] <<
" " << imdc_index[
i] << endl;
172 if((pp_eff_vs_thr_mode!=
"DISABLED")&&(pp_eff_vs_thr_mode!=
"disabled")) {
176 if(pp_eff_vs_thr_eff.IsFloat()) pp_eff = pp_eff_vs_thr_eff.Atoi();
177 if(pp_eff<0) pp_eff=0;
178 if(pp_eff>100) pp_eff=100;
180 TChain
sim(
"waveburst");
189 if((pp_eff_vs_thr_mode!=
"AUTO")&&(pp_eff_vs_thr_mode!=
"auto")) {
190 sprintf(fname,
"%s",pp_eff_vs_thr_mode.Data());
192 FILE* fp = fopen(fname,
"r");
194 cout <<
"cwb_report_sim.C - Error opening file : " << fname << endl;
201 fscanf(fp,
"%s %d %f",xmdc_name, &xmdc_type, &xfactors);
202 cout <<
i <<
" " << xmdc_name <<
" "<< xmdc_type <<
" "<< xfactors << endl;
203 for(
int j=0;j<
ninj;j++) {
212 cout <<
"cwb_report_sim.C - missing factor for : " << imdc_name[
i] << endl << endl;
227 if((pp_eff_vs_thr_mode==
"AUTO")||(pp_eff_vs_thr_mode==
"auto")) {
234 sprintf(sim_cuts,
"%s && %s && %s",
ch2,fac_cuts,rho_cuts);
236 sprintf(mdc_cuts,
"type[0]==%d && %s",(
int)(imdc_type[
i]+1),fac_cuts);
237 sprintf(sim_cuts,
"%s && type[1]==%lu && %s && %s",
ch2,imdc_type[
i]+1,fac_cuts,rho_cuts);
241 mdc.Draw(
"Entry$",mdc_cuts,
"goff");
242 nmdc =
mdc.GetSelectedRows();
243 sim.Draw(
"Entry$",sim_cuts,
"goff");
244 nsim =
sim.GetSelectedRows();
245 double xeff = double(nsim)/double(nmdc);
246 if(
fabs(100*xeff-pp_eff)<
fabs(100*effXX-pp_eff)) {effXX=xeff;factorXX[
i]=
factors[
k];}
252 cout <<
i <<
"\tCompute efficiency vs threshold @ factor : "
253 << factorXX[
i] <<
"\tmdc : " << imdc_name[
i] << endl;
256 sprintf(fac_cuts,
"abs(factor-%f)<0.1*%f",factorXX[
i],factorXX[i]);
259 sprintf(mdc_cuts,
"type[0]==%d && %s",(
int)(imdc_type[
i]+1),fac_cuts);
262 mdc.Draw(
"Entry$",mdc_cuts,
"goff");
263 nmdc =
mdc.GetSelectedRows();
269 if(pp_eff_vs_thr_rho_min.IsFloat()) rho_min = pp_eff_vs_thr_rho_min.Atof();
272 sprintf(fname,
"%s/eff_%d_threshold_%s.txt",
netdir,pp_eff,imdc_name[
i]);
273 cout << fname << endl;
276 if (!fout.good()) {cout <<
"Error Opening File : " << fname << endl;
exit(1);}
286 sprintf(sim_cuts,
"%s && type[1]==%lu && %s && rho[%d]>%f",
ch2,imdc_type[
i]+1,fac_cuts,
pp_irho,rho);
289 sim.Draw(
"Entry$",sim_cuts,
"goff");
290 nsim =
sim.GetSelectedRows();
291 double xeff = double(nsim)/double(nmdc);
294 fout << rho <<
" " << xeff << endl;
299 sprintf(fname,
"%s/eff_%d_threshold_factors.txt",
netdir,pp_eff,imdc_name[
i]);
300 cout << fname << endl;
301 FILE* fp = fopen(fname,
"w");
303 cout <<
"cwb_report_sim.C - Error opening file : " << fname << endl;
307 fprintf(fp,
"%s %d %g\n",imdc_name[
i], (
int)imdc_type[i], factorXX[i]);
312 pp_eff_vs_thr_quit.ToUpper();
313 if(pp_eff_vs_thr_quit==
"TRUE") {
314 cout <<
"efficiency " << pp_eff <<
" vs threshold terminated, exit" << endl;
325 for (
int iset=0;iset<
nset;iset++) {
327 TChain
sim(
"waveburst");
341 cout <<
"cwb_report_sim.C Error : Leaf " << veto_name
342 <<
" is not present in tree" << endl;
348 for(
int j=0;j<
nVDQF;j++) {
353 if(bifo) XDQF[nXDQF++]=
VDQF[
j];
355 for(
int j=0;j<
nXDQF;j++) {
360 for(
int j=0;j<
nXDQF;j++) {
363 ww.fChain->SetBranchAddress(veto_name,&VETO[j]);
371 inf = log10(pp_factor2distance/
factors[nfactor-1]);
372 sup = log10(pp_factor2distance/
factors[0]);
375 TF1 *gfit=
new TF1(
"logNfit",
logNfit,pow(10.0,inf),pow(10.0,sup),5);
376 gfit->SetParameters((inf+sup)/2.,0.3,1.,1.);
377 gfit->SetParNames(
"hrss50",
"sigma",
"betam",
"betap");
378 gfit->FixParameter(4,pp_factor2distance);
379 gfit->SetParLimits(0,inf,sup);
384 gfit->SetParLimits(1,0.05,10.);
385 gfit->SetParLimits(2,0.05,3.);
387 gfit->SetParLimits(1,0.08,10.);
388 gfit->SetParLimits(2,0.00,3.);
390 gfit->SetParLimits(3,0.0,2.50);
393 TLegend *legend =
new TLegend(0.53,0.17,0.99,0.7,
"",
"brNDC");
394 legend->SetLineColor(1);
395 legend->SetTextSize(0.033);
396 legend->SetBorderSize(1);
397 legend->SetLineStyle(1);
398 legend->SetLineWidth(1);
399 legend->SetFillColor(0);
400 legend->SetFillStyle(1001);
404 double a,b,
vED,enrg[5];
418 int nmdc =
mdc.GetEntries();
420 cout<<
" total injections: " << nmdc <<endl;
422 for(
int i=0; i<ninj+1; i++) {
423 if((i<ninj)&&(imdc_iset[i]!=iset))
continue;
427 hT[
i]->SetTitleOffset(1.3,
"Y");
428 hT[
i]->GetXaxis()->SetTitle(
"detected-injected time, sec");
429 hT[
i]->GetYaxis()->SetTitle(
"events");
430 sprintf(fname,
"%s",imdc_name[i]);
431 hT[
i]->SetTitle(fname);
434 hF[
i] =
new TH1F(fname,fname,512,imdc_flow[i],imdc_fhigh[i]);
435 hF[
i]->SetTitleOffset(1.3,
"Y");
436 if(i<ninj) hF[
i]->GetXaxis()->SetTitle(
"frequency, Hz");
437 else hF[
i]->GetXaxis()->SetTitle(
"detected-injected, Hz");
438 hF[
i]->GetYaxis()->SetTitle(
"events");
439 sprintf(fname,
"%s",imdc_name[i]);
440 hF[
i]->SetTitle(fname);
442 for(
int j=0; j<
nIFO; j++) {
443 sprintf(fname,
"%s_hrss_%d",IFO[j],i);
444 HH[
j][
i] =
new TH2F(fname,
"",500,-23.,-19.,500,-23.,-19.);
445 HH[
j][
i]->SetTitleOffset(1.7,
"Y");
446 HH[
j][
i]->SetTitleOffset(1.7,
"Y");
447 HH[
j][
i]->GetXaxis()->SetTitle(
"injected log10(hrss)");
448 HH[
j][
i]->GetYaxis()->SetTitle(
"detected log10(hrss)");
449 HH[
j][
i]->SetStats(kFALSE);
450 sprintf(fname,
"%s: %s",IFO[j],imdc_name[i]);
451 HH[
j][
i]->SetTitle(fname);
452 HH[
j][
i]->SetMarkerColor(2);
454 for(
int j=0; j<
nIFO; j++) {
455 sprintf(fname,
"%s_nre_%d",IFO[j],i);
456 NRE[
j][
i] =
new TH2F(fname,
"",500,0.,7.,500,-3,1.);
457 NRE[
j][
i]->SetTitleOffset(1.7,
"Y");
458 NRE[
j][
i]->SetTitleOffset(1.7,
"Y");
459 NRE[
j][
i]->GetXaxis()->SetTitle(
"injected log10(snr^{2})");
460 NRE[
j][
i]->GetYaxis()->SetTitle(
"log10(normalized residual energy)");
461 NRE[
j][
i]->SetStats(kFALSE);
462 sprintf(fname,
"%s: %s",IFO[j],imdc_name[i]);
463 NRE[
j][
i]->SetTitle(fname);
464 NRE[
j][
i]->SetMarkerColor(2);
468 TH1F* hcor =
new TH1F(
"hcor",
"",100,0,1.);
469 hcor->SetTitleOffset(1.3,
"Y");
470 hcor->GetXaxis()->SetTitle(
"network correlation (ecor)");
471 hcor->GetYaxis()->SetTitle(
"events");
472 hcor->SetStats(kFALSE);
474 TH1F* hcc =
new TH1F(
"hcc",
"",100,0,1.);
475 hcc->SetTitleOffset(1.3,
"Y");
476 hcc->GetXaxis()->SetTitle(
"network correlation");
477 hcc->GetYaxis()->SetTitle(
"events");
478 hcc->SetStats(kFALSE);
480 TH1F* hrho =
new TH1F(
"hrho",
"",500,0,10.);
481 hrho->SetTitleOffset(1.3,
"Y");
482 hrho->GetXaxis()->SetTitle(
"log10(#rho^2)");
483 hrho->GetYaxis()->SetTitle(
"events");
485 TH1F* eSNR =
new TH1F(
"eSNR",
"",500,0,10.);
486 eSNR->SetTitleOffset(1.3,
"Y");
487 eSNR->GetXaxis()->SetTitle(
"log10(SNR/det)");
488 eSNR->GetYaxis()->SetTitle(
"events");
489 eSNR->SetStats(kFALSE);
491 TH2F*
rhocc =
new TH2F(
"rhocc",
"",100,0.,1.,500,0,10.);
492 rhocc->SetTitleOffset(1.3,
"Y");
493 rhocc->GetXaxis()->SetTitle(
"network correlation");
494 rhocc->GetYaxis()->SetTitle(
"log10(rho^2)");
495 rhocc->SetStats(kFALSE);
497 TH2F* skyc =
new TH2F(
"skyc",
"",180,-180,180.,90,-90.,90.);
498 skyc->SetTitleOffset(1.3,
"Y");
499 skyc->GetXaxis()->SetTitle(
"#Delta#phi, deg.");
500 skyc->GetYaxis()->SetTitle(
"#Delta#theta, deg.");
501 skyc->SetStats(kFALSE);
503 TH1F*
hpf =
new TH1F(
"hpf",
"",100,0.,1.);
504 hpf->SetTitleOffset(1.3,
"Y");
506 hpf->GetXaxis()->SetTitle(
"penalty factor");
507 hpf->GetYaxis()->SetTitle(
"events");
509 TH1F*
hvED =
new TH1F(
"hvED",
"",200,0,2.);
510 hvED->SetTitleOffset(1.3,
"Y");
512 hvED->GetXaxis()->SetTitle(
"hvED consistency");
513 hvED->GetYaxis()->SetTitle(
"events");
515 TH2F*
pf_cc =
new TH2F(
"pf_cc",
"",100,0.,1.,100,0,1.);
516 pf_cc->SetTitleOffset(1.3,
"Y");
517 pf_cc->GetXaxis()->SetTitle(
"network correlation");
518 pf_cc->GetYaxis()->SetTitle(
"penalty");
519 pf_cc->SetMarkerStyle(20);
520 pf_cc->SetMarkerColor(1);
521 pf_cc->SetMarkerSize(0.8);
522 pf_cc->SetStats(kFALSE);
524 TH2F*
pf_vED =
new TH2F(
"pf_vED",
"",100,0.,1.,100,0,1.);
525 pf_vED->SetTitleOffset(1.3,
"Y");
526 pf_vED->GetXaxis()->SetTitle(
"H1H2 consistency");
527 pf_vED->GetYaxis()->SetTitle(
"penalty");
528 pf_vED->SetMarkerStyle(20);
529 pf_vED->SetMarkerColor(1);
530 pf_vED->SetMarkerSize(0.8);
531 pf_vED->SetStats(kFALSE);
537 for(
int i=0; i<
nIFO; i++) {
539 sprintf(fname,
"ifoT%s",IFO[i]);
540 ifoT[
i] =
new TH1F(fname,
"",250,-0.05,0.05);
541 ifoT[
i]->SetTitleOffset(1.3,
"Y");
542 ifoT[
i]->GetXaxis()->SetTitle(
"detected-injected time, sec");
543 ifoT[
i]->GetYaxis()->SetTitle(
"events");
544 ifoT[
i]->SetTitle(IFO[i]);
546 sprintf(fname,
"hrs1%s",IFO[i]);
547 hrs1[
i] =
new TH2F(fname,
"",500,-23.,-19.,500,-23.,-19.);
548 hrs1[
i]->SetTitleOffset(1.3,
"Y");
549 hrs1[
i]->GetXaxis()->SetTitle(
"detected");
550 hrs1[
i]->GetYaxis()->SetTitle(
"injected");
551 hrs1[
i]->SetStats(kFALSE);
552 hrs1[
i]->SetTitle(IFO[i]);
553 hrs1[
i]->SetMarkerColor(2);
555 sprintf(fname,
"hrs2%s",IFO[i]);
556 hrs2[
i] =
new TH2F(fname,
"",500,-23.,-19.,500,-23.,-19.);
557 hrs2[
i]->SetTitleOffset(1.3,
"Y");
558 hrs2[
i]->SetMarkerColor(2);
559 hrs2[
i]->SetStats(kFALSE);
571 for(
int i=0; i<
nmdc; i++) {
576 if(imdc_tset[mm.type-1] != iset)
continue;
577 indx = imdc_index[mm.type-1];
579 if(mm.strain<hrss_min[indx])
if(mm.strain>0) hrss_min[indx]=mm.strain;
580 if(mm.strain>hrss_max[indx]) hrss_max[indx]=mm.strain;
585 if(hrss_min[i]==1e20)
continue;
586 if(hrss_max[i]==0)
continue;
587 vhrss[
i][0]=hrss_min[
i];
588 vhrss[
i][nhrss]=hrss_max[
i];
589 double fhrss=exp(
log(vhrss[i][nhrss]/vhrss[i][0])/nhrss);
590 for(
int j=1;j<=nhrss;j++) vhrss[i][j]=fhrss*vhrss[i][j-1];
591 for(
int j=0;j<nhrss;j++) sMDC[i][j]=(vhrss[i][j+1]+vhrss[i][j])/2.;
593 for(
int j=0;j<nhrss;j++)
for(
int k=0; k<
NMDC_MAX; k++) {nMDC[
k].
append(1.);}
597 for(
int i=0; i<
nmdc; i++) {
603 if(imdc_tset[mm.type-1] != iset)
continue;
604 indx = imdc_index[mm.type-1];
609 for(
int j=0; j<
n; j++) {
612 if(pp_factor2distance) {
613 if(
fabs(sMDC[0].data[j]-mm.factor)<0.1*mm.factor) {
614 save =
false; nMDC[indx].
data[
j]+=1.; j=
n;
617 if(
fabs(sMDC[0].data[j]-mm.strain)<0.1*mm.strain) {
618 save =
false; nMDC[indx].
data[
j]+=1.; j=
n;
623 if(
fabs(sMDC[0].data[j]-mm.factor)<0.1*mm.factor) {
624 save =
false; nMDC[indx].
data[
j]+=1.; j=
n;
629 if(mm.strain>vhrss[indx][j] && mm.strain<=vhrss[indx][j+1]) {
630 nMDC[indx].
data[
j]+=1.;
631 mhrss[indx][
j].
append(mm.strain);
641 if(pp_factor2distance) {
655 nMDC[indx].
data[
j]+=1.;
659 cout <<
"cwb_report_sim.C Error : number of detected strains " << sMDC[0].
size()
660 <<
" is greater of nfactor declared in the user_parameters.C : " << nfactor << endl;
668 if(sMDC[k].
size()==0) {
delete [] mhrss[
k];
continue;}
670 for(
int i=0;i<nhrss;i++) {
672 int* hrss_index =
new int[
K];
673 if(K>0) TMath::Sort(K,mhrss[k][i].data,hrss_index,kFALSE);
675 if(K>0) sMDC[
k][
i] = mhrss[
k][
i].
data[hrss_index[K/2]];
676 delete [] hrss_index;
700 if(sMDC[k].
size()==0)
continue;
701 mdc_index[
k] =
new int[sMDC[
k].
size()];
702 TMath::Sort((
int)sMDC[k].
size(),sMDC[k].data,mdc_index[k],kFALSE);
713 cout<<
"total events: " << ntrg << endl;
717 bool ifar_exists=
false;
718 TIter
next(
ww.fChain->GetListOfBranches());
719 while ((branch=(TBranch*)
next())) {
720 if(
TString(
"ifar").CompareTo(branch->GetName())==0) ifar_exists=
true;
723 if (ifar_exists)
ww.fChain->SetBranchAddress(
"ifar",&ifar);
726 for(
int i=0; i<
ntrg; i++) {
729 if(
ww.type[1]<1)
continue;
732 if(imdc_tset[
ww.type[1]-1] != iset)
continue;
737 for(
int j=0;j<
nXDQF;j++) {
738 if((XDQF[j].
cat==
CWB_CAT2)&&(!VETO[j])) bveto=
true;
744 vED =
ww.neted[0]/
ww.ecor;
748 if (
T_ifar>0)
if(ifar_exists) {
if(ifar <
T_ifar)
continue;}
752 for(
int j=0;j<
nFCUT;j++) {
764 indx = imdc_index[
ww.type[1]-1];
767 null = nill = etot = 0.;
768 for(
int j=0; j<
nIFO; j++) {
776 for(
int j=0; j<
nIFO; j++) {
777 if(a2or > etot-enrg[j]) a2or = etot-enrg[
j];
781 acor = sqrt(
ww.ECOR/nIFO);
789 rhocc->Fill(xcor,log10(rho*rho));
790 pf_cc->Fill(xcor,
ww.penalty);
792 if(xcor<
T_cor)
continue;
794 pf_vED->Fill(vED,
ww.penalty);
795 hvED->Fill(
fabs(vED));
796 hpf->Fill(
ww.penalty);
798 if(
ww.netcc[3] <
T_scc)
continue;
799 if(
ww.penalty <
T_pen)
continue;
803 if(rho >
T_cut) save=
true;
806 hrho->Fill(log10(rho*rho));
807 eSNR->Fill(log10(
ww.likelihood/nIFO));
809 a =
ww.phi[0]-
ww.phi[1];
810 if(a> 180) a = 360-
a;
811 if(a< -180) a = -360-
a;
812 b =
ww.theta[0]-
ww.theta[1];
814 if(b< -90) b = -180-b;
817 for(
int j=0; j<
nIFO; j++) {
819 HH[
j][indx]->Fill(log10(
ww.hrss[j+nIFO]),log10(
ratio_hrss*
ww.hrss[j]));
821 double nre =
ww.iSNR[
j] ? (
ww.iSNR[
j]+
ww.oSNR[
j]-2*
ww.ioSNR[
j])/
ww.iSNR[j] : 10;
823 NRE[
j][
ninj]->Fill(log10(
ww.iSNR[j]),log10(nre));
824 NRE[
j][indx]->Fill(log10(
ww.iSNR[j]),log10(nre));
826 ifoT[
j]->Fill(
ww.time[j]-
ww.time[j+nIFO]);
829 hT[indx]->Fill(
ww.time[0]-
ww.time[nIFO]);
831 hT[
ninj]->Fill(
ww.time[0]-
ww.time[nIFO]);
832 hF[
ninj]->Fill(
float(
ww.
frequency[0]-(imdc_flow[indx]+imdc_fhigh[indx])/2));
836 for(
int j=0; j<
n; j++) {
839 if(pp_factor2distance) {
840 if(
fabs(sMDC[0].data[j]-
ww.factor)<0.1*
ww.factor) {
841 nTRG[indx].
data[
j]+=1.; j=
n;
844 if(
fabs(sMDC[0].data[j]-
ww.strain[1])<0.1*
ww.strain[1]) {
845 nTRG[indx].
data[
j]+=1.; j=
n;
850 if(
fabs(sMDC[0].data[j]-
ww.factor)<0.1*
ww.factor) {
851 nTRG[indx].
data[
j]+=1.; j=
n;
856 if(
ww.strain[1]>vhrss[indx][j] &&
ww.strain[1]<=vhrss[indx][j+1]) {
857 nTRG[indx].
data[
j]+=1.; j=
n;
865 for(
int k=0; k<
NMDC_MAX; k++)
for(
int i=0;i<sMDC[
k].size();i++) sMDC[k].data[i] = pp_factor2distance/sMDC[k].data[i];
869 cout <<
"INJ HRSS NUMBER: " << n << endl;
873 for(
int i=0; i<
ninj; i++) {
874 if(imdc_iset[i]!=iset)
continue;
876 in = fopen(f_file,
"w");
877 if(in==NULL) {cout <<
"cwb_report_sim : Error opening file " << f_file << endl;
exit(1);}
878 for(
int k=0; k<
n; k++){
879 int j=mdc_index[
i][
k];
880 if(nMDC[i].data[j] <= 0.)
continue;
883 sqrt(eff[i].data[j]*(1-eff[i].data[j])/nMDC[i].data[j]);
884 fprintf(in,
"%e %i %i %.3f\n", sMDC[i].data[j], (
int)nTRG[i].data[j], (
int)nMDC[i].data[j], eff[i].data[j]);
889 TCanvas *
c1 =
new TCanvas(
"c",
"C",0,0,800,800);
890 c1->SetBorderMode(0);
892 c1->SetBorderSize(2);
896 c1->SetLeftMargin(0.137039);
897 c1->SetRightMargin(0.117039);
898 c1->SetTopMargin(0.0772727);
899 c1->SetBottomMargin(0.143939);
900 c1->ToggleEventStatus();
902 for(
int i=0; i<
ninj; i++) {
903 if(imdc_iset[i]!=iset)
continue;
905 fprintf(evt_all,
"%2d %s ", i,imdc_name[i]);
913 fprintf(evt_all,
"%3.2f %3.2f ", hF[i]->GetMean(), hF[i]->GetRMS());
915 c1->Update(); c1->SaveAs(fname);
920 fprintf(evt_all,
"%3.6f %3.6f ", hT[i]->GetMean(), hT[i]->GetRMS());
921 c1->Update(); c1->SaveAs(fname);
925 for(
int j=0;j<
nIFO;j++){
928 sprintf(fname,
"%s/hrss_%s_%s.gif",
netdir,IFO[j],imdc_name[i]);
929 c1->Update(); c1->SaveAs(fname);
933 sprintf(fname,
"%s/nre_%s_%s.gif",
netdir,IFO[j],imdc_name[i]);
934 c1->Update(); c1->SaveAs(fname);
944 sprintf(fname,
"%s/netcor_%s.gif",
netdir,imdc_set_name[iset].Data());
945 c1->Update(); c1->SaveAs(fname);
952 sprintf(fname,
"%s/cc_%s.gif",
netdir,imdc_set_name[iset].Data());
953 c1->Update(); c1->SaveAs(fname);
959 sprintf(fname,
"%s/rho_%s.gif",
netdir,imdc_set_name[iset].Data());
960 c1->Update(); c1->SaveAs(fname);
964 sprintf(fname,
"%s/penalty_%s.gif",
netdir,imdc_set_name[iset].Data());
965 c1->Update(); c1->SaveAs(fname);
969 sprintf(fname,
"%s/average_snr_%s.gif",
netdir,imdc_set_name[iset].Data());
970 c1->Update(); c1->SaveAs(fname);
976 sprintf(fname,
"%s/rho_cc_%s.gif",
netdir,imdc_set_name[iset].Data());
977 c1->Update(); c1->SaveAs(fname);
982 sprintf(fname,
"%s/penalty_cc_%s.gif",
netdir,imdc_set_name[iset].Data());
983 c1->Update(); c1->SaveAs(fname);
988 sprintf(fname,
"%s/coordinate_%s.gif",
netdir,imdc_set_name[iset].Data());
989 c1->Update(); c1->SaveAs(fname);
994 pf_vED->Draw(
"colz");
995 sprintf(fname,
"%s/penalty_vED_%s.gif",
netdir,imdc_set_name[iset].Data());
996 c1->Update(); c1->SaveAs(fname);
1002 sprintf(fname,
"%s/vED_%s.gif",
netdir,imdc_set_name[iset].Data());
1003 c1->Update(); c1->SaveAs(fname);
1004 c1->SetLogy(kFALSE);
1008 c1 =
new TCanvas(
"c",
"C",0,0,1000,800);
1009 c1->SetBorderMode(0);
1010 c1->SetFillColor(0);
1011 c1->SetBorderSize(2);
1012 c1->SetLogx(kFALSE);
1015 c1->SetRightMargin(0.0517039);
1016 c1->SetTopMargin(0.0772727);
1017 c1->SetBottomMargin(0.143939);
1018 c1->ToggleEventStatus();
1020 double chi2, hrss50,
sigma, betam, betap, hrssEr;
1024 sprintf(fname,
"%s/fit_parameters_%s.txt",
netdir,imdc_set_name[iset].Data());
1025 in = fopen(fname,
"w");
1028 for(
int i=0; i<
ninj; i++)
if(imdc_iset[i]==iset) iset_size++;
1029 double hleg = 0.18+iset_size*0.06;
1031 if(pp_show_eff_fit_curve) legend =
new TLegend(0.45,0.18,0.99,hleg,NULL,
"brNDC");
1032 else legend =
new TLegend(0.55,0.18,0.99,hleg,NULL,
"brNDC");
1034 for(
int i=0; i<
ninj; i++) {
1035 if(imdc_iset[i]!=iset)
continue;
1039 EFF[
i]=
new TGraphErrors(n,sMDC[i].data,eff[i].data,e00[i].data,err[i].data);
1040 if(
simulation==1 && pp_factor2distance) EFF[
i]->GetXaxis()->SetTitle(
"distance (Kpc)");
1041 else if(
simulation==2) EFF[
i]->GetXaxis()->SetTitle(
"snr");
1042 else EFF[
i]->GetXaxis()->SetTitle(
"hrss, #frac{strain}{#sqrt{Hz}}");
1043 EFF[
i]->GetYaxis()->SetTitle(
"Efficiency ");
1044 EFF[
i]->GetXaxis()->SetTitleOffset(1.7);
1045 EFF[
i]->GetXaxis()->SetLabelSize(0.04);
1046 EFF[
i]->GetYaxis()->SetTitleOffset(1.3);
1047 EFF[
i]->GetYaxis()->SetRangeUser(0,1);
1048 EFF[
i]->GetYaxis()->SetLabelSize(0.04);
1052 gfit->SetLineColor(colors[i]);
1057 gfit->SetNpx(100000);
1058 if(pp_show_eff_fit_curve) EFF[
i]->Fit(
"logNfit");
else EFF[
i]->Fit(
"logNfit",
"N");
1059 chi2=gfit->GetChisquare();
1060 hrss50=pow(10.,gfit->GetParameter(0));
1061 hrssEr=(pow(10.,gfit->GetParError(0))-1.)*hrss50;
1062 sigma=gfit->GetParameter(1);
1063 betam=gfit->GetParameter(2);
1064 betap=gfit->GetParameter(3);
1066 printf(
"%d %s %E %E +- %E %E %E %E\n",
1067 i,imdc_name[i],chi2,hrss50,hrssEr,sigma,betam,betap);
1068 fprintf(in,
"%2d %9.3E %9.3E +- %9.3E %9.3E %9.3E %9.3E %s\n",
1069 i,chi2,hrss50,hrssEr,sigma,betam,betap,imdc_name[i]);
1070 fprintf(in_all,
"%2d %9.3E %9.3E +- %9.3E %9.3E %9.3E %9.3E %s\n",
1071 i,chi2,hrss50,hrssEr,sigma,betam,betap,imdc_name[i]);
1073 sprintf(fname,
"%s, hrss50=%5.2E", imdc_name[i],hrss50);
1074 EFF[
i]->SetTitle(fname);
1076 gfit->SetLineColor(colors[i]);
1077 EFF[
i]->SetLineColor(colors[i]);
1078 EFF[
i]->SetMarkerStyle(markers[i]);
1079 EFF[
i]->SetMarkerColor(colors[i]);
1080 EFF[
i]->SetMarkerSize(2);
1082 legend->SetBorderSize(1);
1083 legend->SetTextAlign(22);
1084 legend->SetTextFont(12);
1085 legend->SetLineColor(1);
1086 legend->SetLineStyle(1);
1087 legend->SetLineWidth(1);
1088 legend->SetFillColor(0);
1089 legend->SetFillStyle(1001);
1090 legend->SetTextSize(0.04);
1091 legend->SetLineColor(kBlack);
1092 legend->SetFillColor(kWhite);
1094 if(pp_show_eff_fit_curve)
sprintf(fname,
"%s, %5.2E",imdc_name[i], hrss50);
1095 else sprintf(fname,
"%s",imdc_name[i]);
1096 legend->AddEntry(EFF[i], fname,
"lp");
1097 if(pp_show_eff_fit_curve) EFF[
i]->Draw(
"AP");
else EFF[
i]->Draw(
"APL");
1100 c1->SetLogx(kFALSE);
1106 sprintf(fname,
"%s/eff_%s.eps",
netdir,imdc_set_name[iset].Data());
1107 sprintf(gname,
"%s/eff_%s.gif",
netdir,imdc_set_name[iset].Data());
1108 TPostScript epsfile(fname,113);
1111 if(
simulation==2) c1->SetLogx(kFALSE);
else c1->SetLogx();
1112 c1->SetLogy(kFALSE);
1116 TMultiGraph*
mg=
new TMultiGraph();
1117 sprintf(ptitle,
"%s rho=%3.2f, cc=%3.2f",imdc_set_name[iset].Data(),
T_cut,
T_cor);
1118 mg->SetTitle(ptitle);
1119 double xmin=1.79769e+308;
1121 for(
int i=0;i<
ninj;i++) {
1122 if(imdc_iset[i]!=iset)
continue;
1124 if(sMDC[i].data[mdc_index[i][0]]<xmin) xmin=sMDC[
i].data[mdc_index[
i][0]];
1125 if(sMDC[i].data[mdc_index[i][sMDC[i].
size()-1]]>xmax) xmax=sMDC[
i].data[mdc_index[
i][sMDC[
i].size()-1]];
1127 if(pp_show_eff_fit_curve) mg->Draw(
"AP");
else mg->Draw(
"APL");
1128 if(
simulation==1 && pp_factor2distance) mg->GetHistogram()->GetXaxis()->SetTitle(
"distance (Kpc)");
1129 else if(
simulation==2) mg->GetHistogram()->GetXaxis()->SetTitle(
"snr");
1130 else mg->GetHistogram()->GetXaxis()->SetTitle(
"hrss, #frac{strain}{#sqrt{Hz}}");
1131 mg->GetHistogram()->GetYaxis()->SetTitle(
"Efficiency ");
1132 mg->GetHistogram()->GetXaxis()->SetTitleOffset(1.7);
1133 mg->GetHistogram()->GetXaxis()->SetLabelSize(0.04);
1134 mg->GetHistogram()->GetYaxis()->SetTitleOffset(1.3);
1135 mg->GetHistogram()->GetXaxis()->SetRangeUser(xmin,xmax);
1136 mg->GetHistogram()->GetYaxis()->SetRangeUser(0,1);
1137 mg->GetHistogram()->GetYaxis()->SetLabelSize(0.04);
1138 legend->Draw(
"SAME");
1139 c1->Update(); epsfile.Close();
1141 sprintf(cmd,
"convert %s %s",fname,gname);cout<<cmd<<endl;
1145 for (
int i=0;i<ninj+1;i++) {
1146 if((i<ninj)&&(imdc_iset[i]!=iset))
continue;
1147 for(
int j=0;j<
nIFO;j++)
if(HH[j][i])
delete HH[
j][
i];
1148 for(
int j=0;j<
nIFO;j++)
if(NRE[j][i])
delete NRE[
j][
i];
1149 if(hT[i])
delete hT[
i];
1150 if(hF[i])
delete hF[
i];
1164 for (
int j=0;j<
nIFO;j++) {
1172 if(sMDC[k].
size()==0)
continue;
1173 delete [] mdc_index[
k];
1176 for(
int i=0;i<
NMDC_MAX;i++)
if(EFF[i]) {
delete EFF[
i];EFF[
i]=NULL;}
size_t append(const wavearray< DataType_t > &)
virtual size_t size() const
void Export(TString fname="")
printf("total live time: non-zero lags = %10.1f \n", liveTot)
fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(), x.size()/x.rate(), x.rate())
wavearray< double > a(hp.size())
void Import(TString umacro="")
size_t imdc_iset[NMDC_MAX]
size_t imdc_index[NMDC_MAX]
TIter next(twave->GetListOfBranches())
char imdc_name[NMDC_MAX][128]
double fabs(const Complex &x)
TString pp_eff_vs_thr_quit
strcpy(RunLabel, RUN_LABEL)
size_t imdc_type[NMDC_MAX]
detectorParams detParms[4]
virtual void resize(unsigned int)
double imdc_fcentral[NMDC_MAX]
char imdc_set[NMDC_MAX][128]
void SetSingleDetectorMode()
double imdc_fbandwidth[NMDC_MAX]
sprintf(fname,"%s/eff_%d_threshold_factors.txt", netdir, pp_eff, imdc_name[i])