21 #pragma GCC system_header 28 #include "TObjArray.h" 29 #include "TObjString.h" 54 #define PLOT_DIR "plot" 61 double& enINJ,
double& enREC,
double& xcorINJ_REC);
136 #define SKYMASK_SEED 0 138 #define APPLY_INJ_MRA 163 #define AMP_CAL_ERR 0.1 164 #define PHS_CAL_ERR 10 193 #define WRC_PLUGIN_VERSION "v1.3" 203 cout <<
"ifo " << ifo.Data() << endl;
204 cout <<
"type " << type << endl;
225 cout <<
"-----------------------------------------" << endl;
226 cout <<
"WRC config options " << endl << endl;
227 cout <<
"WRC_RETRY " <<
WRC_RETRY << endl;
229 cout <<
"WRC_NOISE " <<
WRC_NOISE << endl;
244 cout <<
"CWB_Plugin_xWRC.C -> CWB_PLUGIN_ILIKELIHOOD implemented only for 2G" << endl;
249 cout <<
"CWB_Plugin_xWRC.C -> implemented only for zero lag" << endl;
260 cout <<
"CWB_Plugin_xWRC.C -> CWB_PLUGIN_XLIKELIHOOD implemented only for 2G" << endl;
265 int gLRETRY=-1;
IMPORT(
int,gLRETRY)
266 cout <<
"-----> CWB_Plugin_xWRC.C -> " <<
" gLRETRY : " << gLRETRY <<
"/" <<
WRC_RETRY <<
" - WRC : ";
273 #ifdef USE_wINJ // use wINJ (default wREC) (are pixels used as inj+noise in the retry) 276 #ifdef USE_wWHT // use wWHT (default wREC) 279 #if !defined (USE_wINJ) && !defined (USE_wWHT) 288 if(gLRETRY==WRC_RETRY) {
303 if(gLRETRY==WRC_RETRY) {
322 cout <<
"CWB_Plugin_xWRC.C -> CWB_PLUGIN_OLIKELIHOOD implemented only for 2G" << endl;
327 int gLRETRY=-1;
IMPORT(
int,gLRETRY)
348 int ID = size_t(
id.data[
j]+0.5);
350 if(NET->
getwc(zlag)->
sCuts[ID-1]!=-1)
continue;
353 EXPORT(
int,gLRETRY,
"gLRETRY = 0;")
368 cout <<
"-----> Event Detected ID = " << ID << endl;
373 #if defined (SAVE_PLOT_WERR) || defined (SAVE_PLOT_WERR2) || (SAVE_WHT_PLOT) 389 for(
int j=0;
j<wAVR[
n].size();
j++) {
390 if(wTRY[
n][
j]==0)
continue;
391 wAVR[
n][
j]/=wTRY[
n][
j];
392 wRMS[
n][
j]/=wTRY[
n][
j];
393 wRMS[
n][
j]-=wAVR[
n][
j]*wAVR[
n][
j];
394 wRMS[
n][
j] =sqrt(wRMS[
n][
j]);
398 #ifdef SAVE_PLOT_WERR 402 #ifdef SAVE_PLOT_WERR2 409 #ifdef SAVE_MRA_PLOT // monster event display 410 PlotMRA(NET, ID, zlag,
"last_wREC");
418 #ifdef SAVE_MRA_PLOT // monster event display 419 PlotMRA(NET, ID, zlag,
"wREC");
473 double& enINJ,
double& enREC,
double& xcorINJ_REC) {
475 double bINJ = wfINJ->
start();
476 double eINJ = wfINJ->
stop();
477 double bREC = wfREC->
start();
478 double eREC = wfREC->
stop();
483 double R=wfINJ->
rate();
485 int oINJ = bINJ>bREC ? 0 :
int((bREC-bINJ)*R+0.5);
486 int oREC = bINJ<bREC ? 0 :
int((bINJ-bREC)*R+0.5);
489 double startXCOR = bINJ>bREC ? bINJ : bREC;
490 double endXCOR = eINJ<eREC ? eINJ : eREC;
491 int sizeXCOR =
int((endXCOR-startXCOR)*R+0.5);
498 if (sizeXCOR<=0)
return;
502 for (
int i=0;
i<sizeXCOR;
i++) enREC+=wfREC->
data[
i+oREC]*wfREC->
data[
i+oREC];
503 for (
int i=0;
i<sizeXCOR;
i++) xcorINJ_REC+=wfINJ->
data[
i+oINJ]*wfREC->
data[
i+oREC];
505 double erINJ_REC = enINJ+enREC-2*xcorINJ_REC;
515 double bINJ = wfINJ->
start();
516 double eINJ = wfINJ->
stop();
517 double bREC = wfREC->
start();
518 double eREC = wfREC->
stop();
523 double R=wfINJ->
rate();
525 int oINJ = bINJ>bREC ? 0 :
int((bREC-bINJ)*R+0.5);
526 int oREC = bINJ<bREC ? 0 :
int((bINJ-bREC)*R+0.5);
529 double startXCOR = bINJ>bREC ? bINJ : bREC;
530 double endXCOR = eINJ<eREC ? eINJ : eREC;
531 int sizeXCOR =
int((endXCOR-startXCOR)*R+0.5);
535 for (
int i=0;
i<sizeXCOR;
i++) {
549 std::vector<int> sCuts = NET->
getwc(k)->
sCuts;
561 cout <<
"write " << fname << endl;
562 WTS.
plot(pwc, 1, nIFO,
'L', 0, const_cast<char*>(
"COLZ"));
566 cout <<
"write " << fname << endl;
567 WTS.
plot(pwc, 1, nIFO,
'N', 0, const_cast<char*>(
"COLZ"));
578 watplot PTS(const_cast<char*>(
"ptswrc"),200,20,800,500);
585 PTS.
plot(wfINJ, const_cast<char*>(
"ALP"), 1, 0, 0,
true, cfg->
fLow, cfg->
fHigh);
587 PTS.
plot(wfINJ, const_cast<char*>(
"ALP"), 1, 0, 0);
589 PTS.
graph[0]->SetLineWidth(1);
591 PTS.
plot(wfREC, const_cast<char*>(
"SAME"), 2, 0, 0,
true, cfg->
fLow, cfg->
fHigh);
593 PTS.
plot(wfREC, const_cast<char*>(
"SAME"), 2, 0, 0);
595 PTS.
graph[1]->SetLineWidth(2);
600 if(fft)
sprintf(label,
"%s",
"fft");
601 else sprintf(label,
"%s",
"time");
602 if(strain)
sprintf(label,
"%s_%s",label,
"strain");
603 else sprintf(label,
"%s_%s",label,
"white");
608 cout <<
"write : " << fname << endl;
616 gRandom->SetSeed(seed);
643 if(wave_type!=
'j' && wave_type!=
'r' && wave_type!=
's' && wave_type!=
'n' && wave_type!=
'v') {
644 cout <<
"Wave2Sparse - Error : wrong wave type" << endl;
exit(1);
655 if(wave_type==
'v') vSS[
n]=pD->
vSS;
656 if(wave_type==
'r') rSS[
n]=pD->
vSS;
657 if(wave_type==
'j') jSS[
n]=pD->
vSS;
658 if(wave_type==
's') sSS[
n]=pD->
vSS;
660 if(wave_type==
'v')
return;
670 WF[
n].
rate(sSS[
n][0].wdm_rate);
671 WF[
n].
start(sSS[
n][0].wdm_start);
674 for (
int i=0;
i<sINJ[
n].
size();
i++) WF[
n][wos+
i] = sINJ[
n][
i];
678 WF[
n].
rate(jSS[
n][0].wdm_rate);
679 WF[
n].
start(jSS[
n][0].wdm_start);
682 for (
int i=0;
i<wINJ[
n].
size();
i++) WF[
n][wos+
i] = wINJ[
n][
i];
686 WF[
n].
rate(rSS[
n][0].wdm_rate);
687 WF[
n].
start(rSS[
n][0].wdm_start);
690 for (
int i=0;
i<wREC[
n].
size();
i++) WF[
n][wos+
i] = wREC[
n][
i];
698 for (
int i=0;
i<WF[
n].
size();
i++) WF[
n][
i] = gRandom->Gaus(0,1);
706 if(wave_type==
'r') gfile=TString::Format(
"%s/REC_%s.png",
PLOT_DIR,NET->
ifoName[n]);
707 if(wave_type==
'j') gfile=TString::Format(
"%s/INJ_%s.png",
PLOT_DIR,NET->
ifoName[n]);
708 if(wave_type==
's') gfile=TString::Format(
"%s/sINJ_%s.png",
PLOT_DIR,NET->
ifoName[n]);
709 if(wave_type!=
'n') (*plot) >>
gfile;
720 int layers = level>0 ? 1<<level : 0;
721 pwdm = NET->
getwdm(layers+1);
724 if((wave_type==
'n')&&(
i==nRES-1)) {
730 if(wave_type==
'r') size = rSS[
n][
i].sparseMap00.size();
731 if(wave_type==
'j') size = jSS[
n][
i].sparseMap00.size();
732 if(wave_type==
's') size = sSS[
n][
i].sparseMap00.size();
733 if(wave_type==
'n') vN[
n].push_back(w);
737 if(wave_type==
'r') index = rSS[
n][
i].sparseIndex[
j];
738 if(wave_type==
'j') index = jSS[
n][
i].sparseIndex[
j];
739 if(wave_type==
's') index = sSS[
n][
i].sparseIndex[
j];
743 rSS[
n][
i].sparseMap00[
j]=v00;
744 rSS[
n][
i].sparseMap90[
j]=v90;
747 jSS[
n][
i].sparseMap00[
j]=v00;
748 jSS[
n][
i].sparseMap90[
j]=v90;
751 sSS[
n][
i].sparseMap00[
j]=v00;
752 sSS[
n][
i].sparseMap90[
j]=v90;
763 gRandom->SetSeed(seed);
765 TF2 *gaus2d =
new TF2(
"gaus2d",
"exp(-x*x/2)*exp(-y*y/2)", -5,5,-5,5);
774 #ifdef USE_wWHT // use wWHT (default wREC) 802 int size = pD->
vSS[
i].sparseMap00.size();
809 double aa = pD->
vSS[
i].sparseMap00[
j];
810 double AA = pD->
vSS[
i].sparseMap90[
j];
812 pD->
vSS[
i].sparseMap00[
j] = C*aa + S*AA;
813 pD->
vSS[
i].sparseMap90[
j] = -S*aa + C*AA ;
822 int levels = (1<<(nRES-1-
i));
827 itime = levels*sI[itime/levels];
828 index = itime*layers +
ifreq;
830 r00 = vN[
n][
i].pWavelet->pWWS[
index];
831 r90 = vN[
n][
i].pWavelet->pWWS[index+vN[
n][
i].maxIndex()+1];
834 pD->
vSS[
i].sparseMap00[
j]+=r00;
835 pD->
vSS[
i].sparseMap90[
j]+=r90;
844 std::vector<float>* vP;
845 std::vector<int>* vI;
848 if(NET->
wc_List[0].p_Map.size()) {
850 vP = &(NET->
wc_List[0].p_Map[
id-1]);
851 vI = &(NET->
wc_List[0].p_Ind[
id-1]);
856 for(
int j=0;
j<
int(vP->size());
j++) {
858 double th = rec_skyprob.
getTheta(i);
861 rec_skyprob.
set(k,(*vP)[
j]);
868 for(
int l=0;
l<rec_skyprob.
size();
l++) prob+=rec_skyprob.
get(
l);
869 cout <<
"PROB : " << prob << endl;
891 for(
int j=0;
j<geINJ.
size();
j++) geINJ[
j]=sqrt(wINJ[n][
j]*wINJ[n][
j]+geINJ[
j]*geINJ[
j]);
895 for(
int j=0;j<wAVR[
n].size();j++) geAVR[j]=wAVR[n][j];
897 for(
int j=0;j<geAVR.
size();j++) geAVR[j]=sqrt(wAVR[n][j]*wAVR[n][j]+geAVR[j]*geAVR[j]);
901 for(
int j=0;j<wRMS[
n].size();j++) geRMS[j]=wRMS[n][j];
916 for(
int j=1;j<gF.
size();j++) {
920 gF[
j]=atan2(XX[j],xx[j])-atan2(XX[j-1],xx[j-1]);
921 if(gF[j]<-Pi) gF[
j]+=2*
Pi;
922 if(gF[j]> Pi) gF[
j]-=2*
Pi;
928 for(
int j=0;j<yy.
size();j++) yy[j]=wAVR[n][j];
932 for(
int j=1;j<gFF.
size();j++) {
933 gFF[
j]=atan2(YY[j],yy[j])-atan2(YY[j-1],yy[j-1]);
934 if(gFF[j]<-Pi) gFF[
j]+=2*
Pi;
935 if(gFF[j]> Pi) gFF[
j]-=2*
Pi;
941 gfile=TString::Format(
"%s/FREQ_%d.root",
PLOT_DIR,n);
947 for(
int j=0;j<geDIF.
size();j++) geDIF[j]=geAVR[j]-geINJ[j];
952 gfile=TString::Format(
"%s/eDIF_%d.root",
PLOT_DIR,n);
959 for(
int j=0;j<wAVR[
n].size();j++) gwAVR[j]=wAVR[n][j];
961 for(
int j=0;j<wAVR[
n].size();j++) gwAVR[j]=wRMS[n][j];
965 gfile=TString::Format(
"%s/wAVR_%d.root",
PLOT_DIR,n);
970 double bINJ = wINJ[
n].
start();
971 double eINJ = wINJ[
n].
stop();
972 double bREC = wREC[
n].
start();
973 double eREC = wREC[
n].
stop();
975 cout <<
"bINJ : " << bINJ <<
" eINJ : " << eINJ << endl;
976 cout <<
"bREC : " << bREC <<
" eREC : " << eREC << endl;
978 double R=wINJ[
n].
rate();
980 int oINJ = bINJ>bREC ? 0 :
int((bREC-bINJ)*R+0.5);
981 int oREC = bINJ<bREC ? 0 :
int((bINJ-bREC)*R+0.5);
982 cout <<
"oINJ : " << oINJ <<
" oREC : " << oREC << endl;
984 double startXCOR = bINJ>bREC ? bINJ : bREC;
985 double endXCOR = eINJ<eREC ? eINJ : eREC;
986 int sizeXCOR =
int((endXCOR-startXCOR)*R+0.5);
987 cout <<
"sizeXCOR : " << sizeXCOR << endl;
988 cout <<
"wINJ[n].size() : " << wINJ[
n].
size() << endl;
989 cout <<
"wREC[n].size() : " << wREC[
n].
size() << endl;
992 for(
int j=0;j<sizeXCOR;j++) gwREC[j+oINJ]=wREC[n][j+oREC];
995 for(
int j=0;j<wINJ[
n].
size();j++) gwREC[j]=wINJ[n][j];
997 for(
int j=0;j<wAVR[
n].size();j++) gwREC[j]=wRMS[n][j];
1001 gfile=TString::Format(
"%s/wREC_%d.root",
PLOT_DIR,n);
1034 cout <<
"tMin " << tMin <<
" tMax " << tMax << endl;
1038 for(
int j=0;
j<wAVR[
n].size();
j++) gw[
j]=wAVR[
n][
j];
1040 for(
int j=0;
j<wAVR[
n].size();
j++) gw[
j]=wRMS[
n][
j];
1046 sprintf(gtitle,
"%s : %10.3f (gps) - injected(Black) - average(Red) - rms(Blue)",gw.
start(),NET->
ifoName[
n]);
1047 plot->
gtitle(gtitle,
"time(sec)",
"amplitude");
1066 cout <<
"tMin " << tMin <<
" tMax " << tMax << endl;
1076 plot->
gtitle(gtitle,
"time(sec)",
"amplitude");
1094 for(
int i=0;
i<vREC[
n].size();
i++) {
1117 double bINJ = wINJ[
n].
start();
1118 double eINJ = wINJ[
n].
stop();
1119 double bREC = wREC[
n].
start();
1120 double eREC = wREC[
n].
stop();
1122 double R=wINJ[0].
rate();
1124 int oINJ = bINJ>bREC ? 0 :
int((bREC-bINJ)*R+0.5);
1125 int oREC = bINJ<bREC ? 0 :
int((bINJ-bREC)*R+0.5);
1127 double startXCOR = bINJ>bREC ? bINJ : bREC;
1128 double endXCOR = eINJ<eREC ? eINJ : eREC;
1129 int sizeXCOR =
int((endXCOR-startXCOR)*R+0.5);
1131 if (sizeXCOR<=0)
continue;
1133 for (
int j=0;
j<sizeXCOR;
j++) winj[
n][
j+oREC]+=wINJ[
n].data[
j+oINJ];
1142 double bINJ = wREC[
n].
start();
1143 double eINJ = wREC[
n].
stop();
1144 for(
int i=0;
i<vREC[
n].size();
i++) {
1146 double bREC = vREC[
n][
i].start();
1147 double eREC = vREC[
n][
i].stop();
1152 double R=wREC[
n].
rate();
1154 int oINJ = bINJ>bREC ? 0 :
int((bREC-bINJ)*R+0.5);
1155 int oREC = bINJ<bREC ? 0 :
int((bINJ-bREC)*R+0.5);
1158 double startXCOR = bINJ>bREC ? bINJ : bREC;
1159 double endXCOR = eINJ<eREC ? eINJ : eREC;
1160 int sizeXCOR =
int((endXCOR-startXCOR)*R+0.5);
1164 if (sizeXCOR<=0)
continue;
1166 for (
int j=0;
j<sizeXCOR;
j++) wavr[
n][
j+oINJ]+=vREC[
n][
i].data[
j+oREC];
1168 wavr[
n]*=1./vREC[
n].size();
1173 for(
int i=0;
i<vREC[0].size();
i++) {
1179 rnum+=(enINJ[
n]+enREC[
n]-2*xcorINJ_REC[
n]);
1193 rnum+=(enINJ[
n]+enREC[
n]-2*xcorINJ_REC[
n]);
1203 rnum+=(enINJ[
n]+enREC[
n]-2*xcorINJ_REC[
n]);
1215 rnum+=(enINJ[
n]+enREC[
n]-2*xcorINJ_REC[
n]);
1220 #ifdef SAVE_WF_COMPARISON 1226 PlotFinal3(NET, ID,
n, &rINJ[
n], &winj[n],
"rINJvsINJ",
"root");
1227 PlotFinal3(NET, ID, n, &rINJ[n], &wREC[n],
"rINJvsREC",
"root");
1228 PlotFinal3(NET, ID, n, &rINJ[n], &wavr[n],
"rINJvsAVR",
"root");
1230 PlotFinal3(NET, ID, n, &winj[n], &wavr[n],
"INJvsAVR",
"root");
1231 PlotFinal3(NET, ID, n, &wREC[n], &wavr[n],
"RECvsAVR",
"root");
1240 for(
int i=1;
i<10;
i++) erR[
i]=vRES[
index[
i*
int(vRES.
size()/10.)]];
1243 for(
int i=0;
i<vRES.
size();
i++) {
1244 if(vRES[
i]<=jres) erR[10]+=1.;
1246 erR[10]/=vRES.
size();
1250 for(
int i=0;
i<10;
i++) cout << erR[
i] <<
"/";
1251 cout << erR[10] << endl;
1253 #ifdef SAVE_PLOT_RESIDUALS 1256 TCanvas* canvas2=
new TCanvas(
"canvas2",
"canvas2", 200, 20, 600, 600);
1257 TH1F* hres =
new TH1F(
"residuals",
"residuals",10,0,1);
1258 for(
int i=0;
i<vRES.
size();
i++) hres->Fill(vRES[
i]);
1263 float ymax = hres->GetMaximum();
1264 TLine *
line =
new TLine(jres,0,jres,ymax);
1265 line->SetLineColor(kRed);
1268 gStyle->SetLineColor(kWhite);
1269 hres->SetTitle(
"Distribution of residuals");
1270 hres->SetStats(kTRUE);
1271 canvas2->Print(TString::Format(
"%s/residuals_%d.png",
PLOT_DIR,ID));
1288 while(!vSS[
n].empty()) vSS[
n].pop_back();
1289 vSS[
n].clear(); std::vector<SSeries<double> >().swap(vSS[
n]);
1290 while(!sSS[n].empty()) sSS[
n].pop_back();
1291 sSS[
n].clear(); std::vector<SSeries<double> >().swap(sSS[n]);
1292 while(!rSS[n].empty()) rSS[
n].pop_back();
1293 rSS[
n].clear(); std::vector<SSeries<double> >().swap(rSS[n]);
1294 while(!jSS[n].empty()) jSS[
n].pop_back();
1295 jSS[
n].clear(); std::vector<SSeries<double> >().swap(jSS[n]);
1297 while(!vN[n].empty()) vN[
n].pop_back();
1298 vN[
n].clear(); std::vector<WSeries<double> >().swap(vN[n]);
1300 while(!wTRY[n].empty()) wTRY[
n].pop_back();
1301 wTRY[
n].clear(); std::vector<int >().swap(wTRY[n]);
1303 while(!wAVR[n].empty()) wAVR[
n].pop_back();
1304 wAVR[
n].clear(); std::vector<double >().swap(wAVR[n]);
1306 while(!wRMS[n].empty()) wRMS[
n].pop_back();
1307 wRMS[
n].clear(); std::vector<double >().swap(wRMS[n]);
1309 while(!vDREC[n].empty()) vDREC[
n].pop_back();
1310 vDREC[
n].clear(); std::vector<double >().swap(vDREC[n]);
1312 while(!vREC[n].empty()) vREC[
n].pop_back();
1313 vREC[
n].clear(); std::vector<wavearray<double> >().swap(vREC[n]);
1320 TString cwb_inet_options=
TString(gSystem->Getenv(
"CWB_INET_OPTIONS"));
1321 if(cwb_inet_options.CompareTo(
"")!=0) {
1322 TString inet_options = cwb_inet_options;
1323 cout << inet_options << endl;
1324 if(!inet_options.Contains(
"--")) {
1327 for(
int j=0;
j<token->GetEntries();
j++){
1329 TObjString* tok = (TObjString*)token->At(
j);
1330 TString stok = tok->GetString();
1332 if(stok.Contains(
"wrc_id=")) {
1334 wrc_id.Remove(0,wrc_id.Last(
'=')+1);
1335 if(wrc_id.IsDigit())
WRC_ID=wrc_id.Atoi();
1337 cout <<
"WRC_ID : " <<
WRC_ID << endl;
1340 if(stok.Contains(
"wrc_retry=")) {
1342 wrc_retry.Remove(0,wrc_retry.Last(
'=')+1);
1343 if(wrc_retry.IsDigit())
WRC_RETRY=wrc_retry.Atoi();
1345 cout <<
"WRC_RETRY : " <<
WRC_RETRY << endl;
1348 if(stok.Contains(
"wrc_ced_id=")) {
1350 wrc_ced_id.Remove(0,wrc_ced_id.Last(
'=')+1);
1351 if(wrc_ced_id.IsDigit())
WRC_CED_ID=wrc_ced_id.Atoi();
1353 cout <<
"WRC_CED_ID : " <<
WRC_CED_ID << endl;
1356 if(stok.Contains(
"wrc_ced_retry=")) {
1358 wrc_ced_retry.Remove(0,wrc_ced_retry.Last(
'=')+1);
1359 if(wrc_ced_retry.IsDigit())
WRC_CED_RETRY=wrc_ced_retry.Atoi();
1364 if(stok.Contains(
"wrc_skymask=")) {
1366 wrc_skymask.Remove(0,wrc_skymask.Last(
'=')+1);
1372 if(stok.Contains(
"wrc_noise=")) {
1374 wrc_noise.Remove(0,wrc_noise.Last(
'=')+1);
1406 __m128* _m00 = (__m128*) mam;
1407 __m128* _m90 = (__m128*) mAM;
1408 __m128* _amp = (__m128*) amp;
1409 __m128* _AMP = (__m128*) AMP;
1410 __m128* _a00 = (__m128*) NET->
a_00.
data;
1411 __m128* _a90 = (__m128*) NET->
a_90.
data;
1415 for(j=0; j<V; ++j) if(ee[j]>ee[
m]) m=j;
1416 if(ee[m]<=Eo)
break; mm = m*
f;
1427 for(j=0; j<J; j++) {
1433 if(ee[n]>0) cc += c[5];
1449 for(j=0; j<J; j++) {
1451 if(E<E2 && n!=m) {c+=7;
continue;}
1469 float* gSLAGSHIFT=NULL;
IMPORT(
float*,gSLAGSHIFT)
1474 TFile*
froot = NULL;
1475 TList *files = (TList*)gROOT->GetListOfFiles();
1482 while ((file=(TSystemFile*)
next())) {
1483 fname = file->GetName();
1485 if(fname.Contains(
"wave_")) {
1486 froot=(TFile*)file;froot->cd();
1488 outDump.ReplaceAll(
".root.tmp",
".txt");
1493 cout <<
"CWB_Plugin_xWRC.C : Error - output root file not found" << endl;
1497 cout <<
"CWB_Plugin_xWRC.C : Error - output root file not found" << endl;
1501 net_tree = (TTree *) froot->Get(
"waveburst");
1511 net_tree->Branch(
"erR",
erR,TString::Format(
"erR[%i]/F",11));
1529 int gLRETRY=-1;
IMPORT(
int,gLRETRY)
1533 EVT->
output2G(NULL,NET,ID,k,factor);
1535 cout <<
"rec_theta : " << EVT->
theta[0] <<
" rec_phi : " << EVT->
phi[0] << endl;
1541 int idSize = pd->
RWFID.size();
1543 for(
int mm=0; mm<idSize; mm++)
if (pd->
RWFID[mm]==ID) wfIndex=mm;
1546 cout <<
"CWB_Plugin_xWRC.C : Error : Reconstructed waveform not saved !!! : ID -> " 1547 << ID <<
" : detector " << NET->
ifoName[
n] << endl;
1550 if(wfIndex>=0) pwfREC[
n] = pd->
RWFP[wfIndex];
1553 double delay = EVT->
time[
n]-EVT->
time[0];
1558 vREC[
n].push_back(*wfREC);
1559 vDREC[
n].push_back(delay);
1588 int gLRETRY=-1;
IMPORT(
int,gLRETRY)
1592 double recTime = EVT->
time[0];
1597 double injTime = 1.e12;
1599 for(
int m=0;
m<
M;
m++) {
1602 if(T<injTime && INJ.
fill_in(NET,mdcID)) {
1617 pwfINJ[
n] = INJ.
pwf[
n];
1618 if (pwfINJ[
n]==NULL) {
1619 cout <<
"CWB_Plugin_xWRC.C : Error : Injected waveform not saved !!! : detector " 1624 double rFactor = 1.;
1625 rFactor *= factor>0 ?
factor : 1;
1633 double delay = EVT->
time[
n]-EVT->
time[0];
1637 wTRY[
n].resize(wINJ[n].
size());
1638 wAVR[
n].resize(wINJ[n].
size());
1639 wRMS[
n].resize(wINJ[n].
size());
1640 for(
int j=0;
j<wINJ[
n].
size();
j++) {wTRY[
n][
j]=0;wAVR[
n][
j]=0;wRMS[
n][
j]=0;};
1658 std::vector<int>* vint = &(pwc->
cList[ID-1]);
1659 int V = vint->size();
1665 irate = pv!=NULL ? (*pv)[0] : 0;
1667 bool isPCs = !(cfg->
optim&&std::isupper(cfg->
search));
1677 for(
int n=0;
n<V;
n++) {
1680 cout <<
"CWB_Plugin_eWRC.C - Error : is enabled only for WDM (2G)" << endl;
1684 #ifdef APPLY_INJ_MRA 1691 #ifdef APPLY_INJ_MRA 1694 if(!pix->
core)
continue;
1695 if((irate)&&(irate !=
int(pix->
rate+0.5)))
continue;
1703 #ifdef APPLY_INJ_MRA 1706 NET->
rNRG[
n]+=aa*aa+AA*AA;
1712 #ifdef SAVE_MRA_PLOT // monster event display 1716 #ifdef APPLY_INJ_MRA 1718 int V4 = V + (V%4 ? 4 - V%4 : 0);
1727 for(
int n=0;
n<V;
n++) {
1730 cout <<
"CWB_Plugin_eWRC.C - Error : is enabled only for WDM (2G)" << endl;
1737 if(!pix->
core)
continue;
1738 if((irate)&&(irate !=
int(pix->
rate+0.5)))
continue;
1740 double aa = amp[
n*
NIFO+
m];
1741 double AA = AMP[
n*
NIFO+
m];
1747 #ifdef SAVE_MRA_PLOT // monster event display 1748 PlotMRA(NET, ID, k,
"mra_rINJ");
monster wdmMRA
list of pixel pointers for MRA
wavearray< double > gwREC[NIFO_MAX]
std::vector< char * > ifoName
detector * getifo(size_t n)
param: detector index
static float _sse_abs_ps(__m128 *_a)
std::vector< int > vector_int
void AddNoise2Sparse(network *NET, CWB::config *cfg, int seed=0)
std::vector< int > wTRY[NIFO_MAX]
std::vector< vector_int > cRate
void gtitle(TString title="", TString xtitle="", TString ytitle="")
void ComputeErrorWF(wavearray< double > *wfINJ, wavearray< double > *wfREC, int ifoID)
void PlotWaveform(TString ifo, wavearray< double > *wfINJ, wavearray< double > *wfREC, CWB::config *cfg, bool fft=false, bool strain=false)
std::vector< netcluster > wc_List
void SaveSkyProb(network *NET, CWB::config *cfg, int id)
void setSLags(float *slag)
double GetSparseMap(SSeries< double > *SS, bool phase, int index)
std::vector< WSeries< double > > vN[NIFO_MAX]
void PlotFinal2(network *NET, int ID)
std::vector< double > * getmdcTime()
virtual void rate(double r)
static void _sse_zero_ps(__m128 *_p)
wavearray< float > a_90
buffer for cluster sky 00 amplitude
static void PhaseShift(wavearray< double > &x, double pShift=0.)
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<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ 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< vector_float > sArea
wavearray< double > rINJ[NIFO_MAX]
void Wave2Sparse(network *NET, CWB::config *cfg, char wave_type)
std::vector< TGraph * > graph
plot gtitle(gtitle,"frequency (Hz)","strain/#sqrt{Hz}")
static float _sse_rotsub_ps(__m128 *_u, float c, __m128 *_v, float s, __m128 *_a)
double getTheta(size_t i)
void DumpOutput(network *NET, netevent *&EVT, CWB::config *cfg, int ID, int k, int factor)
void Draw(int dpaletteId=1, Option_t *option="colfz")
static void _sse_cpf_ps(float *a, __m128 *_p)
std::vector< vector_int > cList
virtual void start(double s)
bool getMRAwave(size_t ID, size_t lag, char atype='S', int mode=0, bool tof=false)
void plot(wavearray< double > &, char *=NULL, int=1, double=0., double=0., bool=false, float=0., float=0., bool=false, float=0., bool=false)
wavearray< int > sparseIndex
static void _sse_rotadd_ps(__m128 *_u, float c, __m128 *_v, float s, __m128 *_a)
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
void dopen(const char *fname, char *mode, bool header=true)
WDM< double > * getwdm(size_t M)
param: number of wdm layers
size_t getSkyIndex(double th, double ph)
param: theta param: phi
virtual size_t size() const
void output2G(TTree *, network *, size_t, int, double)
void GetInjWaveform(network *NET, netevent *EVT, int ID, int k, int factor)
Int_t Psave
number of detectors
std::vector< double > wRMS[NIFO_MAX]
#define IMPORT(TYPE, VAR)
bool WRC_PHS_CAL_ERR_CONFIG
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...
wavearray< int > sparseLookup
std::vector< SSeries< double > > rSS[NIFO_MAX]
std::vector< SSeries< double > > sSS[NIFO_MAX]
wavearray< double > wREC[NIFO_MAX]
void ComputeErrorStatistic(network *NET, CWB::config *cfg, int ID)
wavearray< float > a_00
wdm multi-resolution analysis
void PlotFinal3(network *NET, int ID, int ifoID, wavearray< double > *w1, wavearray< double > *w2, TString tag, TString gtype="png")
TIter next(twave->GetListOfBranches())
void SetOutput(network *NET, netevent *&EVT, CWB::config *cfg)
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *NET, WSeries< double > *x, TString ifo, int type)
COHERENCE.
bool WRC_SKYMASK_REC_SKY_CONFIG
double GetTimeRange(double &tMin, double &tMax, double efraction=0.9999999)
int _sse_mra_ps(network *NET, float *amp, float *AMP, float Eo, int K)
size_t cpf(const netcluster &, bool=false, int=0)
std::vector< double > wAVR[NIFO_MAX]
void SetTimeRange(double tMin, double tMax)
wavearray< double > sINJ[NIFO_MAX]
std::vector< SSeries< double > > vSS[NIFO_MAX]
netpixel * getPixel(size_t n, size_t i)
static void TimeShift(wavearray< double > &x, double tShift=0.)
#define EXPORT(TYPE, VAR, CMD)
void GetRecWaveform(network *NET, netevent *EVT, int ID, int k, int factor)
void ComputeResidualEnergy(wavearray< double > *wfINJ, wavearray< double > *wfREC, double &enINJ, double &enREC, double &xcorINJ_REC)
Float_t * phi
sqrt(h+*h+ + hx*hx)
void PlotMRA(network *NET, int ID, int k, TString tag)
netevent EVT(itree, nifo)
Float_t * theta
[0]-reconstructed, [1]-injected phi angle, [2]-RA
Double_t * time
beam pattern coefficients for hx
static void _sse_add_ps(__m128 *_a, __m128 *_b)
WSeries< double > waveForm
bool WRC_AMP_CAL_ERR_CONFIG
virtual void stop(double s)
double fabs(const Complex &x)
void SetSkyMask(network *NET, CWB::config *cfg, int seed=0)
wavearray< float > sparseMap00
int SetSkyMask(network *net, CWB::config *cfg, char *options, char skycoord, double skyres=-1)
wavearray< float > sparseMap90
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Meyer< double > S(1024, 2)
netcluster * getwc(size_t n)
param: delay index
std::vector< double > vDREC[NIFO_MAX]
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
void set(size_t i, double a)
param: sky index param: value to set
static void _sse_mul_ps(__m128 *_a, float b)
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
wavearray< double > wINJ[NIFO_MAX]
Bool_t fill_in(network *, int, bool=true)
std::vector< SSeries< double > > vSS
size_t getmdc__ID(size_t n)
for(int i=0;i< 101;++i) Cos2[2][i]=0
pointers to detectors
wavearray< double > ** pwf
[x1,y1,z1,x2,y2,z2] components of spin vector
virtual void resize(unsigned int)
void Print(TString pname)
#define WRC_PLUGIN_VERSION
std::vector< wavearray< double > > vREC[NIFO_MAX]
void Draw(GWAT_DRAW type=GWAT_TIME, TString options="ALP", Color_t color=kBlack)
void GetRInjWaveform(network *NET, netevent *EVT, CWB::config *cfg, int ID, int k)
std::vector< double > vRES
void SetZaxisTitle(TString zAxisTitle)
int binarySearch(int array[], int start, int end, int key)
std::vector< SSeries< double > > jSS[NIFO_MAX]
bool setdata(double a, char type='R', size_t n=0)