3 #pragma GCC system_header
10 #include "TObjArray.h"
11 #include "TObjString.h"
14 #include "TGraphAsymmErrors.h"
103 #define nTRIALS 100 // number of trials for noise statistic (GetNoiseStatistic)
105 #define PE_INDEX "$HOME_CWB/www/pe_index.html" // html template used for pe index report page
108 #define PLOT_TYPE "png"
112 #define PLOT_MEDIAN // plot median,lower50 boud,upper50 bound,lower90 boud,upper90 bound
116 #define SKYMASK_RADIUS 0.1 // used to define skymask
118 #define MAX_TRIALS 1000
120 #define EFRACTION_CUT 0.98 // energy threshold for time and frequency cuts
124 #define PE_TRIALS 100 // number of trials
125 #define PE_SIGNAL 0 // signal used for trials : 0/1/2 reconstructed/injected/(original=reconstructed+noise)
127 #define PE_NOISE 1 // noise used for trials : 0/1/2
131 #define PE_AMP_CAL_ERR 0 // max percentage of amplitude miscalibration if>0 (uniform in -max,+max) else (gaus(max))
132 #define PE_PHS_CAL_ERR 0 // max phase (degrees) of phase miscalibration if>0 (uniform in -max,+max) else (gaus(max))
133 #define PE_SKYMASK 0 // skymask used for trials : 0(def)/1/2
137 #define PE_SEED 0 // 0 : seed used by PE for ramdom generation
138 #define PE_GPS 0 // default 0 - if >0 only gps +/- iwindow is analyzed
140 #define PE_MULTITASK false // if enabled, PE trials are executed by different jobs in multitask : only simulation=0 and single event (gps>0)
142 #define NOISE_SIGMA 1.0 // noise sigma used in AddNoiseAndCalErrToSparse
144 #define PE_CED_DUMP -1 // dump CED at gLRETRY=PE_CED_DUMP
146 #define PE_CED_TFMAP true // Shows the Time-Frequency Maps -> Spectrograms, Scalograms, TF Like,Null
147 #define PE_CED_RDR true // Shows Reconstructed Detector Response tab
148 #define PE_CED_PSD true // Shows Power Spectral Density tab
149 #define PE_CED_SKYMAP true // Shows SkyMaps
150 #define PE_CED_REC true // Shows Reconstructed Waveforms/Instantaneous-Frequency with errors
151 #define PE_CED_INJ true // Shows Injected vs Reconstructed Waveforms/Instantaneous-Frequency with errors
152 #define PE_CED_rINJ false // Shows Reduced Injected vs Reconstructed Waveforms with errors
153 #define PE_CED_CM true // Shows Chirp Mass Value/Error Distributions
154 #define PE_CED_DISTR false // Shows Residuals Distributions
155 #define PE_CED_NULL true // Shows Null Pixels Distributions
156 #define PE_CED_PCA false // Shows PCA likelihood TF plot
158 #define PE_OUTPUT_INJ false // save injected into the output root file
159 #define PE_OUTPUT_REC false // save reconstructed into the output root file
160 #define PE_OUTPUT_WHT false // save whitened data into the output root file
161 #define PE_OUTPUT_DAT false // save rec+nois data into the output root file
162 #define PE_OUTPUT_MED false // save median into the output root file
163 #define PE_OUTPUT_P50 false // save percentile 50 into the output root file
164 #define PE_OUTPUT_P90 false // save percentile 90 into the output root file
165 #define PE_OUTPUT_AVR false // save averaged into the output root file
166 #define PE_OUTPUT_RMS false // save rms into the output root file
233 bool fft=
false,
TString pdir=
"",
double P=0.99, EColor col1=kRed, EColor col2=kBlue, EColor col3=(EColor)418);
375 std::vector<wavearray<double> >
vINJ;
376 std::vector<wavearray<double> >
vREC;
377 std::vector<wavearray<double> >
vWHT;
378 std::vector<wavearray<double> >
vPCA;
379 std::vector<wavearray<double> >
vDAT;
380 std::vector<wavearray<double> >
vNUL;
381 std::vector<wavearray<double> >
cINJ;
383 std::vector<wavearray<double> >
vAVR;
384 std::vector<wavearray<double> >
vRMS;
386 std::vector<wavearray<double> >
vMED;
387 std::vector<wavearray<double> >
vL50;
388 std::vector<wavearray<double> >
vU50;
389 std::vector<wavearray<double> >
vL90;
390 std::vector<wavearray<double> >
vU90;
392 std::vector<wavearray<double> >
fREC;
393 std::vector<wavearray<double> >
fINJ;
394 std::vector<wavearray<double> >
fAVR;
395 std::vector<wavearray<double> >
fRMS;
397 std::vector<wavearray<double> >
fMED;
398 std::vector<wavearray<double> >
fL50;
399 std::vector<wavearray<double> >
fU50;
400 std::vector<wavearray<double> >
fL90;
401 std::vector<wavearray<double> >
fU90;
446 if(gCWB2G->
istage!=
CWB_STAGE_FULL) {cout<<
"CWB_Plugin_PE - Error : PE can be executed only in CWB_STAGE_FULL mode" << endl;
exit(1);}
449 cout <<
"CWB_Plugin_PE Error : PE enable only with pattern>0 !!!" << endl;
459 TString cwb_inet_options=
TString(gSystem->Getenv(
"CWB_INET_OPTIONS"));
462 if(gOPT.ced_inj && gOPT.ced_rinj) gOPT.ced_rinj=
false;
466 cout<<
"CWB_Plugin_PE - Error : when simulation=4, PE multitask can be executed only with nfactors=1" << endl;
exit(1);
472 if(gOPT.output_inj || gOPT.output_rec ||
473 gOPT.output_med || gOPT.output_p50 ||
474 gOPT.output_p90 || gOPT.output_avr ||
475 gOPT.output_rms || gOPT.output_wht || gOPT.output_dat) cfg->
online=
true;
484 if(gtrial.CompareTo(
"")!=0) {
485 if(!gtrial.IsDigit()) {cout<<
"CWB_Plugin_PE - Error : when simulation=0, CWB_MDC_FACTOR must be an interger > 0" << endl;
exit(1);}
486 gMTRIAL=gtrial.Atoi();
487 if(gMTRIAL==0) {cout<<
"CWB_Plugin_PE - Error : when simulation=0, CWB_MDC_FACTOR must be an interger > 0" << endl;
exit(1);}
488 if(gMTRIAL!=gOPT.trials) {
491 cout <<
"CWB_Plugin_PE - MultiTask Mode : new data_label = " << cfg->
data_label << endl;
494 gOPT.output_inj=
false;
495 gOPT.output_rec=
false;
496 gOPT.output_wht=
false;
497 gOPT.output_dat=
false;
498 gOPT.output_med=
false;
499 gOPT.output_p50=
false;
500 gOPT.output_p90=
false;
501 gOPT.output_avr=
false;
502 gOPT.output_rms=
false;
506 bool last_trial = (gMTRIAL==gOPT.trials) ?
true :
false;
507 if(!last_trial) {gOPT.trials=1;gMTRIAL*=-1;}
513 gRandom->SetSeed(gOPT.seed);
519 if(gOPT.signal!=0 && gOPT.signal!=1 && gOPT.signal!=2) {
521 cout <<
"CWB_Plugin_PE Error : option pe_signal not allowed : must be (0/1/2) !!!" << endl;
528 cout <<
"CWB_Plugin_PE Error : option pe_signal=1 is allowed only in simulation mode !!!" << endl;
535 cout <<
"CWB_Plugin_PE Error : option pe_skymask=3 is allowed only in simulation mode !!!" << endl;
546 cout <<
"CWB_Plugin_PE - Error : gps time must be within the segment interval: "
559 cout <<
"CWB_Plugin_PE Error : detected " << NET->
mdcTime.size() <<
" injections."
560 <<
" Must be only one injection per segment !!!" << endl;
576 if(gOPT.multitask)
EXPORT(
int,gLRETRY,
"gLRETRY = 1;")
577 else EXPORT(
int,gLRETRY,TString::Format(
"gLRETRY = %d;",gOPT.trials).Data())
585 if(cid.
size()==0)
return;
588 if(gOPT.trials==0)
return;
591 int gLRETRY=-1;
IMPORT(
int,gLRETRY)
593 cout <<
"-------------------------------------------------------------------" << endl;
594 cout <<
"-----> CWB_Plugin_PE.C -> ID : " << ID
595 <<
" -> gLRETRY : " << gLRETRY <<
"/" << gOPT.trials << endl;
596 cout <<
"-------------------------------------------------------------------" << endl;
602 if(gLRETRY==gOPT.trials) {
611 if(gLRETRY!=gOPT.trials) {
614 if(gOPT.skymask==1) {
616 int isky = (
int)gHSKYPROB.GetRandom();
617 double th = 90-gSKYPROB.
getTheta(isky);
622 if(gOPT.skymask==2) {
626 if(gOPT.skymask==3) {
635 if(gOPT.signal==0) jtype=
'r';
636 if(gOPT.signal==1) jtype=
'j';
637 if(gOPT.signal==2) jtype=
'v';
642 if((gOPT.ced_dump>=0) && (gLRETRY==gOPT.ced_dump)) cfg->
cedDump=
true;
649 if((gLRETRY==0)&&(gMTRIAL==gOPT.trials)) {
674 int gLRETRY=-1;
IMPORT(
int,gLRETRY)
692 for(
int k=0;
k<
K;
k++) {
700 int ID = size_t(
id.data[
j]+0.5);
707 iSNR[
n].push_back(EVT->
iSNR[
n]);
708 oSNR[
n].push_back(EVT->
oSNR[
n]);
709 ioSNR[
n].push_back(EVT->
ioSNR[
n]);
711 for(
int n=0;
n<6;
n++) {
712 vCHIRP[
n].push_back(EVT->
chirp[
n]);
718 cout <<
"event parameters : ID -> " << ID << endl;
720 cout <<
"rec_theta : " << EVT->
theta[0] <<
" rec_phi : " << EVT->
phi[0] << endl;
721 cout <<
"SNRnet : " << sqrt(EVT->
likelihood) <<
" netcc[0] : " << EVT->
netcc[0]
722 <<
" rho[0] : " << EVT->
rho[0] <<
" size : " << EVT->
size[0] << endl;
733 if((gLRETRY==gOPT.trials) || (gOPT.multitask && gLRETRY==1)) {
737 if(vINJ.size()!=
nIFO) {
738 cout <<
"CWB_Plugin_PE Error : Injection Waveform Not Found !!!" << endl;
773 gSEGGPS = EVT->
gps[0];
774 gITHETA = EVT->
theta[1]; gIPHI = EVT->
phi[1];
775 gOTHETA = EVT->
theta[3]; gOPHI = EVT->
phi[3];
784 if(gLRETRY!=gOPT.trials) {
805 cout <<
"dirCED : " << dirCED << endl;
815 gStyle->SetOptFit(1111);
829 gStyle->SetOptFit(0000);
857 int gLRETRY=-1;
IMPORT(
int,gLRETRY)
861 cout <<
"-------------------------------------------------------------------" << endl;
862 cout <<
"-----> CWB_Plugin_PE.C -> RedoAnalysis : "
863 <<
" -> gLRETRY : " << gLRETRY <<
"/" << gOPT.trials << endl;
864 cout <<
"-------------------------------------------------------------------" << endl;
866 EXPORT(
int,gLRETRY,TString::Format(
"gLRETRY = %d;",--gLRETRY).Data())
867 }
while(csize==0 && gLRETRY>0);
868 EXPORT(
int,gLRETRY,TString::Format(
"gLRETRY = %d;",++gLRETRY).Data())
888 std::vector<netpixel> vPIX;
892 if(V>cfg->
BATCH)
return vPIX;
897 __m128* _xi = (__m128*) xi.
data;
898 __m128* _XI = (__m128*) XI.
data;
900 __m128* _aa = (__m128*) NET->
a_00.
data;
901 __m128* _AA = (__m128*) NET->
a_90.
data;
904 for(
int j=0;
j<V;
j++) {
909 if(ee>En) nPC++;
else ee=0.;
916 for(
int j=0;
j<V;
j++) {
918 vPIX.push_back(*pix);
936 if(V>cfg->
BATCH)
return;
937 for(
int j=0;
j<V;
j++) {
942 while(!vPIX->empty()) vPIX->pop_back();
943 vPIX->clear(); std::vector<netpixel>().swap(*vPIX);
948 std::vector<wavearray<double> > xPCA;
950 std::vector<netpixel> vPIX;
951 vPIX =
DoPCA(NET, cfg, lag,
id);
973 std::vector<netpixel> vPIX;
974 vPIX =
DoPCA(NET, cfg, lag,
id);
979 TString xtitle = TString::Format(
"Time (sec) : GPS OFFSET = %.3f",EVT->
gps[0]);
984 sprintf(fname,
"%s/pca_l_tfmap_scalogram.png",pdir.Data());
985 cout <<
"write " << fname << endl;
987 WTS.
hist2D->GetYaxis()->SetRangeUser(EVT->
low[0],EVT->
high[0]);
988 WTS.
hist2D->GetXaxis()->SetTitle(xtitle);
1005 std::vector<wavearray<double> > xINJ;
1009 double recTime = EVT->
time[0];
1014 double injTime = 1.e12;
1016 for(
int m=0;
m<
M;
m++) {
1019 if(T<injTime && INJ.
fill_in(NET,mdcID)) {
1034 pwfINJ[
n] = INJ.
pwf[
n];
1035 if (pwfINJ[
n]==NULL) {
1036 cout <<
"CWB_Plugin_Boostrap.C : Error : Injected waveform not saved !!! : detector "
1041 double rFactor = 1.;
1042 rFactor *= factor>0 ? factor : 1;
1045 xINJ.push_back(*wf);
1056 std::vector<wavearray<double> > xREC;
1064 int idSize = pd->
RWFID.size();
1066 for(
int mm=0; mm<idSize; mm++)
if (pd->
RWFID[mm]==ID) wfIndex=mm;
1069 cout <<
"CWB_Plugin_Boostrap.C : Error : Reconstructed waveform not saved !!! : ID -> "
1070 << ID <<
" : detector " << NET->
ifoName[
n] << endl;
1073 if(wfIndex>=0) pwfREC[
n] = pd->
RWFP[wfIndex];
1076 xREC.push_back(*wf);
1085 if(type!=
'S' && type!=
's' && type!=
'W' && type!=
'w' && type!=
'N' && type!=
'n') {
1086 cout <<
"GetWaveform : not valid type : Abort" << endl;
exit(1);
1089 std::vector<wavearray<double> > vWAV;
1094 if(type==
'S' || type==
's') {
1101 vWAV.push_back(*wf);
1105 if(type==
'W' || type==
'w') {
1112 vWAV.push_back(*wf);
1116 if(type==
'N' || type==
'n') {
1126 vWAV.push_back(*wf);
1136 std::vector<wavearray<double> > vSIG;
1149 int V4 = V + (V%4 ? 4 - V%4 : 0);
1161 std::vector<float*> _DAT;
1162 std::vector<float*> _vtd;
1163 std::vector<float*> _vTD;
1168 for(i=0; i<
NIFO; i++) {
1169 ptmp = (
float*)_mm_malloc((V4*3+8)*
sizeof(float),32);
1170 for(j=0; j<(V4*3+8); j++) ptmp[j]=0; _DAT.push_back(ptmp);
1171 ptmp = (
float*)_mm_malloc(tsize*V4*
sizeof(
float),32);
1172 for(j=0; j<tsize*V4; j++) ptmp[j]=0; _vtd.push_back(ptmp);
1173 ptmp = (
float*)_mm_malloc(tsize*V4*
sizeof(
float),32);
1174 for(j=0; j<tsize*V4; j++) ptmp[j]=0; _vTD.push_back(ptmp);
1177 std::vector<float*> _AVX;
1179 float* p_et = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1180 for(j=0; j<V4; j++) p_et[j]=0; _AVX.push_back(p_et);
1181 float* pMSK = (
float*)_mm_malloc(V44*
sizeof(
float),32);
1182 for(j=0; j<V44; j++) pMSK[j]=0; _AVX.push_back(pMSK); pMSK[V4]=
nIFO;
1183 float* p_fp = (
float*)_mm_malloc(V44*
sizeof(
float),32);
1184 for(j=0; j<V44; j++) p_fp[j]=0; _AVX.push_back(p_fp);
1185 float* p_fx = (
float*)_mm_malloc(V44*
sizeof(
float),32);
1186 for(j=0; j<V44; j++) p_fx[j]=0; _AVX.push_back(p_fx);
1187 float* p_si = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1188 for(j=0; j<V4; j++) p_si[j]=0; _AVX.push_back(p_si);
1189 float* p_co = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1190 for(j=0; j<V4; j++) p_co[j]=0; _AVX.push_back(p_co);
1191 float* p_uu = (
float*)_mm_malloc((V4+16)*
sizeof(float),32);
1192 for(j=0; j<V4+16; j++) p_uu[j]=0; _AVX.push_back(p_uu);
1193 float* p_UU = (
float*)_mm_malloc((V4+16)*
sizeof(float),32);
1194 for(j=0; j<V4+16; j++) p_UU[j]=0; _AVX.push_back(p_UU);
1195 float* p_vv = (
float*)_mm_malloc((V4+16)*
sizeof(float),32);
1196 for(j=0; j<V4+16; j++) p_vv[j]=0; _AVX.push_back(p_vv);
1197 float* p_VV = (
float*)_mm_malloc((V4+16)*
sizeof(float),32);
1198 for(j=0; j<V4+16; j++) p_VV[j]=0; _AVX.push_back(p_VV);
1199 float* p_au = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1200 for(j=0; j<V4; j++) p_au[j]=0; _AVX.push_back(p_au);
1201 float* p_AU = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1202 for(j=0; j<V4; j++) p_AU[j]=0; _AVX.push_back(p_AU);
1203 float* p_av = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1204 for(j=0; j<V4; j++) p_av[j]=0; _AVX.push_back(p_av);
1205 float* p_AV = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1206 for(j=0; j<V4; j++) p_AV[j]=0; _AVX.push_back(p_AV);
1207 float* p_uv = (
float*)_mm_malloc(V4*4*
sizeof(
float),32);
1208 for(j=0; j<V4*4; j++) p_uv[j]=0; _AVX.push_back(p_uv);
1209 float* p_ee = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1210 for(j=0; j<V4; j++) p_ee[j]=0; _AVX.push_back(p_ee);
1211 float* p_EE = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1212 for(j=0; j<V4; j++) p_EE[j]=0; _AVX.push_back(p_EE);
1213 float* pTMP=(
float*)_mm_malloc(V4*4*NIFO*
sizeof(
float),32);
1214 for(j=0; j<V4*4*
NIFO; j++) pTMP[j]=0; _AVX.push_back(pTMP);
1215 float* p_ni = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1216 for(j=0; j<V4; j++) p_ni[j]=0; _AVX.push_back(p_ni);
1217 float* p_ec = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1218 for(j=0; j<V4; j++) p_ec[j]=0; _AVX.push_back(p_ec);
1219 float* p_gn = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1220 for(j=0; j<V4; j++) p_gn[j]=0; _AVX.push_back(p_gn);
1221 float* p_ed = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1222 for(j=0; j<V4; j++) p_ed[j]=0; _AVX.push_back(p_ed);
1223 float* p_rn = (
float*)_mm_malloc(V4*
sizeof(
float),32);
1224 for(j=0; j<V4; j++) p_rn[j]=0; _AVX.push_back(p_rn);
1227 for(
int j=0; j<V; j++) {
1229 if(!pix->
core)
continue;
1242 for(
int i=0; i<
NIFO; i++) {
1260 for(
int j=0; j<V; j++) {
1262 pix->
likelihood = (pMSK[
j]>0) ? (p_ee[j]+p_EE[j])/2 : 0;
1263 if(pMSK[j]>0) pix->
core =
true;
1265 for(
int i=0; i<
nIFO; i++) {
1266 pix->
setdata(
double(pd[i][j]),
'S',i);
1267 pix->
setdata(
double(pD[i][j]),
'P',i);
1272 #ifdef SAVE_REDUCED_INJ
1273 for(
int i=0; i<
nIFO; i++) {
1276 sprintf(title,
"%s (TIME) : vREC(red) - vSIG(blue)",NET->
ifoName[i]);
1277 sprintf(ofname,
"%s_vREC_vSIG_time_id_%d.%s",NET->
ifoName[i],
id,
"root");
1278 PlotWaveform(ofname, title, cfg,&vREC[i], &vSIG[i], NULL, NULL,
false);
1279 sprintf(title,
"%s (TIME) : vINJ(red) - vSIG(blue)",NET->
ifoName[i]);
1280 sprintf(ofname,
"%s_vINJ_vSIG_time_id_%d.%s",NET->
ifoName[i],
id,
"root");
1281 PlotWaveform(ofname, title, cfg,&vINJ[i], &vSIG[i], NULL, NULL,
false);
1295 gPE.ff[0]=0; gPE.ff[1]=0;
1296 gPE.of[0]=0; gPE.of[1]=0;
1298 int size = iSNR[0].size();
1302 double isnr=0;
for(
int n=0;
n<
nIFO;
n++) isnr+= iSNR[
n][
i];
1303 double osnr=0;
for(
int n=0;
n<
nIFO;
n++) osnr+= oSNR[
n][
i];
1304 double iosnr=0;
for(
int n=0;
n<
nIFO;
n++) iosnr+=ioSNR[
n][
i];
1305 double ff = iosnr/sqrt(isnr*isnr);
1306 double of = iosnr/sqrt(isnr*osnr);
1309 gPE.ff[1] += pow(ff,2);
1312 gPE.of[1] += pow(of,2);
1316 gPE.ff[1] = sqrt(gPE.ff[1]/size-gPE.ff[0]*gPE.ff[0]);
1319 gPE.of[1] = sqrt(gPE.of[1]/size-gPE.of[0]*gPE.of[0]);
1327 int size = vCHIRP[1].size();
1331 gPE.mch[0] += vCHIRP[1][
i];
1332 gPE.mch[1] += pow(vCHIRP[1][
i],2);
1336 gPE.mch[1] = sqrt(gPE.mch[1]/size-gPE.mch[0]*gPE.mch[0]);
1344 int size = iSNR[0].size();
1349 for(
int n=0;
n<
nIFO;
n++) osnr+=oSNR[
n][
i];
1350 gPE.snet[0] += sqrt(osnr);
1351 gPE.snet[1] += osnr;
1354 gPE.snet[0] /=
size;
1355 gPE.snet[1] = sqrt(gPE.snet[1]/size-gPE.snet[0]*gPE.snet[0]);
1368 int sCuts = pwc->
sCuts[
id-1];
1369 pwc->
sCuts[
id-1] = 0;
1377 int tsize = pix->
tdAmp[0].size();
1378 if(!tsize || tsize&1) {
1379 cout<<
"GetNullPixels error: wrong pixel TD data\n";
1383 cout <<
"tsize :" << tsize << endl;
1384 int V4 = V + (V%4 ? 4 - V%4 : 0);
1385 std::vector<int>* vtof = &(pwc->
nTofF[
id-1]);
1387 for(
int n=0;
n<2*7*
nIFO;
n++) gPE.nstat[
n]=0;
1390 for(
int m=0;
m<nIFO;
m++)
cnt[
m]=0;
1396 for(
int j=0;
j<V;
j++) {
1402 double SS = a_90[j*
NIFO+
m];
1403 ee += pow(ss,2)+pow(SS,2);
1416 int l = (*vtof)[
m]+tsize/2;
1417 double dd = pix->
tdAmp[
m].data[
l];
1418 double DD = pix->
tdAmp[
m].data[l+tsize];
1424 double SS = a_90[j*
NIFO+
m];
1427 aNUL[
m].push_back(dd-ss);
1430 ANUL[
m].push_back(DD-SS);
1436 pwc->
sCuts[
id-1] = sCuts;
1457 bool fft,
TString pdir,
double P, EColor col1, EColor col2, EColor col3) {
1463 if(wref==NULL)
return;
1464 if(wf1==NULL)
return;
1465 else tmin=wf1->
start();
1466 if(wf2!=NULL)
if(wf2->
start()<tmin) tmin=wf2->
start();
1467 if(wf3!=NULL)
if(wf3->
start()<tmin) tmin=wf3->
start();
1472 if(wf2!=NULL)
if(wf2!=wf1) wf2->
start(wf2->
start()-tmin);
1473 if(wf3!=NULL)
if(wf3!=wf1 && wf3!=wf2) wf3->
start(wf3->
start()-tmin);
1474 if(wref!=wf1 && wref!=wf2 && wref!=wf3) wref->
start(wref->
start()-tmin);
1481 for(
int i=0;
i<wf1->
size();
i++) {
1483 if(time>bT && time<eT) wf[k++]=wf1->
data[
i];
1489 PTS.
graph[0]->GetXaxis()->SetRangeUser(bT, eT);
1491 PTS.
graph[0]->SetLineWidth(1);
1492 PTS.
graph[0]->SetTitle(title);
1494 TString xtitle = TString::Format(
"Time (sec) : GPS OFFSET = %.3f",tmin);
1495 PTS.
graph[0]->GetXaxis()->SetTitle(xtitle);
1501 for(
int i=0;
i<wf2->
size();
i++) {
1503 if(time>bT && time<eT) wf[k++]=wf2->
data[
i];
1511 if(wf2!=NULL) PTS.
graph[1]->SetLineWidth(2);
1518 for(
int i=0;
i<wf3->
size();
i++) {
1520 if(time>bT && time<eT) wf[k++]=wf3->
data[
i];
1528 if(wf3!=NULL) PTS.
graph[2]->SetLineWidth(1);
1532 if(wf2!=NULL)
if(wf2!=wf1) wf2->
start(wf2->
start()+tmin);
1533 if(wf3!=NULL)
if(wf3!=wf1 && wf3!=wf2) wf3->
start(wf3->
start()+tmin);
1534 if(wref!=wf1 && wref!=wf2 && wref!=wf3) wref->
start(wref->
start()+tmin);
1543 PTS.
canvas->Print(ofname);
1544 cout <<
"write : " << ofname << endl;
1550 double R=wf1->
rate();
1552 double b_wf1 = wf1->
start();
1553 double e_wf1 = wf1->
start()+wf1->
size()/R;
1554 double b_wf2 = wf2->
start();
1555 double e_wf2 = wf2->
start()+wf2->
size()/R;
1557 int o_wf1 = b_wf1>b_wf2 ? 0 :
int((b_wf2-b_wf1)*R+0.5);
1558 int o_wf2 = b_wf1<b_wf2 ? 0 :
int((b_wf1-b_wf2)*R+0.5);
1560 double startXCOR = b_wf1>b_wf2 ? b_wf1 : b_wf2;
1561 double endXCOR = e_wf1<e_wf2 ? e_wf1 : e_wf2;
1562 int sizeXCOR =
int((endXCOR-startXCOR)*R+0.5);
1567 for(
int i=0;
i<sizeXCOR;
i++) wf12 += wf1->
data[
i+o_wf1] * wf2->
data[
i+o_wf2];
1568 for(
int i=0;i<wf1->
size();
i++) wf11 += pow(wf1->
data[
i],2);
1569 for(
int i=0;
i<wf2->
size();
i++) wf22 += pow(wf2->
data[
i],2);
1571 double FF = wf12/sqrt(wf11*wf22);
1578 double R=wf1->
rate();
1580 double b_wf1 = wf1->
start();
1581 double e_wf1 = wf1->
start()+wf1->
size()/R;
1582 double b_wf2 = wf2->
start();
1583 double e_wf2 = wf2->
start()+wf2->
size()/R;
1585 int o_wf1 = b_wf1>b_wf2 ? 0 :
int((b_wf2-b_wf1)*R+0.5);
1586 int o_wf2 = b_wf1<b_wf2 ? 0 :
int((b_wf1-b_wf2)*R+0.5);
1588 double startXCOR = b_wf1>b_wf2 ? b_wf1 : b_wf2;
1589 double endXCOR = e_wf1<e_wf2 ? e_wf1 : e_wf2;
1590 int sizeXCOR =
int((endXCOR-startXCOR)*R+0.5);
1595 wfdif.
start(b_wf1+
double(o_wf1)/R);
1597 for(
int i=0;
i<sizeXCOR;
i++) wfdif[
i] = wf1->
data[
i+o_wf1] - wf2->
data[
i+o_wf2];
1616 sprintf(title,
"%s (TIME) : INJ(red) - SM(blue)",NET->
ifoName[ifoID]);
1618 PlotWaveform(ofname, title, cfg, wf, &ss, NULL, NULL,
false);
1624 int n_amp_cal_err=0;
1625 int n_phs_cal_err=0;
1626 if(options.CompareTo(
"")!=0) {
1627 cout << options << endl;
1628 if(!options.Contains(
"--")) {
1631 for(
int j=0;
j<token->GetEntries();
j++){
1633 TObjString* tok = (TObjString*)token->At(
j);
1634 TString stok = tok->GetString();
1636 if(stok.Contains(
"pe_trials=") || stok.Contains(
"pe_retry=")) {
1638 pe_trials.Remove(0,pe_trials.Last(
'=')+1);
1639 if(pe_trials.IsDigit()) gOPT.trials=pe_trials.Atoi();
1641 cout <<
"ReadUserOptions Error : trials must be <= MAX_TRIALS : " <<
MAX_TRIALS << endl;
exit(1);
1645 if(stok.Contains(
"pe_ced_dump=")) {
1647 pe_ced_dump.Remove(0,pe_ced_dump.Last(
'=')+1);
1648 if(pe_ced_dump.IsDigit()) gOPT.ced_dump=pe_ced_dump.Atoi();
1651 if(stok.Contains(
"pe_ced_enable=")) {
1653 pe_ced_enable.Remove(0,pe_ced_enable.Last(
'=')+1);
1654 if(pe_ced_enable==
"tfmap") gOPT.ced_tfmap =
true;
1655 if(pe_ced_enable==
"rdr") gOPT.ced_rdr =
true;
1656 if(pe_ced_enable==
"psd") gOPT.ced_psd =
true;
1657 if(pe_ced_enable==
"skymap") gOPT.ced_skymap =
true;
1658 if(pe_ced_enable==
"rec") gOPT.ced_rec =
true;
1659 if(pe_ced_enable==
"inj") gOPT.ced_inj =
true;
1660 if(pe_ced_enable==
"rinj") gOPT.ced_rinj =
true;
1661 if(pe_ced_enable==
"cm") gOPT.ced_cm =
true;
1662 if(pe_ced_enable==
"distr") gOPT.ced_distr =
true;
1663 if(pe_ced_enable==
"null") gOPT.ced_null =
true;
1664 if(pe_ced_enable==
"pca") gOPT.ced_pca =
true;
1667 if(stok.Contains(
"pe_ced_disable=")) {
1669 pe_ced_disable.Remove(0,pe_ced_disable.Last(
'=')+1);
1670 if(pe_ced_disable==
"tfmap") gOPT.ced_tfmap =
false;
1671 if(pe_ced_disable==
"rdr") gOPT.ced_rdr =
false;
1672 if(pe_ced_disable==
"psd") gOPT.ced_psd =
false;
1673 if(pe_ced_disable==
"skymap") gOPT.ced_skymap =
false;
1674 if(pe_ced_disable==
"rec") gOPT.ced_rec =
false;
1675 if(pe_ced_disable==
"inj") gOPT.ced_inj =
false;
1676 if(pe_ced_disable==
"rinj") gOPT.ced_rinj =
false;
1677 if(pe_ced_disable==
"cm") gOPT.ced_cm =
false;
1678 if(pe_ced_disable==
"distr") gOPT.ced_distr =
false;
1679 if(pe_ced_disable==
"null") gOPT.ced_null =
false;
1680 if(pe_ced_disable==
"pca") gOPT.ced_pca =
false;
1683 if(stok.Contains(
"pe_seed=")) {
1685 pe_seed.Remove(0,pe_seed.Last(
'=')+1);
1686 if(pe_seed.IsDigit()) gOPT.seed=pe_seed.Atoi();
1689 if(stok.Contains(
"pe_gps=")) {
1691 pe_gps.Remove(0,pe_gps.Last(
'=')+1);
1692 if(pe_gps.IsFloat()) gOPT.gps=pe_gps.Atof();
1695 if(stok.Contains(
"pe_amp_cal_err=")) {
1697 pe_amp_cal_err.Remove(0,pe_amp_cal_err.Last(
'=')+1);
1698 if(pe_amp_cal_err.IsFloat()) gOPT.amp_cal_err[n_amp_cal_err]=pe_amp_cal_err.Atof();
1699 if(n_amp_cal_err<(
NIFO_MAX-1)) n_amp_cal_err++;
1702 if(stok.Contains(
"pe_phs_cal_err=")) {
1704 pe_phs_cal_err.Remove(0,pe_phs_cal_err.Last(
'=')+1);
1705 if(pe_phs_cal_err.IsFloat()) gOPT.phs_cal_err[n_phs_cal_err]=pe_phs_cal_err.Atof();
1706 if(n_phs_cal_err<(
NIFO_MAX-1)) n_phs_cal_err++;
1709 if(stok.Contains(
"pe_skymask=")) {
1711 pe_skymask.Remove(0,pe_skymask.Last(
'=')+1);
1712 if(pe_skymask==
"false") gOPT.skymask=0;
1713 if(pe_skymask==
"true") gOPT.skymask=1;
1714 if(pe_skymask.IsDigit()) gOPT.skymask=pe_skymask.Atoi();
1717 if(stok.Contains(
"pe_signal=")) {
1719 pe_signal.Remove(0,pe_signal.Last(
'=')+1);
1720 if(pe_signal.IsDigit()) gOPT.signal=pe_signal.Atoi();
1723 if(stok.Contains(
"pe_noise=")) {
1725 pe_noise.Remove(0,pe_noise.Last(
'=')+1);
1726 if(pe_noise==
"false") gOPT.noise=0;
1727 if(pe_noise==
"true") gOPT.noise=1;
1728 if(pe_noise.IsDigit()) gOPT.noise=pe_noise.Atoi();
1731 if(stok.Contains(
"pe_multitask=")) {
1733 pe_multitask.Remove(0,pe_multitask.Last(
'=')+1);
1734 if(pe_multitask==
"true") gOPT.multitask=
true;
1735 if(pe_multitask==
"false") gOPT.multitask=
false;
1738 if(stok.Contains(
"pe_output_enable=")) {
1739 TString pe_output_enable=stok;
1740 pe_output_enable.Remove(0,pe_output_enable.Last(
'=')+1);
1741 if(pe_output_enable==
"inj") gOPT.output_inj=
true;
1742 if(pe_output_enable==
"rec") gOPT.output_rec=
true;
1743 if(pe_output_enable==
"wht") gOPT.output_wht=
true;
1744 if(pe_output_enable==
"dat") gOPT.output_dat=
true;
1745 if(pe_output_enable==
"med") gOPT.output_med=
true;
1746 if(pe_output_enable==
"p50") gOPT.output_p50=
true;
1747 if(pe_output_enable==
"p90") gOPT.output_p90=
true;
1748 if(pe_output_enable==
"avr") gOPT.output_avr=
true;
1749 if(pe_output_enable==
"rms") gOPT.output_rms=
true;
1752 if(stok.Contains(
"pe_output_disable=")) {
1753 TString pe_output_disable=stok;
1754 pe_output_disable.Remove(0,pe_output_disable.Last(
'=')+1);
1755 if(pe_output_disable==
"inj") gOPT.output_inj=
false;
1756 if(pe_output_disable==
"rec") gOPT.output_rec=
false;
1757 if(pe_output_disable==
"wht") gOPT.output_wht=
false;
1758 if(pe_output_disable==
"dat") gOPT.output_dat=
false;
1759 if(pe_output_disable==
"med") gOPT.output_med=
false;
1760 if(pe_output_disable==
"p50") gOPT.output_p50=
false;
1761 if(pe_output_disable==
"p90") gOPT.output_p90=
false;
1762 if(pe_output_disable==
"avr") gOPT.output_avr=
false;
1763 if(pe_output_disable==
"rms") gOPT.output_rms=
false;
1773 cout <<
"-----------------------------------------" << endl;
1774 cout <<
"PE config options " << endl;
1775 cout <<
"-----------------------------------------" << endl << endl;
1776 cout <<
"PE_TRIALS " << gOPT.trials << endl;
1777 cout <<
"PE_SIGNAL " << gOPT.signal << endl;
1778 cout <<
"PE_NOISE " << gOPT.noise << endl;
1779 for(
int n=0;
n<cfg->
nIFO;
n++) {
1780 cout <<
"PE_AMP_CAL_ERR " << cfg->
ifo[
n] <<
" " << gOPT.amp_cal_err[
n] << endl;
1781 cout <<
"PE_PHS_CAL_ERR " << cfg->
ifo[
n] <<
" " << gOPT.phs_cal_err[
n] << endl;
1783 cout <<
"PE_SKYMASK " << gOPT.skymask << endl;
1784 cout <<
"PE_SEED " << gOPT.seed << endl;
1786 cout <<
"PE_GPS " << gOPT.gps << endl;
1789 cout <<
"PE_MULTITASK " << gOPT.multitask << endl;
1791 cout <<
"PE_CED_DUMP " << gOPT.ced_dump << endl;
1793 cout <<
"PE_CED_TFMAP " << gOPT.ced_tfmap << endl;
1794 cout <<
"PE_CED_RDR " << gOPT.ced_rdr << endl;
1795 cout <<
"PE_CED_PSD " << gOPT.ced_psd << endl;
1796 cout <<
"PE_CED_SKYMAP " << gOPT.ced_skymap << endl;
1797 cout <<
"PE_CED_REC " << gOPT.ced_rec << endl;
1798 cout <<
"PE_CED_INJ " << gOPT.ced_inj << endl;
1799 cout <<
"PE_CED_rINJ " << gOPT.ced_rinj << endl;
1800 cout <<
"PE_CED_CM " << gOPT.ced_cm << endl;
1801 cout <<
"PE_CED_DISTR " << gOPT.ced_distr << endl;
1802 cout <<
"PE_CED_NULL " << gOPT.ced_null << endl;
1803 cout <<
"PE_CED_PCA " << gOPT.ced_pca << endl;
1805 cout <<
"PE_OUTPUT_INJ " << gOPT.output_inj << endl;
1806 cout <<
"PE_OUTPUT_REC " << gOPT.output_rec << endl;
1807 cout <<
"PE_OUTPUT_WHT " << gOPT.output_wht << endl;
1808 cout <<
"PE_OUTPUT_DAT " << gOPT.output_dat << endl;
1809 cout <<
"PE_OUTPUT_MED " << gOPT.output_med << endl;
1810 cout <<
"PE_OUTPUT_P50 " << gOPT.output_p50 << endl;
1811 cout <<
"PE_OUTPUT_P90 " << gOPT.output_p90 << endl;
1812 cout <<
"PE_OUTPUT_AVR " << gOPT.output_avr << endl;
1813 cout <<
"PE_OUTPUT_RMS " << gOPT.output_rms << endl;
1824 if(!out.good()) {cout <<
"DumpUserOptions : Error Opening File : " << ofName << endl;
exit(1);}
1833 out <<
"--------------------------------" << endl;
1834 out <<
"PE version " << endl;
1835 out <<
"--------------------------------" << endl;
1837 out << version << endl;
1840 out <<
"--------------------------------" << endl;
1841 out <<
"PE config options " << endl;
1842 out <<
"--------------------------------" << endl;
1845 out <<
"pe_trials " << gOPT.trials << endl;
1846 info =
"// number of trials (def=100)";
1847 out << tabs << info << endl;
1849 out <<
"pe_signal " << gOPT.signal << endl;
1850 info =
"// signal used for trials : 0(def)/1/2";
1851 out << tabs << info << endl;
1852 info =
"// 0 - reconstructed waveform";
1853 out << tabs << info << endl;
1854 info =
"// 1 - injected waveform (available only in simulation mode)";
1855 out << tabs << info << endl;
1856 info =
"// 2 - original waveform = (reconstructed+null)";
1857 out << tabs << info << endl;
1859 out <<
"pe_noise " << gOPT.noise << endl;
1860 info =
"// noise used for trials : 0/1(def)/2 ";
1861 out << tabs << info << endl;
1862 info =
"// 0 - waveform is injected in the whitened HoT and apply Coherent+SuperCluster+Likelihood stages";
1863 out << tabs << info << endl;
1864 info =
"// 1 - add gaussian noise to likelihood sparse maps and apply Likelihood stage";
1865 out << tabs << info << endl;
1866 info =
"// 2 - add whitened HoT to likelihood sparse maps and apply Likelihood stage";
1867 out << tabs << info << endl;
1869 for(
int n=0;
n<cfg->
nIFO;
n++) {
1870 out <<
"pe_amp_cal_err " << cfg->
ifo[
n] <<
" " << gOPT.amp_cal_err[
n] << endl;
1872 info =
"// max percentage of amplitude miscalibration : def(0) -> disabled";
1873 out << tabs << info << endl;
1874 info =
"// if>0 -> uniform in [1-pe_amp_cal_err,1+pe_amp_cal_err] else gaus(1,pe_amp_cal_err)";
1875 out << tabs << info << endl;
1877 for(
int n=0;
n<cfg->
nIFO;
n++) {
1878 out <<
"pe_phs_cal_err " << cfg->
ifo[
n] <<
" " << gOPT.phs_cal_err[
n] << endl;
1880 info =
"// max phase (degrees) miscalibration : def(0) -> disabled";
1881 out << tabs << info << endl;
1882 info =
"// if>0 -> uniform in [-pe_phs_cal_err,+pe_phs_cal_err] else gaus(pe_phs_cal_err)";
1883 out << tabs << info << endl;
1885 out <<
"pe_multitask " << (gOPT.multitask ?
"enable" :
"disable") << endl;
1886 info =
"// if enabled, PE trials are executed by different jobs in multitask : only simulation=0 and single event (gps>0)";
1887 out << tabs << info << endl;
1889 out <<
"pe_ced_dump " << gOPT.ced_dump << endl;
1890 info =
"// dump CED at gLRETRY=PE_CED_DUMP (def(-1) -> disabled)";
1891 out << tabs << info << endl;
1893 out <<
"pe_skymask " << gOPT.skymask << endl;
1894 info =
"// skymask used for trials : 0(def)/1/2/3 ";
1895 out << tabs << info << endl;
1896 info =
"// 0 - disable -> search over all sky";
1897 out << tabs << info << endl;
1898 out << tabs <<
"// 1 - skymask with radius " <<
SKYMASK_RADIUS <<
"(deg) and source selected according the reconstructed skymap probability " << endl;
1899 info =
"// 2 - skymask select only the sky position where the waveform is reconstructed (the maximum of detection statistic)";
1900 out << tabs << info << endl;
1901 info =
"// 3 - skymask select only the sky position where the waveform has been injected";
1902 out << tabs << info << endl;
1904 out <<
"pe_seed " << gOPT.seed << endl;
1905 info =
"// seed used by PE for random generation - 0(def) -> random seed";
1906 out << tabs << info << endl;
1909 out <<
"pe_gps " << gOPT.gps << endl;
1910 info =
"// if >0 only gps +/- iwindow is analyzed - 0(def) -> disabled";
1911 out << tabs << info << endl;
1915 out <<
"pe_ced tfmap " << (gOPT.ced_tfmap ?
"enable" :
"disable") << endl;
1916 info =
"// tfmap - pe_ced_enable(def)/disable : Shows/Hides the Time-Frequency Maps Tab -> Spectrograms, Scalograms, TF Likelihood,Null";
1917 out << tabs << info << endl;
1918 out <<
"pe_ced rdr " << (gOPT.ced_rdr ?
"enable" :
"disable") << endl;
1919 info =
"// rdr - pe_ced_enable(def)/disable : Shows/Hides Reconstructed Detector Response Tab";
1920 out << tabs << info << endl;
1921 out <<
"pe_ced psd " << (gOPT.ced_psd ?
"enable" :
"disable") << endl;
1922 info =
"// psd - pe_ced_enable(def)/disable : Shows/Hides Power Spectral Density Tab";
1923 out << tabs << info << endl;
1924 out <<
"pe_ced skymap " << (gOPT.ced_skymap ?
"enable" :
"disable") << endl;
1925 info =
"// skymap - pe_ced_enable(def)/disable : Shows/Hides SkyMaps Tab";
1926 out << tabs << info << endl;
1927 out <<
"pe_ced rec " << (gOPT.ced_rec ?
"enable" :
"disable") << endl;
1928 info =
"// rec - pe_ced_enable(def)/disable : Shows/Hides Reconstructed Waveforms/Instantaneous-Frequency with errors Tab";
1929 out << tabs << info << endl;
1930 out <<
"pe_ced inj " << (gOPT.ced_inj ?
"enable" :
"disable") << endl;
1931 info =
"// inj - pe_ced_enable(def)/disable : Showsi/Hides Injection' tab reports the comparison with the injected whitened waveform";
1932 out << tabs << info << endl;
1933 out <<
"pe_ced rinj " << (gOPT.ced_rinj ?
"enable" :
"disable") << endl;
1934 info =
"// rinj - pe_ced_enable/disable(def) : Shows/Hides Injection' tab reports the comparison with the injected whitened waveform";
1935 out << tabs << info << endl;
1936 info =
" in time-frequency subspace of the PointEstimated waveform";
1937 out << tabs << info << endl;
1938 out <<
"pe_ced cm " << (gOPT.ced_cm ?
"enable" :
"disable") << endl;
1939 info =
"// cm - pe_ced_enable(def)/disable : Shows/Hides Chirp Mass Value/Error Distributions Tab";
1940 out << tabs << info << endl;
1941 out <<
"pe_ced distr " << (gOPT.ced_distr ?
"enable" :
"disable") << endl;
1942 info =
"// distr - pe_ced_enable/disable(def) : Shows/Hides Residuals Distributions Tab";
1943 out << tabs << info << endl;
1944 out <<
"pe_ced null " << (gOPT.ced_null ?
"enable" :
"disable") << endl;
1945 info =
"// null - pe_ced_enable/disable(def) : Shows/Hides Null Pixels Distributions Tab";
1946 out << tabs << info << endl;
1947 out <<
"pe_ced pca " << (gOPT.ced_pca ?
"enable" :
"disable") << endl;
1948 info =
"// pca - pe_ced_enable/disable(def) : Shows/Hides PCA likelihood TF plot in the Time-Frequency Maps Tab";
1949 out << tabs << info << endl;
1995 float* gSLAGSHIFT=NULL;
IMPORT(
float*,gSLAGSHIFT)
2000 TFile*
froot = NULL;
2001 TList *
files = (TList*)gROOT->GetListOfFiles();
2008 while ((file=(TSystemFile*)
next())) {
2009 fname = file->GetName();
2011 if(fname.Contains(
"wave_")) {
2012 froot=(TFile*)file;froot->cd();
2014 gOUTPUT.ReplaceAll(
".root.tmp",
".txt");
2019 cout <<
"CWB_Plugin_Boostrap.C : Error - output root file not found" << endl;
2023 cout <<
"CWB_Plugin_Boostrap.C : Error - output root file not found" << endl;
2029 if(gOPT.output_rec)
for(
int n=0;
n<
nIFO;
n++) gPE.wREC[
n] = &vREC[
n];
2030 if(gOPT.output_wht)
for(
int n=0;
n<
nIFO;
n++) gPE.wWHT[
n] = &vWHT[
n];
2031 if(gOPT.output_dat)
for(
int n=0;
n<
nIFO;
n++) gPE.wDAT[
n] = &vDAT[
n];
2032 if(gOPT.output_med)
for(
int n=0;
n<
nIFO;
n++) gPE.wMED[
n] = &vMED[
n];
2033 if(gOPT.output_p50)
for(
int n=0;
n<
nIFO;
n++) gPE.wL50[
n] = &vL50[
n];
2034 if(gOPT.output_p50)
for(
int n=0;
n<
nIFO;
n++) gPE.wU50[
n] = &vU50[
n];
2035 if(gOPT.output_p90)
for(
int n=0;
n<
nIFO;
n++) gPE.wL90[
n] = &vL90[
n];
2036 if(gOPT.output_p90)
for(
int n=0;
n<
nIFO;
n++) gPE.wU90[
n] = &vU90[
n];
2037 if(gOPT.output_avr)
for(
int n=0;
n<
nIFO;
n++) gPE.wAVR[
n] = &vAVR[
n];
2038 if(gOPT.output_rms)
for(
int n=0;
n<
nIFO;
n++) gPE.wRMS[
n] = &vRMS[
n];
2041 gTREE = (TTree *) froot->Get(
"waveburst");
2045 gTREE->Branch(
"pe_trials",&gPE.trials,
"pe_trials/I");
2046 gTREE->Branch(
"pe_erR",gPE.erR,TString::Format(
"pe_erR[%i]/F",11));
2047 gTREE->Branch(
"pe_erF",gPE.erF,TString::Format(
"pe_erF[%i]/F",11));
2048 gTREE->Branch(
"pe_nstat",gPE.nstat,TString::Format(
"pe_nstat[%i]/F",2*7*cfg->
nIFO));
2049 gTREE->Branch(
"pe_snet",gPE.snet,TString::Format(
"pe_snet[%i]/F",2));
2050 gTREE->Branch(
"pe_ff",gPE.ff,TString::Format(
"pe_ff[%i]/F",2));
2051 gTREE->Branch(
"pe_of",gPE.of,TString::Format(
"pe_of[%i]/F",2));
2052 gTREE->Branch(
"pe_mch",gPE.mch,TString::Format(
"pe_mch[%i]/F",2));
2055 if(gOPT.output_inj)
if(cfg->
simulation) gTREE->Branch(TString::Format(
"pe_wINJ_%d",
n).Data(),
"wavearray<double>",&gPE.wINJ[
n],32000,0);
2056 if(gOPT.output_rec) gTREE->Branch(TString::Format(
"pe_wREC_%d",
n).Data(),
"wavearray<double>",&gPE.wREC[
n],32000,0);
2057 if(gOPT.output_wht) gTREE->Branch(TString::Format(
"pe_wWHT_%d",
n).Data(),
"wavearray<double>",&gPE.wWHT[
n],32000,0);
2058 if(gOPT.output_dat) gTREE->Branch(TString::Format(
"pe_wDAT_%d",
n).Data(),
"wavearray<double>",&gPE.wDAT[
n],32000,0);
2059 if(gOPT.output_med) gTREE->Branch(TString::Format(
"pe_wMED_%d",
n).Data(),
"wavearray<double>",&gPE.wMED[
n],32000,0);
2060 if(gOPT.output_p50) gTREE->Branch(TString::Format(
"pe_wL50_%d",
n).Data(),
"wavearray<double>",&gPE.wL50[
n],32000,0);
2061 if(gOPT.output_p50) gTREE->Branch(TString::Format(
"pe_wU50_%d",
n).Data(),
"wavearray<double>",&gPE.wU50[
n],32000,0);
2062 if(gOPT.output_p90) gTREE->Branch(TString::Format(
"pe_wL90_%d",
n).Data(),
"wavearray<double>",&gPE.wL90[
n],32000,0);
2063 if(gOPT.output_p90) gTREE->Branch(TString::Format(
"pe_wU90_%d",
n).Data(),
"wavearray<double>",&gPE.wU90[
n],32000,0);
2064 if(gOPT.output_avr) gTREE->Branch(TString::Format(
"pe_wAVR_%d",
n).Data(),
"wavearray<double>",&gPE.wAVR[
n],32000,0);
2065 if(gOPT.output_rms) gTREE->Branch(TString::Format(
"pe_wRMS_%d",
n).Data(),
"wavearray<double>",&gPE.wRMS[
n],32000,0);
2068 if(gOPT.multitask&&(gMTRIAL!=gOPT.trials)) {
2069 for(
int n=0;
n<
nIFO;
n++) gTREE->Branch(TString::Format(
"mtpe_wREC_%d",
n).Data(),
"wavearray<double>",&wREC[0][
n],32000,0);
2070 gTREE->Branch(
"mtpe_skyprob",
"skymap",&wSKYPROB[0], 32000,0);
2083 if(cfg->
dump) EVT->
dopen(gOUTPUT.Data(),
const_cast<char*
>(
"a"),
false);
2084 EVT->
output2G(gTREE,NET,ID,k,factor);
2091 while(!vINJ.empty()) vINJ.pop_back();
2092 vINJ.clear(); std::vector<wavearray<double> >().swap(vINJ);
2093 while(!vREC.empty()) vREC.pop_back();
2094 vREC.clear(); std::vector<wavearray<double> >().swap(vREC);
2095 while(!vWHT.empty()) vWHT.pop_back();
2096 vWHT.clear(); std::vector<wavearray<double> >().swap(vWHT);
2097 while(!vPCA.empty()) vPCA.pop_back();
2098 vPCA.clear(); std::vector<wavearray<double> >().swap(vPCA);
2099 while(!vDAT.empty()) vDAT.pop_back();
2100 vDAT.clear(); std::vector<wavearray<double> >().swap(vDAT);
2101 while(!vNUL.empty()) vNUL.pop_back();
2102 vNUL.clear(); std::vector<wavearray<double> >().swap(vNUL);
2103 while(!cINJ.empty()) cINJ.pop_back();
2104 cINJ.clear(); std::vector<wavearray<double> >().swap(cINJ);
2106 while(!vAVR.empty()) vAVR.pop_back();
2107 vAVR.clear(); std::vector<wavearray<double> >().swap(vAVR);
2108 while(!vRMS.empty()) vRMS.pop_back();
2109 vRMS.clear(); std::vector<wavearray<double> >().swap(vRMS);
2111 while(!vMED.empty()) vMED.pop_back();
2112 vMED.clear(); std::vector<wavearray<double> >().swap(vMED);
2113 while(!vL50.empty()) vL50.pop_back();
2114 vL50.clear(); std::vector<wavearray<double> >().swap(vL50);
2115 while(!vU50.empty()) vU50.pop_back();
2116 vU50.clear(); std::vector<wavearray<double> >().swap(vU50);
2117 while(!vL90.empty()) vL90.pop_back();
2118 vL90.clear(); std::vector<wavearray<double> >().swap(vL90);
2119 while(!vU90.empty()) vU90.pop_back();
2120 vU90.clear(); std::vector<wavearray<double> >().swap(vU90);
2122 while(!fREC.empty()) fREC.pop_back();
2123 fREC.clear(); std::vector<wavearray<double> >().swap(fREC);
2124 while(!fINJ.empty()) fINJ.pop_back();
2125 fINJ.clear(); std::vector<wavearray<double> >().swap(fINJ);
2126 while(!fAVR.empty()) fAVR.pop_back();
2127 fAVR.clear(); std::vector<wavearray<double> >().swap(fAVR);
2128 while(!fRMS.empty()) fRMS.pop_back();
2129 fRMS.clear(); std::vector<wavearray<double> >().swap(fRMS);
2131 while(!fMED.empty()) fMED.pop_back();
2132 fMED.clear(); std::vector<wavearray<double> >().swap(fMED);
2133 while(!fL50.empty()) fL50.pop_back();
2134 fL50.clear(); std::vector<wavearray<double> >().swap(fL50);
2135 while(!fU50.empty()) fU50.pop_back();
2136 fU50.clear(); std::vector<wavearray<double> >().swap(fU50);
2137 while(!fL90.empty()) fL90.pop_back();
2138 fL90.clear(); std::vector<wavearray<double> >().swap(fL90);
2139 while(!fU90.empty()) fU90.pop_back();
2140 fU90.clear(); std::vector<wavearray<double> >().swap(fU90);
2142 while(!vRES.empty()) vRES.pop_back();
2143 vRES.clear(); std::vector<double >().swap(vRES);
2144 while(!fRES.empty()) fRES.pop_back();
2145 fRES.clear(); std::vector<double >().swap(fRES);
2148 while(!wREC[
n].empty()) wREC[
n].pop_back();
2149 wREC[
n].clear(); std::vector<wavearray<double> >().swap(wREC[
n]);
2153 while(!iSNR[
n].empty()) iSNR[
n].pop_back();
2154 iSNR[
n].clear(); std::vector<double>().swap(iSNR[
n]);
2155 while(!oSNR[n].empty()) oSNR[
n].pop_back();
2156 oSNR[
n].clear(); std::vector<double>().swap(oSNR[n]);
2157 while(!ioSNR[n].empty()) ioSNR[
n].pop_back();
2158 ioSNR[
n].clear(); std::vector<double>().swap(ioSNR[n]);
2161 for(
int n=0;
n<6;
n++) {
2162 while(!vCHIRP[
n].empty()) vCHIRP[
n].pop_back();
2163 vCHIRP[
n].clear(); std::vector<double>().swap(vCHIRP[
n]);
2166 while(!vLIKELIHOOD.empty()) vLIKELIHOOD.pop_back();
2167 vLIKELIHOOD.clear(); std::vector<double >().swap(vLIKELIHOOD);
2170 while(!aNUL[
n].empty()) aNUL[
n].pop_back();
2171 aNUL[
n].clear(); std::vector<double>().swap(aNUL[
n]);
2172 while(!ANUL[n].empty()) ANUL[
n].pop_back();
2173 ANUL[
n].clear(); std::vector<double>().swap(ANUL[n]);
2177 while(!aNSE[
n].empty()) aNSE[
n].pop_back();
2178 aNSE[
n].clear(); std::vector<double>().swap(aNSE[
n]);
2179 while(!ANSE[n].empty()) ANSE[
n].pop_back();
2180 ANSE[
n].clear(); std::vector<double>().swap(ANSE[n]);
2184 while(!vSS[
n].empty()) vSS[
n].pop_back();
2185 vSS[
n].clear(); std::vector<SSeries<double> >().swap(vSS[
n]);
2186 while(!rSS[n].empty()) rSS[
n].pop_back();
2187 rSS[
n].clear(); std::vector<SSeries<double> >().swap(rSS[n]);
2188 while(!jSS[n].empty()) jSS[
n].pop_back();
2189 jSS[
n].clear(); std::vector<SSeries<double> >().swap(jSS[n]);
2190 while(!dSS[n].empty()) dSS[
n].pop_back();
2191 dSS[
n].clear(); std::vector<SSeries<double> >().swap(dSS[n]);
2201 char title[256];
char ofname[256];
2227 #ifdef SET_WAVEFORM_CUTS
2230 sprintf(title,
"%s (TIME) : cINJ(red) - cINJ-vAVR(blue) - vRMS(green)",NET->
ifoName[n]);
2232 PlotWaveform(ofname, title, cfg, &cINJ[n], &wDIF, &vRMS[n], NULL,
false, pdir);
2241 #ifdef SET_WAVEFORM_CUTS
2243 sprintf(title,
"%s (TIME) : cINJ(red) - vREC(blue) - vAVR(green)",NET->
ifoName[n]);
2245 PlotWaveform(ofname, title, cfg, &cINJ[n], &vREC[n], &vAVR[n], NULL,
false, pdir);
2248 #ifdef SET_WAVEFORM_CUTS
2250 sprintf(title,
"%s (TIME) : cINJ(red) - cINJ(blue)",NET->
ifoName[n]);
2252 PlotWaveform(ofname, title, cfg, &cINJ[n], &cINJ[n], NULL, NULL,
false, pdir);
2259 double scale = sqrt(1<<cfg->
levelR);
2263 PlotWaveform(ofname,
"", cfg, &vWHT[n], &vREC[n], NULL, &vREC[n],
false, pdir,0.999,(EColor)16,kRed);
2265 PlotWaveform(ofname,
"", cfg, &vWHT[n], &vREC[n], NULL, &vREC[n],
true, pdir,0.999,(EColor)16,kRed);
2270 for(
int j=0;
j<2;
j++) {
2273 double P =
j==0 ? 0.995 : 0.9;
2281 sprintf(title,
"%s (TIME) : vREC (red) - vRMS:2vRMS (grey:light_gray)",NET->
ifoName[n]);
2284 PlotWaveformAsymmErrors(ofname,
"", cfg, &vREC[n], &vMED[n], &vL50[n], &vU50[n], &vL90[n], &vU90[n], &vREC[n], pdir, P);
2286 PlotWaveformErrors(ofname,
"", cfg, &vREC[n], &vAVR[n], &vRMS[n], &vREC[n], pdir, P);
2295 sprintf(title,
"%s (TIME) : vREC-vAVR(red) - 2*vRMS(gray)",NET->
ifoName[n]);
2298 PlotWaveformAsymmErrors(ofname,
"", cfg, &wDIF, &wDIF, &vL50[n], &vU50[n], &vL90[n], &vU90[n], &vREC[n], pdir, P);
2304 sprintf(title,
"%s (FREQ) : fREC (red) - fRMS:2fRMS (grey:light_gray)",NET->
ifoName[n]);
2307 PlotWaveformAsymmErrors(ofname,
"", cfg, &fREC[n], &fMED[n], &fL50[n], &fU50[n], &fL90[n], &fU90[n], &vREC[n], pdir, P,
true);
2314 sprintf(title,
"%s (TIME) : vREC(red) - vDAT-vREC(blue) - vRMS(green)",NET->
ifoName[n]);
2316 PlotWaveform(ofname,
"", cfg, &vREC[n], &wDIF, &vRMS[n], &vREC[n],
false, pdir, P);
2326 sprintf(title,
"%s (TIME) : vINJ (red) - vRMS:2vRMS (grey:light_gray)",NET->
ifoName[n]);
2329 PlotWaveformAsymmErrors(ofname,
"", cfg, &wALI, &vMED[n], &vL50[n], &vU50[n], &vL90[n], &vU90[n], &vREC[n], pdir, P);
2335 sprintf(title,
"%s (FREQ) : fINJ (red) - fRMS:2fRMS (grey:light_gray)",NET->
ifoName[n]);
2338 PlotWaveformAsymmErrors(ofname,
"", cfg, &fINJ[n], &fMED[n], &fL50[n], &fU50[n], &fL90[n], &fU90[n], &vREC[n], pdir, P,
true);
2350 sprintf(title,
"%s (TIME) : vAVR-vINJ(red) - 2*vRMS(gray)",NET->
ifoName[n]);
2353 PlotWaveformAsymmErrors(ofname,
"", cfg, &wDIF, &wDIF, &vL50[n], &vU50[n], &vL90[n], &vU90[n], &vREC[n], pdir, P);
2370 sprintf(title,
"%s (TIME) : vAVR-cINJ(red) - 2*vRMS(gray)",NET->
ifoName[n]);
2373 PlotWaveformAsymmErrors(ofname,
"", cfg, &wDIF, &wDIF, &vL50[n], &vU50[n], &vL90[n], &vU90[n], &vREC[n], pdir, P);
2379 sprintf(title,
"%s (TIME) : cINJ (red) - vRMS:2vRMS (grey:light_gray)",NET->
ifoName[n]);
2382 PlotWaveformAsymmErrors(ofname,
"", cfg, &cINJ[n], &vMED[n], &vL50[n], &vU50[n], &vL90[n], &vU90[n], &vREC[n], pdir, P);
2384 PlotWaveformErrors(ofname,
"", cfg, &cINJ[n], &vAVR[n], &vRMS[n], &vREC[n], pdir, P);
2388 sprintf(title,
"%s (FREQ) : fINJ (red) - fRMS:2fRMS (grey:light_gray)",NET->
ifoName[n]);
2391 PlotWaveformAsymmErrors(ofname,
"", cfg, &fINJ[n], &fMED[n], &fL50[n], &fU50[n], &fL90[n], &fU90[n], &vREC[n], pdir, P,
true);
2397 sprintf(title,
"%s (TIME) : cINJ(red) - vREC(blue) - vRMS(green)",NET->
ifoName[n]);
2399 PlotWaveform(ofname,
"", cfg, &cINJ[n], &vREC[n], &vRMS[n], &vREC[n],
false, pdir, P);
2405 sprintf(title,
"%s (TIME) : vINJ(red) - cINJ(blue) - vINJ-cINJ(green)",NET->
ifoName[n]);
2407 PlotWaveform(ofname,
"", cfg, &wALI, &cINJ[n], &wDIF, &vREC[n],
false, pdir, P);
2425 cout << endl <<
"SetSkyMask : " << cfg->
skyMaskFile << endl << endl;
2432 if(type!=
'j' && type!=
'r' && type!=
'v') {
2433 cout <<
"AddNoiseAndCalErrToSparse - Error : wrong wave type" << endl;
exit(1);
2459 switch(gOPT.noise) {
2461 for(
int i=0;
i<WF[
n].
size();
i++) {
2464 cout <<
"AddNoiseAndCalErrToSparse - Info : add gaussian noise to waveform type : " << type << endl;
2473 k = gRandom->Uniform(jS+jW,jE-jW);
2475 for(
int i=jS;
i<jE;
i++) {
2483 ts = gRandom->Uniform(0.,dt);
2484 for (
int j=0;
j<size/2;
j++) {
2485 TComplex X(WF[
n].data[2*
j],WF[
n].data[2*j+1]);
2486 X=X*C.Exp(TComplex(0.,-2*
PI*j*df*ts));
2488 WF[
n].
data[2*j+1]=X.Im();
2491 printf(
"Info : %s add data noise with time shift : %.5f sec to waveform type : %c\n",NET->
getifo(
n)->
Name,(double)k/rate+ts,type);
2502 if(gOPT.amp_cal_err[
n]) {
2503 double amp_cal_err = 1;
2504 if(gOPT.amp_cal_err[
n]>0) amp_cal_err = gRandom->Uniform(1-gOPT.amp_cal_err[
n],1+gOPT.amp_cal_err[
n]);
2505 else amp_cal_err = gRandom->Gaus(1,
fabs(gOPT.amp_cal_err[
n]));
2506 cout << NET->
ifoName[
n] <<
" -> amp_cal_err : " << amp_cal_err << endl;
2510 if(gOPT.phs_cal_err[
n]) {
2511 double phs_cal_err = 0;
2512 if(gOPT.phs_cal_err[
n]>0) phs_cal_err = gRandom->Uniform(-gOPT.phs_cal_err[
n],gOPT.phs_cal_err[
n]);
2513 else phs_cal_err = gRandom->Gaus(0,
fabs(gOPT.phs_cal_err[
n]));
2514 cout << NET->
ifoName[
n] <<
" -> phs_cal_err : " << phs_cal_err << endl;
2515 C = cos(-phs_cal_err*d2r);
2516 S = sin(-phs_cal_err*d2r);
2526 aa = rSS[
n][
k].sparseMap00[
l];
2527 AA = rSS[
n][
k].sparseMap90[
l];
2530 aa = jSS[
n][
k].sparseMap00[
l];
2531 AA = jSS[
n][
k].sparseMap90[
l];
2534 aa = vSS[
n][
k].sparseMap00[
l];
2535 AA = vSS[
n][
k].sparseMap90[
l];
2540 double bb = C*aa +
S*AA;
2541 double BB = -
S*aa + C*AA ;
2553 if(type!=
'j' && type!=
'r' && type!=
'v' && type!=
'd') {
2554 cout <<
"Wave2Sparse - Error : wrong wave type" << endl;
exit(1);
2561 if(type==
'v') vSS[
n]=pD->
vSS;
2562 if(type==
'r') rSS[
n]=pD->
vSS;
2563 if(type==
'j') jSS[
n]=pD->
vSS;
2564 if(type==
'd') dSS[
n]=pD->
vSS;
2566 if(type==
'v')
return;
2574 WF[
n].
rate(vSS[
n][0].wdm_rate);
2575 WF[
n].
start(vSS[
n][0].wdm_start);
2579 for (
int i=0;
i<vINJ[
n].size();
i++) WF[
n][wos+
i] = vINJ[
n][
i];
2583 for (
int i=0;
i<vREC[
n].size();
i++) WF[
n][wos+
i] = vREC[
n][
i];
2587 for (
int i=0;
i<vDAT[
n].size();
i++) WF[
n][wos+
i] = vDAT[
n][
i];
2590 #ifdef SAVE_WAVE2SPARSE_PLOT
2595 if(type==
'r') gfile=TString::Format(
"%s/WAVE2SPARSE_REC_%s.root",
".",NET->
ifoName[n]);
2596 if(type==
'j') gfile=TString::Format(
"%s/WAVE2SPARSE_INJ_%s.root",
".",NET->
ifoName[n]);
2597 if(type==
'd') gfile=TString::Format(
"%s/WAVE2SPARSE_DAT_%s.root",
".",NET->
ifoName[n]);
2612 if(type==
'r') size = rSS[
n][
k].sparseMap00.size();
2613 if(type==
'j') size = jSS[
n][
k].sparseMap00.size();
2614 if(type==
'd') size = dSS[
n][
k].sparseMap00.size();
2618 if(type==
'r') index = rSS[
n][
k].sparseIndex[
j];
2619 if(type==
'j') index = jSS[
n][
k].sparseIndex[
j];
2620 if(type==
'd') index = dSS[
n][
k].sparseIndex[
j];
2624 rSS[
n][
k].sparseMap00[
j]=v00;
2625 rSS[
n][
k].sparseMap90[
j]=v90;
2628 jSS[
n][
k].sparseMap00[
j]=v00;
2629 jSS[
n][
k].sparseMap90[
j]=v90;
2632 dSS[
n][
k].sparseMap00[
j]=v00;
2633 dSS[
n][
k].sparseMap90[
j]=v90;
2643 std::vector<int> nrec(nIFO);
2645 #ifdef SET_WAVEFORM_CUTS
2653 if(bT<gBT) gBT=bT;
if(eT>gET) gET=eT;
2654 if(bF<gBF) gBF=bF;
if(eF>gEF) gEF=eF;
2657 cout <<
"CUT TIME RANGE : " << gBT-gBT <<
" " << gET-gBT << endl;
2658 cout <<
"CUT FREQ RANGE : " << gBF <<
" " << gEF << endl;
2669 #ifdef SET_WAVEFORM_CUTS
2681 for(
int i=0;
i<gOPT.trials;
i++) {
2682 if(wREC[
i].
size()) {
2685 #ifdef SET_WAVEFORM_CUTS
2690 wrec[
i][
n] = vREC[
n]; wrec[
i][
n]=0;
2691 frec[
i][
n] = wrec[
i][
n];
2694 if(nrec[
n]==0) {nrec.clear();
return nrec;}
2697 gPE.trials = nrec[0];
2698 cout << endl <<
"ComputeStatisticalErrors - detected/trials : " << nrec[0] <<
"/" << gOPT.trials << endl << endl;
2704 wrms[
n] = wrec[0][
n]; wrms[
n] = 0;
2705 wavr[
n] = wrec[0][
n]; wavr[
n] = 0;
2706 for(
int i=0;
i<gOPT.trials;
i++) {
2707 wavr[
n] += wrec[
i][
n];
2708 for(
int j=0;
j< wrec[
i][
n].
size();
j++) wrms[
n][
j] += wrec[
i][
n][
j]*wrec[
i][
n][
j];
2710 wavr[
n] *= 1./nrec[
n];
2712 wrms[
n] *= 1./nrec[
n];
2713 for(
int j=0;
j< wrms[
n].
size();
j++) {
2714 wrms[
n][
j] -= wavr[
n][
j]*wavr[
n][
j];
2715 wrms[
n][
j] = sqrt(wrms[
n][
j]);
2718 vAVR.push_back(wavr[
n]);
2719 vRMS.push_back(wrms[n]);
2729 wmed[
n] = wrec[0][
n]; wmed[
n] = 0;
2730 wl50[
n] = wrec[0][
n]; wl50[
n] = 0;
2731 wu50[
n] = wrec[0][
n]; wu50[
n] = 0;
2732 wl90[
n] = wrec[0][
n]; wl90[
n] = 0;
2733 wu90[
n] = wrec[0][
n]; wu90[
n] = 0;
2735 int ntry = gPE.trials;
2736 int *
index =
new int[ntry];
2737 float *
value =
new float[ntry];
2738 for(
int j=0;
j<wrec[0][
n].
size();
j++) {
2740 int k=0;
for(
int i=0;
i<gOPT.trials;
i++)
if(wREC[
i].
size()) value[k++] = wrec[
i][
n][
j];
2741 TMath::Sort(ntry,value,index,
false);
2743 int imed = (ntry*50.)/100.;
if(imed>=ntry) imed=ntry-1;
2744 wmed[
n][
j] = value[index[imed]];
2746 int il50 = (ntry*25.)/100.;
if(il50>=ntry) il50=ntry-1;
2747 int iu50 = (ntry*75.)/100.;
if(iu50>=ntry) iu50=ntry-1;
2748 int il90 = (ntry*5.)/100.;
if(il90>=ntry) il90=ntry-1;
2749 int iu90 = (ntry*95.)/100.;
if(iu90>=ntry) iu90=ntry-1;
2751 wl50[
n][
j] = wmed[
n][
j]>0 ? wmed[
n][
j]-value[index[il50]] : value[index[iu50]]-wmed[
n][
j];
2752 wu50[
n][
j] = wmed[
n][
j]>0 ? value[index[iu50]]-wmed[
n][
j] : wmed[
n][
j]-value[index[il50]];
2753 wl90[
n][
j] = wmed[
n][
j]>0 ? wmed[
n][
j]-value[index[il90]] : value[index[iu90]]-wmed[
n][
j];
2754 wu90[
n][
j] = wmed[
n][
j]>0 ? value[index[iu90]]-wmed[
n][
j] : wmed[
n][
j]-value[index[il90]];
2759 vMED.push_back(wmed[
n]);
2760 vL50.push_back(wl50[n]);
2761 vU50.push_back(wu50[n]);
2762 vL90.push_back(wl90[n]);
2763 vU90.push_back(wu90[n]);
2770 frms[
n] = frec[0][
n]; frms[
n] = 0;
2771 favr[
n] = frec[0][
n]; favr[
n] = 0;
2772 for(
int i=0;
i<gOPT.trials;
i++) {
2773 favr[
n] += frec[
i][
n];
2774 for(
int j=0;
j< frec[
i][
n].
size();
j++) frms[
n][
j] += frec[
i][
n][
j]*frec[
i][
n][
j];
2776 favr[
n] *= 1./nrec[
n];
2778 frms[
n] *= 1./nrec[
n];
2779 for(
int j=0;
j< frms[
n].
size();
j++) {
2780 frms[
n][
j] -= favr[
n][
j]*favr[
n][
j];
2781 frms[
n][
j] = sqrt(frms[
n][
j]);
2784 fAVR.push_back(favr[
n]);
2785 fRMS.push_back(frms[n]);
2795 fmed[
n] = frec[0][
n]; fmed[
n] = 0;
2796 fl50[
n] = frec[0][
n]; fl50[
n] = 0;
2797 fu50[
n] = frec[0][
n]; fu50[
n] = 0;
2798 fl90[
n] = frec[0][
n]; fl90[
n] = 0;
2799 fu90[
n] = frec[0][
n]; fu90[
n] = 0;
2801 int ntry = gPE.trials;
2802 int *
index =
new int[ntry];
2803 float *
value =
new float[ntry];
2804 for(
int j=0;
j<frec[0][
n].
size();
j++) {
2806 int k=0;
for(
int i=0;
i<gOPT.trials;
i++)
if(wREC[
i].
size()) value[k++] = frec[
i][
n][
j];
2807 TMath::Sort(ntry,value,index,
false);
2809 int imed = (ntry*50.)/100.;
if(imed>=ntry) imed=ntry-1;
2810 fmed[
n][
j] = value[index[imed]];
2812 int il50 = (ntry*25.)/100.;
if(il50>=ntry) il50=ntry-1;
2813 int iu50 = (ntry*75.)/100.;
if(iu50>=ntry) iu50=ntry-1;
2814 int il90 = (ntry*5.)/100.;
if(il90>=ntry) il90=ntry-1;
2815 int iu90 = (ntry*95.)/100.;
if(iu90>=ntry) iu90=ntry-1;
2817 fl50[
n][
j] = fmed[
n][
j]>0 ? fmed[
n][
j]-value[index[il50]] : value[index[iu50]]-fmed[
n][
j];
2818 fu50[
n][
j] = fmed[
n][
j]>0 ? value[index[iu50]]-fmed[
n][
j] : fmed[
n][
j]-value[index[il50]];
2819 fl90[
n][
j] = fmed[
n][
j]>0 ? fmed[
n][
j]-value[index[il90]] : value[index[iu90]]-fmed[
n][
j];
2820 fu90[
n][
j] = fmed[
n][
j]>0 ? value[index[iu90]]-fmed[
n][
j] : fmed[
n][
j]-value[index[il90]];
2825 fMED.push_back(fmed[
n]);
2826 fL50.push_back(fl50[n]);
2827 fU50.push_back(fu50[n]);
2828 fL90.push_back(fl90[n]);
2829 fU90.push_back(fu90[n]);
2837 for(
int j=0;
j< wavr[
n].
size();
j++) eavr += pow(wavr[
n][
j],2);
2839 for(
int i=0;
i<gOPT.trials;
i++) {
2840 if(wREC[
i].
size()) {
2843 for(
int j=0;
j< wavr[
n].
size();
j++) {
2844 eres += pow(wrec[
i][
n][
j]-wavr[
n][
j],2);
2848 vRES.push_back(eres/eavr);
2854 for(
int j=0;
j< wavr[
n].
size();
j++) {
2855 jres += pow(winj[
n][
j]-wavr[
n][
j],2);
2861 vRES.push_back(jres);
2864 int size = vRES.size()-1;
2867 for(
int i=0;
i<
size;
i++) eres[
i] = vRES[
i];
2868 TMath::Sort(size,const_cast<double*>(eres.
data),index.
data,
false);
2870 for(
int i=1;
i<10;
i++) gPE.erR[
i]=eres[index[
int(
i*size/10.)]];
2873 if(eres[
i]<=jres) gPE.erR[10]+=1.;
2878 cout << endl <<
"gPE.erR : ";
2879 for(
int i=0;
i<10;
i++) cout << gPE.erR[
i] <<
"/";
2880 cout << gPE.erR[10] << endl << endl;
2886 for(
int n=0;
n<nIFO;
n++) {
2887 for(
int j=0;
j< wavr[
n].
size();
j++) eavr += pow(wavr[
n][
j],4)*pow(favr[
n][j],2);
2889 for(
int i=0;
i<gOPT.trials;
i++) {
2890 if(wREC[
i].
size()) {
2893 for(
int j=0;
j<wavr[
n].
size();
j++) {
2894 eres += pow(wavr[
n][
j],4)*pow(frec[
i][
n][j]-favr[
n][j],2);
2898 fRES.push_back(eres/eavr);
2904 for(
int j=0;
j< wavr[
n].
size();
j++) {
2905 jres += pow(wavr[
n][
j],4)*pow(finj[
n][j]-favr[
n][j],2);
2911 fRES.push_back(jres);
2914 for(
int i=0;
i<
size;
i++) eres[
i] = fRES[
i];
2915 TMath::Sort(size,const_cast<double*>(eres.
data),index.
data,
false);
2917 for(
int i=1;
i<10;
i++) gPE.erF[
i]=eres[index[
int(
i*size/10.)]];
2920 if(eres[
i]<=jres) gPE.erF[10]+=1.;
2925 cout << endl <<
"gPE.erF : ";
2926 for(
int i=0;
i<10;
i++) cout << gPE.erF[
i] <<
"/";
2927 cout << gPE.erF[10] << endl << endl;
2937 if(wf1==NULL)
return wf;
2938 if(wf1->
size()==0)
return wf;
2940 double R=wf1->
rate();
2942 double b_wf1 = wf1->
start();
2943 double e_wf1 = wf1->
start()+wf1->
size()/R;
2944 double b_wf2 = wf2->
start();
2945 double e_wf2 = wf2->
start()+wf2->
size()/R;
2947 int o_wf1 = b_wf1>b_wf2 ? 0 :
int((b_wf2-b_wf1)*R+0.5);
2948 int o_wf2 = b_wf1<b_wf2 ? 0 :
int((b_wf1-b_wf2)*R+0.5);
2950 double startXCOR = b_wf1>b_wf2 ? b_wf1 : b_wf2;
2951 double endXCOR = e_wf1<e_wf2 ? e_wf1 : e_wf2;
2952 int sizeXCOR =
int((endXCOR-startXCOR)*R+0.5);
2954 for(
int i=0;
i<sizeXCOR;
i++) wf[
i+o_wf2] = wf1->
data[
i+o_wf1];
2963 if(wf==NULL)
return wfq;
2964 if(wf->
size()==0)
return wfq;
2969 for(
int i=0;
i<wf->
size();
i++) wfq[
i] = sqrt(pow(wf->
data[
i],2)+pow(wfq[
i],2));
2976 std::vector<float>* vP;
2977 std::vector<int>* vI;
2983 if(NET->
wc_List[0].p_Map.size()) {
2985 vP = &(NET->
wc_List[0].p_Map[
id-1]);
2986 vI = &(NET->
wc_List[0].p_Ind[
id-1]);
2991 for(
int j=0;
j<
int(vP->size());
j++) {
2996 skyprob.
set(k,(*vP)[
j]);
3005 std::vector<float>* vP;
3006 std::vector<int>* vI;
3009 if(NET->
wc_List[0].p_Map.size()) {
3011 vP = &(NET->
wc_List[0].p_Map[
id-1]);
3012 vI = &(NET->
wc_List[0].p_Ind[
id-1]);
3017 for(
int j=0;
j<
int(vP->size());
j++) {
3022 gSKYPROB.
set(k,(*vP)[
j]);
3026 gHSKYPROB.SetBins(gSKYPROB.
size(),0,gSKYPROB.
size()-1);
3029 for(
int l=0;
l<gSKYPROB.
size();
l++) PROB+=gSKYPROB.
get(
l);
3030 cout <<
"PROB : " << PROB << endl;
3031 for(
int l=0;
l<gSKYPROB.
size();
l++) gHSKYPROB.SetBinContent(
l,gSKYPROB.
get(
l));
3061 gSM.
Print(TString::Format(
"%s/gSkyprob_%d.png",
".",
id));
3067 std::vector<double > xRES = type==
'f' ? fRES :
vRES;
3069 int size = xRES.size()-1;
3070 double jres = xRES[
size];
3075 for(
int i=0;
i<
size;
i++) {
if(xRES[
i]>hmax) hmax=xRES[
i];
if(xRES[
i]<hmin) hmin=xRES[
i];}
3076 hmin-=0.4*(hmax-hmin);
3077 hmax+=0.4*(hmax-hmin);
3078 hmax*=1.2;
if(hmax>1) hmax=1;
3081 TCanvas* cres =
new TCanvas(
"residuals",
"residuals", 200, 20, 600, 600);
3082 TH1F* hres =
new TH1F(
"residuals",
"residuals",100,hmin,hmax);
3083 for(
int i=0;
i<
size;
i++) hres->Fill(xRES[
i]);
3086 double integral = hres->ComputeIntegral();
3087 if (integral==0) {cout <<
"PlotResiduals : Empty histogram !!!" << endl;
return;}
3088 double* cumulative = hres->GetIntegral();
3089 for (
int i=0;i<hres->GetNbinsX();i++) hres->SetBinContent(i,cumulative[i]/integral);
3091 hres->SetLineWidth(2);
3095 hres->GetXaxis()->SetTitle(
"residuals");
3096 hres->GetXaxis()->SetTitleOffset(1.4);
3097 hres->GetXaxis()->SetLabelOffset(0.02);
3098 hres->GetXaxis()->SetNoExponent();
3099 hres->GetXaxis()->SetMoreLogLabels();
3101 hres->GetXaxis()->SetRangeUser(0.001,1.0);
3103 hres->GetXaxis()->SetRangeUser(0.01,1.0);
3106 hres->GetYaxis()->SetTitle(
"probability");
3107 hres->GetYaxis()->SetTitleOffset(1.6);
3108 hres->GetYaxis()->SetLabelOffset(0.02);
3109 hres->GetYaxis()->SetNoExponent();
3110 hres->GetYaxis()->SetRangeUser(0.01,1.0);
3118 float ymax = hres->GetMaximum();
3119 TLine *
line =
new TLine(jres,0,jres,ymax);
3120 line->SetLineWidth(2);
3121 line->SetLineColor(kRed);
3122 if(sim) line->Draw();
3127 cout << endl <<
"---------> gPE.erF : ";
3128 for(
int i=0;i<10;i++) cout << gPE.erF[i] <<
"/";
3129 cout << gPE.erF[10] << endl << endl;
3131 cout << endl <<
"---------> gPE.erR : ";
3132 for(
int i=0;i<10;i++) cout << gPE.erR[i] <<
"/";
3133 cout << gPE.erR[10] << endl << endl;
3136 gStyle->SetLineColor(kBlack);
3139 sprintf(title,
"Cumulative distribution of freq residuals errors (entries=%d - prob inj = %2.2g)",size,gPE.erF[10]);
3141 sprintf(title,
"Cumulative distribution of residuals errors (entries=%d - prob inj = %2.2g)",size,gPE.erR[10]);
3143 hres->SetTitle(title);
3144 hres->SetStats(kFALSE);
3146 cres->Print(TString::Format(
"%s/fresiduals_errors.%s",pdir.Data(),
PLOT_TYPE));
3148 cres->Print(TString::Format(
"%s/residuals_errors.%s",pdir.Data(),
PLOT_TYPE));
3158 int size = sim ? iSNR[0].size() : vLIKELIHOOD.size();
3165 if(sim)
for(
int n=0;
n<
nIFO;
n++) osnr+=oSNR[
n][
i];
3166 else osnr+=vLIKELIHOOD[
i];
3168 if(osnr>hmax) hmax=osnr;
3169 if(osnr<hmin) hmin=osnr;
3171 hmin-=0.4*(hmax-hmin);
3172 hmax+=0.4*(hmax-hmin);
3176 for(
int n=0;
n<
nIFO;
n++) isnr+=iSNR[
n][0];
3178 if(isnr<hmin) hmin=0.8*isnr;
3179 if(isnr>hmax) hmax=1.2*isnr;
3183 TCanvas* csnr =
new TCanvas(
"SNRnet",
"SNRnet", 200, 20, 600, 600);
3184 TH1F* hsnr =
new TH1F(
"SNRnet",
"SNRnet",10,hmin,hmax);
3187 if(sim)
for(
int n=0;
n<
nIFO;
n++) osnr+=oSNR[
n][
i];
3188 else osnr+=vLIKELIHOOD[
i];
3189 hsnr->Fill(sqrt(osnr));
3191 hsnr->SetLineWidth(2);
3195 hsnr->GetXaxis()->SetTitle(
"SNRnet");
3196 hsnr->GetXaxis()->SetTitleOffset(1.4);
3197 hsnr->GetXaxis()->SetLabelOffset(0.02);
3198 hsnr->GetXaxis()->SetNoExponent();
3200 hsnr->GetYaxis()->SetTitle(
"counts");
3201 hsnr->GetYaxis()->SetTitleOffset(1.6);
3202 hsnr->GetYaxis()->SetLabelOffset(0.02);
3203 hsnr->GetYaxis()->SetNoExponent();
3210 float ymax = hsnr->GetMaximum();
3211 TLine *
line =
new TLine(isnr,0,isnr,ymax);
3212 line->SetLineWidth(2);
3213 line->SetLineColor(kRed);
3214 if(sim) line->Draw();
3216 hsnr->Fit(
"gaus",
"q0");
3217 TF1* gauss = (TF1*)hsnr->GetFunction(
"gaus");
3219 gauss->Draw(
"SAME");
3220 gauss->SetLineColor(kGreen);
3221 gauss->SetLineWidth(2);
3225 if(sim)
sprintf(title,
"SNRnet Distribution (entries=%d - SNRnet inj = %2.2g)",size,isnr);
3226 else sprintf(title,
"SNRnet Distribution (entries=%d)",size);
3227 hsnr->SetTitle(title);
3228 hsnr->SetStats(kTRUE);
3229 csnr->Print(TString::Format(
"%s/SNRnet.%s",pdir.Data(),
PLOT_TYPE));
3230 hsnr->SetStats(kFALSE);
3245 if(gtype<1 || gtype>5)
return;
3248 if(gtype==1) gname=
"chirp_mass";
3249 if(gtype==2) gname=
"chirp_mass_err";
3250 if(gtype==3) gname=
"chirp_ell";
3251 if(gtype==4) gname=
"chirp_pfrac";
3252 if(gtype==5) gname=
"chirp_efrac";
3254 int size = vCHIRP[gtype].size();
3260 double omchirp = vCHIRP[gtype][
i];
3261 if(omchirp>hmax) hmax=omchirp;
3262 if(omchirp<hmin) hmin=omchirp;
3264 hmin-=0.4*(hmax-hmin);
3265 hmax+=0.4*(hmax-hmin);
3267 double imchirp=vCHIRP[0][0];
3268 if(sim && gtype==1) {
3269 if(imchirp<hmin) hmin=0.8*imchirp;
3270 if(imchirp>hmax) hmax=1.2*imchirp;
3274 TCanvas* cmchirp =
new TCanvas(
"mchirp",
"mchirp", 200, 20, 600, 600);
3275 TH1F* hmchirp =
new TH1F(
"mchirp",
"mchirp",10,hmin,hmax);
3276 for(
int i=0;
i<
size;
i++) hmchirp->Fill(vCHIRP[gtype][
i]);
3277 hmchirp->SetLineWidth(2);
3281 hmchirp->GetXaxis()->SetTitle(gname);
3282 hmchirp->GetXaxis()->SetTitleOffset(1.4);
3283 hmchirp->GetXaxis()->SetLabelOffset(0.02);
3284 hmchirp->GetXaxis()->SetNoExponent();
3286 hmchirp->GetYaxis()->SetTitle(
"counts");
3287 hmchirp->GetYaxis()->SetTitleOffset(1.6);
3288 hmchirp->GetYaxis()->SetLabelOffset(0.02);
3289 hmchirp->GetYaxis()->SetNoExponent();
3292 cmchirp->SetGridx();
3293 cmchirp->SetGridy();
3296 float ymax = hmchirp->GetMaximum();
3297 TLine *
line =
new TLine(imchirp,0,imchirp,ymax);
3298 line->SetLineWidth(2);
3299 line->SetLineColor(kRed);
3300 if(sim && gtype==1) line->Draw();
3310 gStyle->SetLineColor(kBlack);
3312 if(sim && gtype==1)
sprintf(title,
"%s distribution (entries=%d - inj chirp mass = %2.2g)",gname.Data(),
size,imchirp);
3313 else sprintf(title,
"%s distribution (entries=%d)",gname.Data(),
size);
3314 hmchirp->SetTitle(title);
3315 hmchirp->SetStats(kTRUE);
3316 cmchirp->Print(TString::Format(
"%s/%s.%s",pdir.Data(),gname.Data(),
PLOT_TYPE));
3317 hmchirp->SetStats(kFALSE);
3329 if(gtype!=0 && gtype!=1)
return;
3331 int size = iSNR[0].size();
3337 double isnr=0;
for(
int n=0;
n<
nIFO;
n++) isnr+= iSNR[
n][
i];
3338 double osnr=0;
for(
int n=0;
n<
nIFO;
n++) osnr+= oSNR[
n][
i];
3339 double iosnr=0;
for(
int n=0;
n<
nIFO;
n++) iosnr+=ioSNR[
n][
i];
3340 double ff = iosnr/sqrt(isnr*isnr);
3341 double of = iosnr/sqrt(isnr*osnr);
3342 double factor = gtype==0 ? ff : of;
3343 if(factor>hmax) hmax=
factor;
if(factor<hmin) hmin=
factor;
3345 hmin-=0.4*(hmax-hmin);
3346 hmax+=0.4*(hmax-hmin);
3349 TCanvas* cf =
new TCanvas(
"Factors",
"Factors", 200, 20, 600, 600);
3350 TH1F* hf =
new TH1F(
"Factors",
"Factors",10,hmin,hmax);
3352 double isnr=0;
for(
int n=0;
n<
nIFO;
n++) isnr+= iSNR[
n][
i];
3353 double osnr=0;
for(
int n=0;
n<
nIFO;
n++) osnr+= oSNR[
n][
i];
3354 double iosnr=0;
for(
int n=0;
n<
nIFO;
n++) iosnr+=ioSNR[
n][
i];
3355 double ff = iosnr/sqrt(isnr*isnr);
3356 double of = iosnr/sqrt(isnr*osnr);
3357 double factor = gtype==0 ? ff : of;
3360 hf->SetLineWidth(2);
3364 if(gtype==0) hf->GetXaxis()->SetTitle(
"Fitting Factor");
3365 if(gtype==1) hf->GetXaxis()->SetTitle(
"Overlap Factor");
3366 hf->GetXaxis()->SetTitleOffset(1.4);
3367 hf->GetXaxis()->SetLabelOffset(0.02);
3368 hf->GetXaxis()->SetNoExponent();
3370 hf->GetYaxis()->SetTitle(
"counts");
3371 hf->GetYaxis()->SetTitleOffset(1.6);
3372 hf->GetYaxis()->SetLabelOffset(0.02);
3373 hf->GetYaxis()->SetNoExponent();
3380 if(gtype==0)
sprintf(title,
"Fitting Factor Distribution (entries=%d)",size);
3381 if(gtype==1)
sprintf(title,
"Overlap Factor Distribution (entries=%d)",size);
3382 hf->SetTitle(title);
3383 gStyle->SetLineColor(kBlack);
3384 hf->SetStats(kTRUE);
3385 if(gtype==0) cf->Print(TString::Format(
"%s/FittingFactor.%s",pdir.Data(),
PLOT_TYPE));
3386 if(gtype==1) cf->Print(TString::Format(
"%s/OverlapFactor.%s",pdir.Data(),
PLOT_TYPE));
3387 hf->SetStats(kFALSE);
3395 if(vNUL->size()==0 || vNSE->size()==0)
return;
3399 for(
int i=0;
i<vNUL->size();
i++) {
3400 if((*vNUL)[
i]>xmax) xmax=(*vNUL)[
i];
if((*vNUL)[
i]<xmin) xmin=(*vNUL)[
i];
3402 xmin *= xmin>0 ? 0.5 : 1.5;
3403 xmax *= xmax>0 ? 1.5 : 0.5;
3406 TH1F* hnul =
new TH1F(
"null",
"null",10,xmin,xmax);
3407 for(
int i=0;
i<vNUL->size();
i++) hnul->Fill((*vNUL)[
i]);
3408 hnul->SetLineWidth(2);
3409 hnul->SetLineColor(kRed);
3411 TH1F* hnse =
new TH1F(
"noise",
"noise",10,xmin,xmax);
3412 for(
int i=0;i<vNSE->size();i++) hnse->Fill((*vNSE)[i]);
3413 double scale = double(hnul->GetEntries())/hnse->GetEntries();
3415 hnse->SetLineWidth(2);
3416 hnse->SetLineColor(kBlack);
3420 for(
int i=0;i<=hnul->GetNbinsX();i++) {
3421 double binc = hnul->GetBinContent(i);
3422 if(binc>ymax) ymax=binc;
3424 for(
int i=0;i<=hnse->GetNbinsX();i++) {
3425 double binc = hnse->GetBinContent(i);
3426 if(binc>ymax) ymax=binc;
3429 hnse->GetYaxis()->SetRangeUser(ymin,ymax);
3431 double KS = hnul->Integral() ? hnul->KolmogorovTest(hnse,
"N") : 0.;
3432 double AD = hnul->Integral() ? hnul->AndersonDarlingTest(hnse) : 0;
3434 hnul->Fit(
"gaus",
"q0");
3435 TF1* gaus = (TF1*)hnul->GetFunction(
"gaus");
3436 double chi2 = gaus ? gaus->GetChisquare()/gaus->GetNDF() : 0;
3440 TCanvas* cnul =
new TCanvas(
"cnull",
"cnull", 200, 20, 600, 600);
3445 hnse->GetXaxis()->SetTitle(gtype);
3446 hnse->GetXaxis()->SetTitleOffset(1.4);
3447 hnse->GetXaxis()->SetLabelOffset(0.02);
3448 hnse->GetXaxis()->SetNoExponent();
3450 hnse->GetYaxis()->SetTitle(
"counts");
3451 hnse->GetYaxis()->SetTitleOffset(1.6);
3452 hnse->GetYaxis()->SetLabelOffset(0.02);
3453 hnse->GetYaxis()->SetNoExponent();
3459 hnse->Fit(
"gaus",
"q0");
3460 gaus = (TF1*)hnse->GetFunction(
"gaus");
3463 gaus->SetLineColor(kGreen);
3464 gaus->SetLineWidth(2);
3468 sprintf(title,
"%s : %s (entries=%d) : KS=%0.2f : AD=%0.2f : CHI2=%0.2f",ifoName.Data(),gtype.Data(),vNUL->size(),KS,AD,
chi2);
3469 hnse->SetTitle(title);
3470 hnse->SetStats(kTRUE);
3471 cnul->Print(TString::Format(
"%s/%s_%s_dist.%s",pdir.Data(),ifoName.Data(),gtype.Data(),
PLOT_TYPE));
3472 hnse->SetStats(kFALSE);
3481 if(gtype==
"null_90") I+=7;
3483 gPE.nstat[I+0] = hnul->GetMean();
3484 gPE.nstat[I+1] = hnul->GetRMS();
3485 gPE.nstat[I+2] = hnse->GetMean();
3486 gPE.nstat[I+3] = hnse->GetRMS();
3487 gPE.nstat[I+4] =
chi2;
3488 gPE.nstat[I+5] = KS;
3489 gPE.nstat[I+6] = AD;
3492 cout <<
" Pixels Statistic : " << ifoName <<
" " << gtype << endl;
3493 cout <<
" Null Pixels Mean : " << gPE.nstat[I+0] << endl;
3494 cout <<
" Null Pixels RMS : " << gPE.nstat[I+1] << endl;
3495 cout <<
" Noise Pixels Mean : " << gPE.nstat[I+2] << endl;
3496 cout <<
" Noise Pixels RMS : " << gPE.nstat[I+3] << endl;
3497 cout <<
" Noise Pixels Chi2 : " << gPE.nstat[I+4] << endl;
3498 cout <<
" KolmogorovTest : " << gPE.nstat[I+5] << endl;
3499 cout <<
" AndersonDarlingTest : " << gPE.nstat[I+6] << endl;
3518 T = E>0 ?
T/E : 0.5*size/
rate;
3532 for(
int i=0;
i<
N;
i++) {avr+=
i*x[
i]*x[
i]; E+=x[
i]*x[
i];}
3539 double sum = ((M>=0)&&(M<N)) ? x[M]*x[M] : 0.;
3540 for(
int j=1;
j<
N;
j++) {
3541 a = ((M-
j>=0)&&(M-j<N)) ? x[M-j] : 0.;
3542 b = ((M+j>=0)&&(M+j<N)) ? x[M+j] : 0.;
3546 if(sum/E > P)
break;
3562 double dF=(rate/(double)size)/2.;
3564 a = x[
j]*x[
j]+x[
j+1]*x[
j+1];
3568 F = E>0 ?
F/E : 0.5*
rate;
3584 for(
int i=0;
i<
N;
i++) {avr+=
i*x[
i]*x[
i]; E+=x[
i]*x[
i];}
3591 double sum = ((M>=0)&&(M<N)) ? x[M]*x[M] : 0.;
3592 for(
int j=1;
j<
N;
j++) {
3593 a = ((M-
j>=0)&&(M-j<N)) ? x[M-j] : 0.;
3594 b = ((M+j>=0)&&(M+j<N)) ? x[M+j] : 0.;
3598 if(sum/E > P)
break;
3601 double dF=(x.
rate()/(double)x.
size())/2.;
3621 if(T<bT || T>eT) x[
j]=0;
3626 double dF=(x.
rate()/(double)x.
size())/2.;
3630 if(F<bF || F>eF) {x[
j]=0;x[
j+1]=0;}
3646 double dT = 1./x->
rate();
3649 if(T<bT || T>eT) x->
data[
j]=0;
3654 double dF=(x->
rate()/(double)x->
size())/2.;
3658 if(F<bF || F>eF) {x->
data[
j]=0;x->
data[
j+1]=0;}
3674 if(factor==0)
sprintf(sfactor,
"z%g",factor);
3675 if(factor>0)
sprintf(sfactor,
"p%g",factor);
3676 }
else sprintf(sfactor,
"%g",factor);
3680 char sim_label[512];
3682 sprintf(out_CED,
"%s/ced_%s_%d",cfg->
nodedir,sim_label,gSystem->GetPid());
3684 char prod_label[512];
3685 sprintf(prod_label,
"%d_%d_%s_slag%d_lag%lu_%lu_job%d",
3687 sprintf(out_CED,
"%s/ced_%s_%d",cfg->
nodedir,prod_label,gSystem->GetPid());
3689 cout << out_CED << endl;
3699 cout<<
"DumpCED : dump ced to disk ..." <<endl;
3715 ced.
Write(ofactor,fullCED);
3725 sprintf(dirCED,
"%s/%s", out_CED, ifostr);
3749 if(gSystem->Getenv(
"HOME_CED_WWW")==NULL) {
3750 cout <<
"Error : environment HOME_CED_WWW is not defined!!! " << endl;
exit(1);}
3753 sprintf(ofileName,
"%s/pe_index.html",dirCED.Data());
3754 cout <<
"ofileName : " << ofileName << endl;
3761 if(!out.good()) {cout <<
"Error Opening File : " << ofileName << endl;
exit(1);}
3763 WriteBodyHTML(pe_index,
"PE_HEADER_BEG",
"PE_HEADER_END", &out);
3765 WriteBodyHTML(pe_index,
"PE_BODY_BEG",
"PE_BODY_END", &out);
3767 int sbody_height = 0;
3768 out <<
"<div class=\"tabber\">" << endl;
3772 if(gOPT.ced_tfmap) {
3773 out <<
"<div class=\"tabbertab\">" << endl;
3774 out <<
" <h2>Time-Frequency Maps</h2>" << endl;
3775 sbody_height = 1000+nIFO*400;
3776 if(gOPT.ced_pca) sbody_height += 600;
3777 out <<
" <iframe src=\"tfmap_body.html\" width=\"100%\" "
3778 <<
" height=\"" << sbody_height <<
"px\" frameborder=\"0\"></iframe>" << endl;
3779 out <<
"</div>" << endl;
3785 out <<
"<div class=\"tabbertab\">" << endl;
3786 out <<
" <h2>Reconstructed Detector Responses</h2>" << endl;
3787 if(!sim) sbody_height = 350+nIFO*400;
3788 else sbody_height = 350+nIFO*1340;
3789 out <<
" <iframe src=\"rec_signal_body.html\" width=\"100%\" "
3790 <<
" height=\"" << sbody_height <<
"px\" frameborder=\"0\"></iframe>" << endl;
3791 out <<
"</div>" << endl;
3797 out <<
"<div class=\"tabbertab\">" << endl;
3798 out <<
" <h2>PSD</h2>" << endl;
3799 sbody_height = 350+(nIFO+1)*750;
3800 out <<
" <iframe src=\"psd_body.html\" width=\"100%\" "
3801 <<
" height=\"" << sbody_height <<
"px\" frameborder=\"0\"></iframe>" << endl;
3802 out <<
"</div>" << endl;
3807 sbody_height = 2050;
3809 out <<
"<div class=\"tabbertab\">" << endl;
3810 out <<
" <h2>SkyMaps</h2>" << endl;
3811 out <<
" <iframe src=\"skymap_body.html\" width=\"100%\" "
3812 <<
" height=\"" << sbody_height <<
"px\" frameborder=\"0\"></iframe>" << endl;
3813 out <<
"</div>" << endl;
3819 out <<
"<div class=\"tabbertab\">" << endl;
3820 out <<
" <h2>Waveform Errors</h2>" << endl;
3821 sbody_height = 370+nIFO*880;
3822 out <<
" <iframe src=\"rec_avr_signal_body.html\" width=\"100%\" "
3823 <<
" height=\"" << sbody_height <<
"px\" frameborder=\"0\"></iframe>" << endl;
3824 out <<
"</div>" << endl;
3829 if(sim && gOPT.ced_inj) {
3830 out <<
"<div class=\"tabbertab\">" << endl;
3831 out <<
" <h2>Injection</h2>" << endl;
3832 sbody_height = 390+nIFO*880;
3833 out <<
" <iframe src=\"rec_inj_signal_body.html\" width=\"100%\" "
3834 <<
" height=\"" << sbody_height <<
"px\" frameborder=\"0\"></iframe>" << endl;
3835 out <<
"</div>" << endl;
3841 if(sim && gOPT.ced_rinj) {
3842 out <<
"<div class=\"tabbertab\">" << endl;
3843 out <<
" <h2>Injection</h2>" << endl;
3844 sbody_height = 390+nIFO*320;
3845 out <<
" <iframe src=\"rec_inj_signal_body.html\" width=\"100%\" "
3846 <<
" height=\"" << sbody_height <<
"px\" frameborder=\"0\"></iframe>" << endl;
3847 out <<
"</div>" << endl;
3853 sbody_height = 2100;
3854 out <<
"<div class=\"tabbertab\">" << endl;
3855 out <<
" <h2>Chirp Mass</h2>" << endl;
3856 out <<
" <iframe src=\"mchirp_body.html\" width=\"100%\" "
3857 <<
" height=\"" << sbody_height <<
"px\" frameborder=\"0\"></iframe>" << endl;
3858 out <<
"</div>" << endl;
3863 if(gOPT.ced_distr) {
3864 sbody_height = 2000;
3865 out <<
"<div class=\"tabbertab\">" << endl;
3866 out <<
" <h2>Distributions</h2>" << endl;
3867 out <<
" <iframe src=\"distribution_body.html\" width=\"100%\" "
3868 <<
" height=\"" << sbody_height <<
"px\" frameborder=\"0\"></iframe>" << endl;
3869 out <<
"</div>" << endl;
3875 sbody_height = 100+nIFO*1300;
3876 out <<
"<div class=\"tabbertab\">" << endl;
3877 out <<
" <h2>Null Pixels</h2>" << endl;
3878 out <<
" <iframe src=\"nstat_body.html\" width=\"100%\" "
3879 <<
" height=\"" << sbody_height <<
"px\" frameborder=\"0\"></iframe>" << endl;
3880 out <<
"</div>" << endl;
3885 out <<
"</div>" << endl;
3889 if(gOPT.ced_tfmap) {
3890 out.open(TString::Format(
"%s/tfmap_body.html",dirCED.Data()).Data(),
ios::out);
3892 WriteBodyHTML(pe_index,
"PE_JSCRIPT_BEG",
"PE_JSCRIPT_END", &out);
3893 WriteBodyHTML(pe_index,
"PE_TFMAP_HEADER_BEG",
"PE_TFMAP_HEADER_END", &out);
3894 WriteBodyHTML(pe_index,
"PE_TFMAP_BEG",
"PE_TFMAP_END", &out, nIFO, ifo);
3896 out <<
"<p><br /></p> " << endl;
3897 out <<
"<p><br /></p> " << endl;
3898 WriteBodyHTML(pe_index,
"PE_LIKELIHOOD_BEG",
"PE_LIKELIHOOD_END", &out);
3903 if(gOPT.ced_pca)
WriteBodyHTML(pe_index,
"PE_PCA_BEG",
"PE_PCA_END", &out);
3909 out.open(TString::Format(
"%s/rec_signal_body.html",dirCED.Data()).Data(),
ios::out);
3910 WriteBodyHTML(pe_index,
"PE_JSCRIPT_BEG",
"PE_JSCRIPT_END", &out);
3911 WriteBodyHTML(pe_index,
"PE_REC_SIGNAL_HEADER_BEG",
"PE_REC_SIGNAL_HEADER_END", &out);
3913 WriteBodyHTML(pe_index,
"PE_IFO_SIGNAL_BEG",
"PE_IFO_SIGNAL_END", &out, 1, &ifo[
n]);
3915 WriteBodyHTML(pe_index,
"PE_INJ_SIGNAL_BEG",
"PE_INJ_SIGNAL_END", &out, 1, &ifo[n]);
3917 WriteBodyHTML(pe_index,
"PE_REC_SIGNAL_BEG",
"PE_REC_SIGNAL_END", &out, 1, &ifo[n]);
3918 out <<
"<p><br /></p> " << endl;
3919 out <<
"<p><br /></p> " << endl;
3926 out.open(TString::Format(
"%s/psd_body.html",dirCED.Data()).Data(),
ios::out);
3927 WriteBodyHTML(pe_index,
"PE_JSCRIPT_BEG",
"PE_JSCRIPT_END", &out);
3928 WriteBodyHTML(pe_index,
"PE_PSD_HEADER_BEG",
"PE_PSD_HEADER_END", &out);
3930 WriteBodyHTML(pe_index,
"PE_IFO_PSD_BEG",
"PE_IFO_PSD_END", &out, 1, &ifo[
n]);
3931 WriteBodyHTML(pe_index,
"PE_PSD_BEG",
"PE_PSD_END", &out, 1, &ifo[n]);
3932 out <<
"<p><br /></p> " << endl;
3933 out <<
"<p><br /></p> " << endl;
3940 out.open(TString::Format(
"%s/skymap_body.html",dirCED.Data()).Data(),
ios::out);
3941 WriteBodyHTML(pe_index,
"PE_JSCRIPT_BEG",
"PE_JSCRIPT_END", &out);
3942 WriteBodyHTML(pe_index,
"PE_SKYMAP_BEG",
"PE_SKYMAP_END", &out);
3948 out.open(TString::Format(
"%s/rec_avr_signal_body.html",dirCED.Data()).Data(),
ios::out);
3949 WriteBodyHTML(pe_index,
"PE_JSCRIPT_BEG",
"PE_JSCRIPT_END", &out);
3950 WriteBodyHTML(pe_index,
"PE_REC_AVR_SIGNAL_HEADER_BEG",
"PE_REC_AVR_SIGNAL_HEADER_END", &out);
3952 WriteBodyHTML(pe_index,
"PE_REC_AVR_SIGNAL_BEG",
"PE_REC_AVR_SIGNAL_END", &out, 1, &ifo[
n]);
3953 out <<
"<p><br /></p> " << endl;
3954 out <<
"<p><br /></p> " << endl;
3960 if(sim && (gOPT.ced_inj || gOPT.ced_rinj)) {
3961 out.open(TString::Format(
"%s/rec_inj_signal_body.html",dirCED.Data()).Data(),
ios::out);
3962 WriteBodyHTML(pe_index,
"PE_JSCRIPT_BEG",
"PE_JSCRIPT_END", &out);
3963 if(gOPT.ced_inj)
WriteBodyHTML(pe_index,
"PE_REC_INJ_SIGNAL_0_HEADER_BEG",
"PE_REC_INJ_SIGNAL_0_HEADER_END", &out);
3964 if(gOPT.ced_rinj)
WriteBodyHTML(pe_index,
"PE_REC_INJ_SIGNAL_HEADER_BEG",
"PE_REC_INJ_SIGNAL_HEADER_END", &out);
3967 if(gOPT.ced_inj)
WriteBodyHTML(pe_index,
"PE_REC_INJ_SIGNAL_0_BEG",
"PE_REC_INJ_SIGNAL_0_END", &out, 1, &ifo[
n]);
3968 if(gOPT.ced_rinj)
WriteBodyHTML(pe_index,
"PE_REC_INJ_SIGNAL_BEG",
"PE_REC_INJ_SIGNAL_END", &out, 1, &ifo[n]);
3969 out <<
"<p><br /></p> " << endl;
3970 out <<
"<p><br /></p> " << endl;
3978 if(gOPT.ced_distr) {
3979 out.open(TString::Format(
"%s/distribution_body.html",dirCED.Data()).Data(),
ios::out);
3980 WriteBodyHTML(pe_index,
"PE_JSCRIPT_BEG",
"PE_JSCRIPT_END", &out);
3981 WriteBodyHTML(pe_index,
"PE_DISTRIBUTION_BEG",
"PE_DISTRIBUTION_END", &out);
3982 if(sim)
WriteBodyHTML(pe_index,
"PE_FACTORS_DISTRIBUTION_BEG",
"PE_FACTORS_DISTRIBUTION_END", &out);
3988 out.open(TString::Format(
"%s/nstat_body.html",dirCED.Data()).Data(),
ios::out);
3989 WriteBodyHTML(pe_index,
"PE_JSCRIPT_BEG",
"PE_JSCRIPT_END", &out);
3990 WriteBodyHTML(pe_index,
"PE_NULPIX_DISTRIBUTION_HEADER_BEG",
"PE_NULPIX_DISTRIBUTION_HEADER_END", &out);
3992 WriteBodyHTML(pe_index,
"PE_NULPIX_DISTRIBUTION_BEG",
"PE_NULPIX_DISTRIBUTION_END", &out, 1, &ifo[
n]);
3993 out <<
"<p><br /></p> " << endl;
3994 out <<
"<p><br /></p> " << endl;
4001 out.open(TString::Format(
"%s/mchirp_body.html",dirCED.Data()).Data(),
ios::out);
4002 WriteBodyHTML(pe_index,
"PE_JSCRIPT_BEG",
"PE_JSCRIPT_END", &out);
4003 WriteBodyHTML(pe_index,
"PE_MCHIRP_DISTRIBUTION_HEADER_BEG",
"PE_MCHIRP_DISTRIBUTION_HEADER_END", &out);
4004 WriteBodyHTML(pe_index,
"PE_MCHIRP_DISTRIBUTION_BEG",
"PE_MCHIRP_DISTRIBUTION_END", &out);
4010 sprintf(cmd,
"cp ${HOME_WAT}//html/etc/html/tabber.css %s",dirCED.Data());
4012 sprintf(cmd,
"cp ${HOME_WAT}//html/etc/html/tabber.js %s",dirCED.Data());
4015 sprintf(cmd,
"mv %s/pe_index.html %s/index.html",dirCED.Data(),dirCED.Data());
4030 if(gSystem->Getenv(
"CWB_DOC_URL")!=NULL) {
4031 cwb_doc_url=
TString(gSystem->Getenv(
"CWB_DOC_URL"));
4035 TString home_www=
"~waveburst/LSC/waveburst";
4036 if(gSystem->Getenv(
"HOME_WWW")!=NULL) {
4037 home_www=
TString(gSystem->Getenv(
"HOME_WWW"));
4044 in.open(html_template.Data(),
ios::in);
4045 if(!in.good()) {cout <<
"Error Opening File : " << html_template.Data() << endl;
exit(1);}
4047 in.getline(line,1024);
4048 if(!in.good())
break;
4050 if(sline.Contains(html_tag_beg)&&!
output) output=
true;
4051 if(sline.Contains(html_tag_end)&&
output) {output=
false;
break;}
4052 if(ifo!=NULL) sline.ReplaceAll(
"IFO",ifo[
n]);
4053 sline.ReplaceAll(
"xNTRIALS",TString::Format(
"%d",(
int)vCHIRP[0].
size()-1).Data());
4054 sline.ReplaceAll(
"xINRATE",TString::Format(
"%d",(
int)gINRATE));
4055 sline.ReplaceAll(
"xRATEANA",TString::Format(
"%d",(
int)gRATEANA));
4057 sline.ReplaceAll(
"xMEAN",
"Median");
4058 sline.ReplaceAll(
"xSIGMA",
"50% Percentile");
4059 sline.ReplaceAll(
"x2SIGMA",
"90% Percentile");
4061 sline.ReplaceAll(
"xMEAN",
"Mean");
4062 sline.ReplaceAll(
"xSIGMA",
"RMS");
4063 sline.ReplaceAll(
"x2SIGMA",
"2RMS");
4065 sline.ReplaceAll(
"HOME_CED_WWW",home_ced_www.Data());
4066 sline.ReplaceAll(
"CWB_DOC_URL",cwb_doc_url.Data());
4067 sline.ReplaceAll(
"HOME_WWW",home_www.Data());
4068 if(cwb_doc_url!=
"") {
4069 sline.ReplaceAll(
"<!--CWB_DOC_URL",
"");
4070 sline.ReplaceAll(
"CWB_DOC_URL-->",
"");
4071 sline.ReplaceAll(
"XCWB_DOC_URL",cwb_doc_url.Data());
4073 if(output) (*out) << sline.Data() << endl;
4084 n = ifo->
IWFP.size();
4085 for (
int i=0;
i<
n;
i++) {
4092 n = ifo->
RWFP.size();
4093 for (
int i=0;
i<
n;
i++) {
4112 TString xtitle = TString::Format(
"Time (sec) : GPS OFFSET = %.3f",gSEGGPS);
4118 for(
int i=0;
i<frec->
size();
i++) {
4120 if(xtime>bT && xtime<eT)
continue;
4127 TGraphErrors* e2gr =
new TGraphErrors(size,time.
data,favr->
data,etime.
data,f2err.
data);
4128 e2gr->SetLineColor(17);
4129 e2gr->SetFillStyle(1001);
4130 e2gr->SetFillColor(17);
4131 e2gr->GetXaxis()->SetTitle(xtitle);
4132 e2gr->GetYaxis()->SetTitle(
"frequency (hz)");
4133 e2gr->SetTitle(title);
4134 e2gr->GetXaxis()->SetTitleFont(42);
4135 e2gr->GetXaxis()->SetLabelFont(42);
4136 e2gr->GetXaxis()->SetLabelOffset(0.012);
4137 e2gr->GetXaxis()->SetTitleOffset(1.5);
4138 e2gr->GetYaxis()->SetTitleFont(42);
4139 e2gr->GetYaxis()->SetLabelFont(42);
4140 e2gr->GetYaxis()->SetLabelOffset(0.01);
4141 e2gr->GetYaxis()->SetTitleOffset(1.4);
4144 TGraphErrors* egr =
new TGraphErrors(size,time.
data,favr->
data,etime.
data,ferr->
data);
4145 egr->SetLineColor(15);
4146 egr->SetFillStyle(1001);
4147 egr->SetFillColor(15);
4149 TGraph* agr =
new TGraph(size,time.
data,favr->
data);
4150 agr->SetLineWidth(1);
4151 agr->SetLineColor(kWhite);
4152 agr->SetLineStyle(1);
4154 TGraph*
gr =
new TGraph(size,time.
data,frec->
data);
4155 gr->SetLineWidth(1);
4156 gr->SetLineColor(2);
4158 TCanvas*
canvas =
new TCanvas(
"frec",
"frec",200,20,800,500);
4163 e2gr->GetXaxis()->SetRangeUser(bT, eT);
4165 egr->Draw(
"P4same");
4171 canvas->Print(ofname);
4172 cout <<
"write : " << ofname << endl;
4193 TString xtitle = TString::Format(
"Time (sec) : GPS OFFSET = %.3f",gSEGGPS);
4195 TGraphErrors* e2gr =
new TGraphErrors(size,time.
data,wavr->
data,etime.
data,w2err.
data);
4196 e2gr->SetLineColor(17);
4197 e2gr->SetFillStyle(1001);
4198 e2gr->SetFillColor(17);
4199 e2gr->GetXaxis()->SetTitle(xtitle);
4200 e2gr->GetYaxis()->SetTitle(
"magnitude");
4201 e2gr->SetTitle(title);
4202 e2gr->GetXaxis()->SetTitleFont(42);
4203 e2gr->GetXaxis()->SetLabelFont(42);
4204 e2gr->GetXaxis()->SetLabelOffset(0.012);
4205 e2gr->GetXaxis()->SetTitleOffset(1.5);
4206 e2gr->GetYaxis()->SetTitleFont(42);
4207 e2gr->GetYaxis()->SetLabelFont(42);
4208 e2gr->GetYaxis()->SetLabelOffset(0.01);
4209 e2gr->GetYaxis()->SetTitleOffset(1.4);
4211 TGraphErrors* egr =
new TGraphErrors(size,time.
data,wavr->
data,etime.
data,werr->
data);
4212 egr->SetLineColor(15);
4213 egr->SetFillStyle(1001);
4214 egr->SetFillColor(15);
4216 TGraph* agr =
new TGraph(size,time.
data,wavr->
data);
4217 agr->SetLineWidth(1);
4218 agr->SetLineColor(kWhite);
4219 agr->SetLineStyle(1);
4221 TGraph*
gr =
new TGraph(size,time.
data,wrec->
data);
4222 gr->SetLineWidth(1);
4223 gr->SetLineColor(2);
4225 TCanvas*
canvas =
new TCanvas(
"wrec",
"wrec",200,20,800,500);
4234 e2gr->GetXaxis()->SetRangeUser(bT, eT);
4242 canvas->Print(ofname);
4243 cout <<
"write : " << ofname << endl;
4271 for(
int i=0;
i<wrec->
size();
i++) {
4272 if(time[
i]>bT && time[
i]<eT)
continue;
4277 TString xtitle = TString::Format(
"Time (sec) : GPS OFFSET = %.3f",gSEGGPS);
4279 TGraphAsymmErrors* egr90 =
new TGraphAsymmErrors(size,time.
data,wmed->
data,etime.
data,etime.
data,wl90->
data,wu90->
data);
4280 egr90->SetLineColor(17);
4281 egr90->SetFillStyle(1001);
4282 egr90->SetFillColor(17);
4283 egr90->GetXaxis()->SetTitle(xtitle);
4284 if(freq) egr90->GetYaxis()->SetTitle(
"frequency (hz)");
4285 else egr90->GetYaxis()->SetTitle(
"magnitude");
4286 egr90->SetTitle(title);
4287 egr90->GetXaxis()->SetTitleFont(42);
4288 egr90->GetXaxis()->SetLabelFont(42);
4289 egr90->GetXaxis()->SetLabelOffset(0.012);
4290 egr90->GetXaxis()->SetTitleOffset(1.5);
4291 egr90->GetYaxis()->SetTitleFont(42);
4292 egr90->GetYaxis()->SetLabelFont(42);
4293 egr90->GetYaxis()->SetLabelOffset(0.01);
4294 egr90->GetYaxis()->SetTitleOffset(1.4);
4296 TGraphAsymmErrors* egr50 =
new TGraphAsymmErrors(size,time.
data,wmed->
data,etime.
data,etime.
data,wl50->
data,wu50->
data);
4297 egr50->SetLineColor(15);
4298 egr50->SetFillStyle(1001);
4299 egr50->SetFillColor(15);
4301 TGraph* agr =
new TGraph(size,time.
data,wmed->
data);
4302 agr->SetLineWidth(1);
4303 agr->SetLineColor(kWhite);
4304 agr->SetLineStyle(1);
4306 TGraph*
gr =
new TGraph(size,time.
data,wrec->
data);
4307 gr->SetLineWidth(1);
4308 gr->SetLineColor(2);
4310 TCanvas*
canvas =
new TCanvas(
"wrec",
"wrec",200,20,800,500);
4315 egr90->GetXaxis()->SetRangeUser(bT, eT);
4317 egr50->Draw(
"3same");
4323 canvas->Print(ofname);
4324 cout <<
"write : " << ofname << endl;
4346 if (!out.good()) {cout <<
"Error Opening Output File : " << ofname << endl;
exit(1);}
4347 cout <<
"Create Output File : " << ofname << endl;
4351 out <<
"#whitened data : time," <<
4352 " amp_point, amp_mean, amp_rms," <<
4353 " amp_median, amp_lower_50_perc, amp_lower_90_perc, amp_upper_50_perc, amp_upper_90_perc," <<
4354 " frq_point, frq_mean, frq_rms," <<
4355 " frq_median, frq_lower_50_perc, frq_lower_90_perc, frq_upper_50_perc, frq_upper_90_perc" << endl;
4358 int size = vREC[
n].size();
4359 double dt=1./vREC[
n].rate();
4364 double time =
i*dt+vREC[
n].start();
4366 double vl50 = vMED[
n][
i]-
fabs(vL50[
n][
i]);
4367 double vu50 = vMED[
n][
i]+
fabs(vU50[
n][i]);
4368 double vl90 = vMED[
n][
i]-
fabs(vL90[
n][i]);
4369 double vu90 = vMED[
n][
i]+
fabs(vU90[
n][i]);
4371 double fl50 = fMED[
n][
i]-
fabs(fL50[
n][i]);
4372 double fu50 = fMED[
n][
i]+
fabs(fU50[
n][i]);
4373 double fl90 = fMED[
n][
i]-
fabs(fL90[
n][i]);
4374 double fu90 = fMED[
n][
i]+
fabs(fU90[
n][i]);
4377 <<
" " << vREC[
n][
i] <<
" " << vAVR[
n][
i] <<
" " << vRMS[
n][
i]
4378 <<
" " << vMED[
n][
i] <<
" " << vl50 <<
" " << vl90 <<
" " << vu50 <<
" " << vu90
4379 <<
" " << fREC[
n][
i] <<
" " << fAVR[
n][
i] <<
" " << fRMS[
n][
i]
4380 <<
" " << fMED[
n][
i] <<
" " << fl50 <<
" " << fl90 <<
" " << fu50 <<
" " << fu90
4400 if (!out.good()) {cout <<
"Error Opening Output File : " << ofname << endl;
exit(1);}
4401 cout <<
"Create Output File : " << ofname << endl;
4405 out <<
"# time white_amp" << endl;
4408 int size = vINJ[
n].size();
4409 double dt=1./vINJ[
n].rate();
4411 double time =
i*dt+vINJ[
n].start();
4412 out << time <<
" " << vINJ[
n][
i] << endl;
4422 int gLRETRY=-1;
IMPORT(
int,gLRETRY)
4426 if(gLRETRY==0)
return 1;
4428 bool use_original_data = (gLRETRY==1) ?
true :
false;
4430 if(gOPT.multitask) use_original_data = (gMTRIAL==gOPT.trials) ?
true :
false;
4432 if(use_original_data) cout << endl <<
"Last Trial : Analysis of the original data" << endl << endl;
4437 vector<TString> delObjList;
4439 delObjList.push_back(
"supercluster");
4440 delObjList.push_back(
"sparse");
4442 jname.ReplaceAll(
":/",
"");
4454 if(gOPT.signal==0) wf=vREC[
n];
4456 if(gOPT.signal==2) wf=vDAT[
n];
4459 if((!use_original_data) && gOPT.amp_cal_err[n]) {
4460 double amp_cal_err = 1;
4461 if(gOPT.amp_cal_err[n]>0) amp_cal_err = gRandom->Uniform(1-gOPT.amp_cal_err[n],1+gOPT.amp_cal_err[n]);
4462 else amp_cal_err = gRandom->Gaus(1,
fabs(gOPT.amp_cal_err[n]));
4463 cout <<
"Apply Amplitude Mis-Calibration (" << gOPT.amp_cal_err[
n] <<
"%) "
4464 << NET->
ifoName[
n] <<
" -> amp_cal_err : " << amp_cal_err << endl;
4468 if((!use_original_data) && gOPT.phs_cal_err[n]) {
4469 double phs_cal_err = 0;
4470 if(gOPT.phs_cal_err[n]>0) phs_cal_err = gRandom->Uniform(-gOPT.phs_cal_err[n],gOPT.phs_cal_err[n]);
4471 else phs_cal_err = gRandom->Gaus(0,
fabs(gOPT.phs_cal_err[n]));
4472 cout <<
"Apply Phase Mis-Calibration (" << gOPT.phs_cal_err[
n] <<
" deg) "
4473 << NET->
ifoName[
n] <<
" -> phs_cal_err : " << phs_cal_err <<
" deg" << endl;
4478 size = gHOT[
n].
size();
4479 rate = gHOT[
n].
rate();
4484 k = gRandom->Uniform(jS+jW,jE-jW);
4485 if(use_original_data) {
4486 *gCWB2G->
hot[
n] = gHOT[
n];
4489 for(
int i=jS;
i<jE;
i++) {
4509 jfile =
new TFile(jname);
4513 cout <<
"Loading sparse TF map ... " << endl;
4521 else sprintf(swname,
"sparse/%s-level:%d",ifo.Data(),
i+cfg->
l_low);
4524 cout <<
"CWB_Plugin_PE - Likelihood : sparse map " << swname
4525 <<
" not exist in job file" << endl;
exit(1);
4528 pD->
vSS.push_back(SS);
4538 vector<int> clist = pwc->
read(jfile,
"supercluster",
"clusters",0,cycle);
4539 return clist.size();
4546 if(wf1==NULL)
return wf;
4547 if(wf1->
size()==0)
return wf;
4549 double R=wf1->
rate();
4551 double b_wf1 = wf1->
start();
4552 double e_wf1 = wf1->
start()+wf1->
size()/R;
4553 double b_wf2 = wf2->
start();
4554 double e_wf2 = wf2->
start()+wf2->
size()/R;
4556 int o_wf1 = b_wf1>b_wf2 ? 0 :
int((b_wf2-b_wf1)*R+0.5);
4557 int o_wf2 = b_wf1<b_wf2 ? 0 :
int((b_wf1-b_wf2)*R+0.5);
4559 double startXCOR = b_wf1>b_wf2 ? b_wf1 : b_wf2;
4560 double endXCOR = e_wf1<e_wf2 ? e_wf1 : e_wf2;
4561 int sizeXCOR =
int((endXCOR-startXCOR)*R+0.5);
4563 for(
int i=0;
i<sizeXCOR;
i++) wf[
i+o_wf1] += wf2->
data[
i+o_wf2];
4575 for(
int n=0;
n<cfg->
nIFO;
n++) {
4589 cout <<
"Write file : " << cfg->
DQF[cfg->
nDQF-1].
file << endl;
4590 if (!out.good()) {cout <<
"Error Opening File : " << cfg->
DQF[cfg->
nDQF-1].
file << endl;
exit(1);}
4594 out <<
"1 " << istart <<
" " << istop <<
" " << 2*cfg->
iwindow << endl;
4605 int ntry=0;
for(
int i=0;
i<gOPT.trials;
i++)
if(wREC[
i].
size()) ntry++;
4607 int *
index =
new int[ntry];
4608 float *
value =
new float[ntry];
4610 int L = skyprob.
size();
4612 for(
int l=0;
l<
L;
l++) {
4614 int k=0;
for(
int i=0;
i<gOPT.trials;
i++)
if(wREC[
i].
size()) value[k++] = wSKYPROB[
i].
get(
l);
4615 TMath::Sort(ntry,value,index,
false);
4617 int imed = (ntry*50.)/100.;
if(imed>=ntry) imed=ntry-1;
4619 skyprob.
set(
l,value[index[imed]]);
4621 sum+=value[index[imed]];
4625 if(sum>0)
for(
int l=0;
l<
L;
l++) skyprob.
set(
l,skyprob.
get(
l)/sum);
4638 int L = skyprobcc.
size();
4639 for(
int l=0;
l<
L;
l++) {
4645 skyprobcc.
set(k,skyprob->
get(
l));
4651 sprintf(fname,
"%s/mskyprobcc.%s", odir.Data(),
"fits");
4655 if(NET->
MRA) analysis+=
":MRA";
4656 if(NET->
pattern==0) analysis+=
":Packet(0)";
4660 char configur[64]=
"";
4662 if (search==
'r')
sprintf(configur,
"%s un-modeled",analysis.Data());
4663 if (search==
'i')
sprintf(configur,
"%s iota-wave",analysis.Data());
4664 if (search==
'p')
sprintf(configur,
"%s psi-wave",analysis.Data());
4665 if((search==
'l')||(search==
's'))
sprintf(configur,
"%s linear",analysis.Data());
4666 if((search==
'c')||(search==
'g'))
sprintf(configur,
"%s circular",analysis.Data());
4667 if((search==
'e')||(search==
'b'))
sprintf(configur,
"%s elliptical",analysis.Data());
4677 char beginString[1024];
4679 char endString[1024];
4680 sprintf(endString,
"_job%d.root",runID);
4681 char containString[1024];
4687 for(
int i=0;i<nIFO;i++) mtpe_wREC[i] = new wavearray<double>;
4689 float*
chirp =
new float[6];
4690 float* isnr =
new float[
nIFO];
4691 float* osnr =
new float[
nIFO];
4692 float* iosnr =
new float[
nIFO];
4696 int nfile = fileList.size();
4697 for(
int n=0;
n<nfile;
n++) {
4698 cout <<
n <<
" " << fileList[
n].Data() << endl;
4700 TFile*
froot =
new TFile(fileList[
n].Data(),
"READ");
4702 cout <<
"CWB_Plugin_PE Error : Failed to open file : " << fileList[
n].Data() << endl;
4705 TTree*
itree = (TTree*)froot->Get(
"waveburst");
4707 cout <<
"CWB_Plugin_PE Error : Failed to open tree waveburst from file : " << fileList[
n].Data() << endl;
4711 for(
int i=0;
i<
nIFO;
i++) itree->SetBranchAddress(TString::Format(
"mtpe_wREC_%d",
i).Data(),&mtpe_wREC[
i]);
4712 itree->SetBranchAddress(
"mtpe_skyprob",&mtpe_skyprob);
4713 itree->SetBranchAddress(
"chirp",chirp);
4714 itree->SetBranchAddress(
"likelihood",&likelihood);
4720 if(mtpe_wREC[
i]==NULL) {
4721 cout <<
"CWB_Plugin_PE Error : Object wavearray not exist !!! " << endl;
4725 if(mtpe_skyprob==NULL) {
4726 cout <<
"CWB_Plugin_PE Error : Object skymap not exist !!! " << endl;
4731 wSKYPROB[gOPT.trials-
n-1] = *mtpe_skyprob;
4734 wREC[gOPT.trials-
n-1].push_back(*mtpe_wREC[
i]);
4736 for(
int i=0;i<
nIFO;i++) {
4737 iSNR[
i].push_back(isnr[i]);
4738 oSNR[
i].push_back(osnr[i]);
4739 ioSNR[
i].push_back(iosnr[i]);
4744 vLIKELIHOOD.push_back(likelihood);
4749 sprintf(command,
"/bin/rm %s",fileList[
n].Data());
4750 cout << command << endl;
4751 gSystem->Exec(command);
4754 for(
int i=0;
i<
nIFO;
i++)
delete mtpe_wREC[
i];
4755 delete mtpe_skyprob;
4782 int layers_high = 1<<cfg->
l_high;
4785 int index_min=2*gHOT[0].
size();
4786 for(
int k=0;
k<V;
k++) {
4790 if(index<index_min) index_min=
index;
4793 index_min = index_min-(index_min%layers_high);
4805 int index_off = gRandom->Uniform(jS,jE-jW);
4806 index_off = index_off-(index_off%layers_high);
4807 for(
int k=0;
k<V;
k++) {
4811 int ilayers = 1<<(ires+cfg->
l_low);
4812 index = index_off+(index-index_min);
4820 aNSE[
n].push_back(dd);
4821 ANSE[
n].push_back(DD);
4832 if(type!=
"data" && type!=
"null" && type!=
"signal" && type!=
"injection") {
4833 cout <<
"CWB_Plugin_PE Error : wrong PlotSpectrogram type, must be data/null/signal" << endl;
4847 if(type==
"signal") {
4852 if(type==
"injection") {
4858 *pTF*=1./sqrt(1<<cfg->
levelR);
4862 if(pTF->
size()==0)
continue;
4868 double fparm=nfact*6;
4873 int ysize=ystop-ystart;
4886 tstart+=0.9;tstop-=0.9;
4889 sprintf(fname,
"%s/%s_%s_spectrogram_0.png", pdir.Data(), NET->
ifoName[
n], type.Data());
4892 canvas->SetLogy(
true);
4894 sprintf(fname,
"%s/%s_%s_spectrogram_logy_0.png", pdir.Data(), NET->
ifoName[
n], type.Data());
4904 std::vector<wavearray<double> > xWHT;
4928 vector<int> clist = pwc->
read(jfile,
"supercluster",
"clusters",0,cycle);
4931 if(clist.size()==0)
return;
4937 for(
int k=0;
k<(
int)clist.size();
k++) {
4939 pwc->
read(jfile,
"supercluster",
"clusters",-2,cycle,0,clist[
k]);
4940 int psize = pwc->
size()-npixels;
4941 if(psize>msize) {maxk=
k;msize=psize;}
4949 for(
int k=0;
k<ctime.
size();
k++) {
4953 int kindex = gps>0 ? mink : maxk;
4957 pwc->
read(jfile,
"supercluster",
"clusters",-1,cycle,0,clist[kindex]);
4960 vector<TString> delObjList;
4961 delObjList.push_back(
"supercluster");
4963 jname.ReplaceAll(
":/",
"");
4968 jfile =
new TFile(jname,
"UPDATE");
4970 if(jfile==NULL||!jfile->IsOpen())
4971 {cout <<
"CWB_Plugin_PE.C - Error opening : " << jname << endl;
exit(1);}
4972 pwc->
write(jfile,
"supercluster",
"clusters",0,cycle);
4973 pwc->
write(jfile,
"supercluster",
"clusters",-1,cycle);
monster wdmMRA
list of pixel pointers for MRA
std::vector< char * > ifoName
std::vector< wavearray< double > > fL50
detector * getifo(size_t n)
param: detector index
static float _sse_abs_ps(__m128 *_a)
skymap wSKYPROB[MAX_TRIALS]
double FittingFactor(wavearray< double > *wf1, wavearray< double > *wf2)
virtual size_t size() const
std::vector< SSeries< double > > dSS[NIFO_MAX]
size_t write(const char *file, int app=0)
wavearray< double > GetWaveformEnvelope(wavearray< double > *wf)
void GetNullPixels(std::vector< double > *aNUL, std::vector< double > *ANUL, network *NET, CWB::config *cfg, int lag, int id)
Float_t * rho
biased null statistics
std::vector< wavearray< double > > fINJ
std::vector< wavearray< double > > vREC
Float_t * high
min frequency
std::vector< wavearray< double > > GetSigWaveform(network *NET, CWB::config *cfg, int lag, int id)
void FileGarbageCollector(TString ifName, TString ofName="", vector< TString > delObjList=vector< TString >())
std::vector< wavearray< double > > wREC[MAX_TRIALS]
std::vector< double > vCHIRP[6]
std::vector< netcluster > wc_List
void setSLags(float *slag)
void PlotWaveforms(network *NET, CWB::config *cfg, int ID, TString pdir="")
void SetSkyMask(network *NET, CWB::config *cfg, double theta, double phi, double radius)
printf("total live time: non-zero lags = %10.1f \n", liveTot)
void GetFactorsStatistic(int nIFO)
bool singleDetector
used for the stage stuff
void DumpRecWavePE(network *NET, TString pdir="")
Double_t * start
cluster duration = stopW-startW
WSeries< double > * pTF[NIFO_MAX]
pointers to detectors
std::vector< double > * getmdcTime()
WDM< double > * pwdm[NRES_MAX]
virtual void rate(double r)
void PlotChirpMass(int gtype, TString pdir, int sim=true)
std::vector< wavearray< float > > tdAmp
Float_t * low
average center_of_snr frequency
void Print(TString pname)
wavearray< double > a(hp.size())
std::vector< wavearray< double > > fU50
std::vector< wavearray< double > > fREC
std::vector< wavearray< double > > vDAT
wavearray< double > GetCutWaveform(wavearray< double > x, double bT, double eT, double bF, double eF)
wavearray< double > AddWaveforms(wavearray< double > *wf1, wavearray< double > *wf2)
wavearray< double > gHOT[NIFO_MAX]
static void _sse_zero_ps(__m128 *_p)
void PlotFactors(int gtype, int nIFO, TString pdir)
cout<< "skymap size : "<< L<< endl;for(int l=0;l< L;l++) sm.set(l, l);sm > const_cast< char * >("skymap.dat")
wavearray< float > a_90
buffer for cluster sky 00 amplitude
char * watversion(char c='s')
size_t setIndexMode(size_t=0)
void SetEventWindow(CWB::config *cfg, double gps)
TFile * jfile
output root file
std::vector< wavearray< double > > vRMS
Int_t * size
cluster volume
void ResetPCA(network *NET, CWB::config *cfg, netcluster *pwc, std::vector< netpixel > *vPIX, int ID)
std::vector< double > aNUL[NIFO_MAX]
static void PhaseShift(wavearray< double > &x, double pShift=0.)
void DumpSkyProb(skymap *skyprob, network *NET, netevent *&EVT, TString odir)
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
std::vector< pixdata > data
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...
std::vector< wavearray< double > * > RWFP
std::vector< wavearray< double > > GetInjWaveform(network *NET, netevent *EVT, int id, double factor)
int _sse_mra_ps(float *, float *, float, int)
void PlotWaveform(TString ofname, TString title, CWB::config *cfg, wavearray< double > *wf1, wavearray< double > *wf2, wavearray< double > *wf3, wavearray< double > *wref, bool fft=false, TString pdir="", double P=0.99, EColor col1=kRed, EColor col2=kBlue, EColor col3=(EColor) 418)
std::vector< wavearray< double > > vMED
std::vector< TGraph * > graph
std::vector< SSeries< double > > vSS[NIFO_MAX]
std::vector< wavearray< double > > cINJ
std::vector< wavearray< double > > vINJ
void Coherence(int ifactor)
double getTheta(size_t i)
Float_t * ioSNR
reconstructed snr waveform
std::vector< wavearray< double > > vNUL
TString DumpCED(network *NET, netevent *&EVT, CWB::config *cfg, double factor)
std::vector< double > fRES
size_t setcore(bool core, int id=0)
WSeries< double > waveBand
void SetOptions(int simulation, double rho, double inRate, bool useSparse=false, char *gtype=const_cast< char * >("png"), int paletteId=0)
std::vector< wavearray< double > > fL90
std::vector< double > oSNR[NIFO_MAX]
virtual void start(double s)
bool getMRAwave(size_t ID, size_t lag, char atype='S', int mode=0, bool tof=false)
std::vector< double > mdcTime
void setTDFilter(int nCoeffs, int L=1)
std::vector< double > ANUL[NIFO_MAX]
void plot(wavearray< double > &, char *=NULL, int=1, double=0., double=0., bool=false, float=0., float=0., bool=false, float=0., bool=false)
void Draw(int dpaletteId=0, Option_t *option="colfz")
wavearray< int > sparseIndex
std::vector< SSeries< double > > jSS[NIFO_MAX]
std::vector< detector * > ifoList
std::vector< double > ioSNR[NIFO_MAX]
void dopen(const char *fname, char *mode, bool header=true)
void Draw(double t1=0.0, double t2=0.0, double f1=0.0, double f2=0.0, double z1=0.0, double z2=0.0, int dpaletteId=DUMMY_PALETTE_ID, Option_t *option="colfz")
cout<< "Selected Pixels : "<< nPix<< endl;wc.cluster(1, 1);SSeries< double > ss
static float _avx_loadata_ps(float **p, float **q, float **u, float **v, float En, std::vector< float * > &pAVX, int I)
std::vector< wavearray< double > > vL90
size_t getSkyIndex(double th, double ph)
param: theta param: phi
std::vector< wavearray< double > > GetRecWaveform(network *NET, netevent *EVT, int id)
wavearray< float > pNRG
buffers for cluster residual energy
void PlotWaveformAsymmErrors(TString ofname, TString title, CWB::config *cfg, wavearray< double > *wrec, wavearray< double > *wmed, wavearray< double > *wl50, wavearray< double > *wu50, wavearray< double > *wl90, wavearray< double > *wu90, wavearray< double > *wref, TString pdir, double P, bool freq=false)
std::vector< wavearray< double > > vU50
void output2G(TTree *, network *, size_t, int, double)
Int_t Psave
number of detectors
std::vector< double > wRMS[NIFO_MAX]
void PlotSparse(int ifoID, network *NET, CWB::config *cfg, int ID, wavearray< double > *wf)
#define IMPORT(TYPE, VAR)
xtalk getXTalk(int nLay1, size_t indx1, int nLay2, size_t indx2)
param: numbers of layers identifying the resolution of the first pixel param: TF map index of the fir...
Double_t * gps
average center_of_gravity time
void SuperCluster(int ifactor)
wavearray< int > sparseLookup
double GetTimeBoundaries(wavearray< double > x, double P, double &bT, double &eT)
void setDelayIndex(double rate)
param: MDC log file
void PlotFrequencyErrors(TString ofname, TString title, CWB::config *cfg, wavearray< double > *frec, wavearray< double > *favr, wavearray< double > *ferr, wavearray< double > *wref, TString pdir, double P)
std::vector< netpixel > DoPCA(network *NET, CWB::config *cfg, int lag, int id)
Double_t * stop
GPS start time of the cluster.
wavearray< float > a_00
wdm multi-resolution analysis
std::vector< WDM< double > * > wdmList
WSeries< double > pTF[nRES]
static void _avx_free_ps(std::vector< float * > &v)
TIter next(twave->GetListOfBranches())
std::vector< wavearray< double > > fRMS
void GetNoisePixels(std::vector< double > *aNSE, std::vector< double > *ANSE, network *NET, CWB::config *cfg, int lag, int id)
int Write(double factor, bool fullCED=true)
std::vector< wavearray< double > > fAVR
void SetOutputFile(network *NET, netevent *&EVT, CWB::config *cfg, bool dump_pe)
static float _avx_setAMP_ps(float **p, float **q, std::vector< float * > &pAVX, int I)
char channelNamesRaw[NIFO_MAX][50]
void SetChannelName(char *chName)
std::vector< double > iSNR[NIFO_MAX]
Float_t * netcc
effective correlated SNR
void ClearWaveforms(detector *ifo)
std::vector< wavearray< double > * > IWFP
void PlotSNRnet(int nIFO, TString pdir, int sim=true)
std::vector< double > wAVR[NIFO_MAX]
void DumpOutputFile(network *NET, netevent *&EVT, CWB::config *cfg, int ID, int k, int factor)
void PlotResiduals(int ID, TString pdir="", int sim=true, char type='w')
size_t loadTDampSSE(network &net, char c, size_t BATCH=10000, size_t LOUD=0)
void DumpInjWavePE(network *NET, TString pdir="")
double GetFrequencyBoundaries(wavearray< double > x, double P, double &bF, double &eF)
skymap GetMedianSkyProb(network *NET)
netpixel * getPixel(size_t n, size_t i)
void WriteBodyHTML(TString html_template, TString html_tag_beg, TString html_tag_end, ofstream *out, int nIFO=1, TString *ifo=NULL)
std::vector< wavearray< double > > vWHT
std::vector< wavearray< double > > vAVR
void PrintUserOptions(CWB::config *cfg)
#define EXPORT(TYPE, VAR, CMD)
double GetCentralFrequency(wavearray< double > x)
std::vector< double > vLIKELIHOOD
Float_t * phi
sqrt(h+*h+ + hx*hx)
netevent EVT(itree, nifo)
double GetCentralTime(wavearray< double > x)
Float_t * theta
[0]-reconstructed, [1]-injected phi angle, [2]-RA
Double_t * time
beam pattern coefficients for hx
void Expand(bool bcore=true)
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *NET, WSeries< double > *x, TString ifo, int type)
COHERENCE.
void PlotWaveformErrors(TString ofname, TString title, CWB::config *cfg, wavearray< double > *wrec, wavearray< double > *wavr, wavearray< double > *werr, wavearray< double > *wref, TString pdir="", double P=0.99)
WSeries< double > waveForm
void Wave2Sparse(network *NET, CWB::config *cfg, char type)
wavearray< double > ts(N)
void GetChirpMassStatistic(std::vector< double > *vCHIRP)
std::vector< wavearray< double > > GetWhitenedData(network *NET, CWB::config *cfg)
WSeries< double > * getTFmap()
param: no parameters
cout<< "Starting reading output directory ..."<< endl;vector< TString > fileList
wavearray< float > _avx_norm_ps(float **, float **, std::vector< float * > &, int)
double fabs(const Complex &x)
double GetSparseMapData(SSeries< double > *SS, bool phase, int index)
void GetNullStatistic(std::vector< double > *vNUL, std::vector< double > *vNSE, int ifoId, TString ifoName, TString gtype, TString pdir="")
wavearray< float > sparseMap00
int SetSkyMask(network *net, CWB::config *cfg, char *options, char skycoord, double skyres=-1)
wavearray< float > sparseMap90
std::vector< wavearray< double > > vL50
std::vector< int > ComputeStatisticalErrors(network *NET, CWB::config *cfg, int ID)
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
strcpy(RunLabel, RUN_LABEL)
Meyer< double > S(1024, 2)
netcluster * getwc(size_t n)
param: delay index
void DumpUserOptions(TString odir, CWB::config *cfg)
std::vector< SSeries< double > > rSS[NIFO_MAX]
void PlotSpectrogram(TString type, network *NET, netevent *&EVT, CWB::config *cfg, TString pdir)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
void CreateIndexHTML(TString dirCED, int nIFO, TString *ifo, bool sim=false)
std::vector< clusterdata > cData
wavearray< double > GetDifWaveform(wavearray< double > *wf1, wavearray< double > *wf2)
void SaveSkyProb(network *NET, CWB::config *cfg, int id)
void set(size_t i, double a)
param: sky index param: value to set
std::vector< double > aNSE[NIFO_MAX]
std::vector< wavearray< double > > vPCA
std::vector< wavearray< double > > GetPCAWaveform(network *NET, CWB::config *cfg, int lag, int id)
int RedoAnalysis(TFile *jfile, CWB::config *cfg, network *NET)
wavearray< float > rNRG
buffer for cluster sky 90 amplitudes
double factors[FACTORS_MAX]
double get(size_t i)
param: sky index
WaveDWT< DataType_t > * pWavelet
skymap GetSkyProb(network *NET, int id)
wavearray< double > * hot[NIFO_MAX]
wavelet pointers: pwdm[0] - l_high, wdm[nRES-1] l_low
wavearray< double > wINJ[NIFO_MAX]
void ReplaceSuperclusterData(TFile *&jfile, CWB::config *cfg, network *NET, double gps=0)
SSeries< double > * getSTFmap(size_t n)
wavearray< double > GetAlignedWaveform(wavearray< double > *wf1, wavearray< double > *wf2)
Bool_t fill_in(network *, int, bool=true)
std::vector< SSeries< double > > vSS
std::vector< wavearray< double > > vU90
size_t read(const char *)
WSeries< double > waveNull
size_t getmdc__ID(size_t n)
std::vector< wavearray< double > > GetWaveform(network *NET, int lag, int id, char type, bool shift=true)
void AddNoiseAndCalErrToSparse(network *NET, CWB::config *cfg, char type)
wavearray< double > ** pwf
[x1,y1,z1,x2,y2,z2] components of spin vector
virtual void resize(unsigned int)
void Print(TString pname)
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
static float _avx_packet_ps(float **p, float **q, std::vector< float * > &pAVX, int I)
std::vector< wavearray< double > > fU90
void PlotTimeFrequencyPCA(network *NET, netevent *EVT, CWB::config *cfg, int id, int lag, TString pdir)
void SetOptions(TString projection="hammer", TString coordinate="Geographic", double resolution=1, bool goff=false)
std::vector< double > ANSE[NIFO_MAX]
Float_t * chirp
range to source: [0/1]-rec/inj
void GetSNRnetStatistic(int nIFO)
void Draw(GWAT_DRAW type=GWAT_TIME, TString options="ALP", Color_t color=kBlack)
void ReadUserOptions(TString options)
std::vector< double > vRES
void LoadFromMultiTaskJobOutput(int runID, CWB::config *cfg)
std::vector< wavearray< double > > fMED
void SetZaxisTitle(TString zAxisTitle)
int binarySearch(int array[], int start, int end, int key)
std::vector< vector_int > nTofF
bool setdata(double a, char type='R', size_t n=0)
Float_t * iSNR
energy of reconstructed responses Sk*Sk
void SetWaveformCuts(wavearray< double > *x, double bT, double eT, double bF, double eF)
Float_t * oSNR
injected snr waveform