28 #define MIN_SKYRES_HEALPIX 4 29 #define MIN_SKYRES_ANGLE 3 32 #define REGRESSION_FILTER_LENGTH 8 33 #define REGRESSION_MATRIX_FRACTION 0.95 34 #define REGRESSION_SOLVE_EIGEN_THR 0. 35 #define REGRESSION_SOLVE_EIGEN_NUM 10 36 #define REGRESSION_SOLVE_REGULATOR 'h' 37 #define REGRESSION_APPLY_THR 0.8 40 #define SELECT_SUBRHO 5.0 41 #define SELECT_SUBNET 0.1 44 #define WDM_BETAORDER 6 // beta function order for Meyer 45 #define WDM_PRECISION 10 // wavelet precision 47 #define EXIT(ERR) gSystem->Exit(ERR) // better exit handling for ROOT stuff 57 if(pwdm[
i]!=NULL)
delete pwdm[
i];
69 sprintf(cmd,
"gCWB2G = (void*)%p;",
this);
EXPORT(
void*,gCWB2G,cmd);
73 {cout <<
"cwb2G::Init - Error : analysis must be 2G" << endl;
EXIT(1);}
75 {cout <<
"cwb2G::Init - Error : if initial stage!=CWB_STAGE_FULL " 76 <<
"then input file name must be a job file " << endl;
EXIT(1);}
84 char MRAcatalog[1024];
86 cout <<
"cwb2G::Init - Loading catalog of WDM cross-talk coefficients ... " << endl;
87 cout << MRAcatalog << endl;
104 int layers = level>0 ? 1<<level : 0;
110 cout <<
"cwb2G::Init : Error - filter length must be <= segEdge !!!" << endl;
111 cout <<
"filter length : " << wdmFLen <<
" sec" << endl;
112 cout <<
"cwb scratch : " <<
cfg.
segEdge <<
" sec" << endl;
115 cout <<
"Filter length = " << wdmFLen <<
" (sec)" << endl;
122 cout <<
"cwb2G::Init : Error - segEdge must be > " 123 <<
"1.5x the length for time delay amplitudes!!!" << endl;
124 cout <<
"TD length : " <<
cfg.
TDSize/rate <<
" sec" << endl;
125 cout <<
"segEdge : " <<
cfg.
segEdge <<
" sec" << endl;
126 cout <<
"Select segEdge > " <<
int(1.5*(
cfg.
TDSize/rate)+0.5) << endl << endl;
138 int layers = level>0 ? 1<<level : 0;
142 if(check_layers!=
nRES) {
143 cout <<
"cwb2G::Init - Error : analysis layers do not match the MRA catalog" << endl;
144 cout << endl <<
"analysis layers : " << endl;
146 int layers = level>0 ? 1<<level : 0;
147 cout <<
"level : " << level <<
" layers : " << layers << endl;
149 cout << endl <<
"MRA catalog layers : " << endl;
157 int layers = level>0 ? 1<<level : 0;
159 cout <<
"level : " << level <<
"\t rate(hz) : " << rate <<
"\t layers : " << layers
160 <<
"\t df(hz) : " <<
rateANA/2./double(1<<level)
161 <<
"\t dt(ms) : " << 1000./rate << endl;
171 double dt_max = 1./rate_min;
172 if(fmod(rate_min,1.)) {
173 cout <<
"cwb2G::Init - Error : rate min=" << rate_min <<
"(Hz) is not integer" << endl << endl;
177 cout <<
"cwb2G::Init - Error : lagStep=" <<
cfg.
lagStep <<
"(sec)" 178 <<
" is not a multple of 2*max_time_resolution=" << 2*dt_max <<
"(sec)" << endl << endl;
182 cout <<
"cwb2G::Init - Error : segEdge=" <<
cfg.
segEdge <<
"(sec)" 183 <<
" is not a multple of 2*max_time_resolution=" << 2*dt_max <<
"(sec)" << endl << endl;
187 cout <<
"cwb2G::Init - Error : segMLS=" <<
cfg.
segMLS <<
"(sec)" 188 <<
" is not a multple of 2*max_time_resolution=" << 2*dt_max <<
"(sec)" << endl << endl;
226 std::vector<wavearray<double> >
y;
232 WDM<double> WDMwhite(layers_high,layers_high,6,10);
235 double wdmFLen = double(WDMwhite.
m_H)/
rateANA;
238 cout <<
"cwb2G::ReadData : Error - filter scratch must be <= cwb scratch!!!" << endl;
239 cout <<
"filter length : " << wdmFLen <<
" sec" << endl;
240 cout <<
"cwb scratch : " <<
cfg.
segEdge <<
" sec" << endl;
243 cout <<
"WDM filter length for regression = " << wdmFLen <<
" (sec)" << endl;
253 {cout <<
"cwb2G::ReadData - Error opening root file : " <<
jname << endl;
EXIT(1);}
254 TDirectory* cdstrain = (TDirectory*)
jfile->Get(
"strain");
255 if(cdstrain==NULL) cdstrain =
jfile->mkdir(
"strain");
271 if(TMath::IsNaN(x.
mean()))
272 {cout <<
"cwb2G::ReadData - Error : found NaN in strain data !!!" << endl;
EXIT(1);}
275 {cout <<
"cwb2G::ReadData - input rate from frame " << x.
rate()
276 <<
" do not match the one defined in config : " <<
cfg.
inRate << endl;
EXIT(1);}
295 if(
i>0 && xstart != x.
start()) {
296 cout <<
"cwb2G::ReadData - Error : ifo noise data not synchronized" << endl;
297 cout <<
ifo[
i] <<
" " << x.
start() <<
" != " <<
ifo[0] <<
" " << xstart << endl;
300 if(
i>0 && xrate != x.
rate()) {
301 cout <<
"cwb2G::ReadData - Error : ifo noise data have different rates" << endl;
302 cout <<
ifo[
i] <<
" " << x.
rate() <<
" != " <<
ifo[0] <<
" " << xrate << endl;
307 TDirectory* cdmdc = (TDirectory*)
jfile->Get(
"mdc");
308 if(cdmdc==NULL) cdmdc =
jfile->mkdir(
"mdc");
316 {cout <<
"cwb2G::ReadData - input rate from frame " << x.
rate()
317 <<
" do not match the one defined in config : " <<
cfg.
inRate << endl;
EXIT(1);}
320 if(TMath::IsNaN(x.
mean()))
321 {cout <<
"cwb2G::ReadData - Error : found NaN in MDC data !!!" << endl;
EXIT(1);}
345 double layer =
j+0.01;
358 if(xstart != x.
start()) {
359 cout <<
"cwb2G::ReadData - Error : mdc/noise data with different start time" << endl;
360 printf(
"start time : noise = %10.6f - mdc = %10.6f\n",xstart,x.
start());
363 if(xrate != x.
rate()) {
364 cout <<
"cwb2G::ReadData - Error : mdc/noise data with different rate" << endl;
365 printf(
"rate : noise = %10.6f - mdc = %10.6f\n",xrate,x.
rate());
368 if(xsize != x.
size()) {
369 cout <<
"cwb2G::ReadData - Error : mdc/noise data with different buffer size" << endl;
370 printf(
"buffer size : noise = %d - mdc = %lu\n",xsize,x.
size());
379 TDirectory* cdmdc = (TDirectory*)
jfile->Get(
"mdc");
380 if(cdmdc==NULL) cdmdc =
jfile->mkdir(
"mdc");
383 std::vector<double> mdcFactor;
384 for (
int k=0;
k<(
int)
pD[0]->ISNR.size();
k++) {
387 for(
int i=0;
i<
nIFO;
i++) snr+=
pD[0]->ISNR.data[
k];
389 for(
int i=0;
i<
nIFO;
i++) snr+=
pD[
i]->ISNR.data[
k];
392 if(snr>0) mdcFactor.push_back(1./snr);
else mdcFactor.push_back(0.);
395 size_t K = mdcFactor.size();
396 for(
int k=0;
k<(
int)K;
k++) {
397 if(mdcFactor[
k]) cout << k <<
" mdcFactor : " << mdcFactor[
k] << endl;
410 for(
int k=0;
k<(
int)
pD[i]->IWFP.size();
k++) {
415 for(
int k=0;
k<(
int)K;
k++) {
423 for(
int k=0;
k<(
int)K;
k++) {
424 int ilog[5] = {1,3,12,13,14};
425 for(
int l=0;
l<5;
l++) {
426 double mfactor =
l<2 ? mdcFactor[
k] : mdcFactor[
k]*mdcFactor[
k];
480 TDirectory* cdrms=NULL;
481 TDirectory* cdcstrain=NULL;
485 WDM<double> WDMwhite(layers_high,layers_high,6,10);
490 double wdmFLen = double(WDMwhite.
m_H)/
rateANA;
493 cout <<
"cwb2G::DataConditioning : Error - filter scratch must be <= cwb scratch!!!" << endl;
494 cout <<
"filter length : " << wdmFLen <<
" sec" << endl;
495 cout <<
"cwb scratch : " <<
cfg.
segEdge <<
" sec" << endl;
498 cout <<
"WDM filter max length = " << wdmFLen <<
" (sec)" << endl;
506 {cout <<
"cwb2G::DataConditioning - Error : file " <<
jname <<
" not found" << endl;
EXIT(1);}
510 TDirectory* dcstrain = (TDirectory*)
jfile->Get(
"cstrain");
511 if(dcstrain==NULL) dcstrain=
jfile->mkdir(
"cstrain");
512 char cdcstrain_name[32];
sprintf(cdcstrain_name,
"cstrain-f%d",ifactor);
513 cdcstrain=dcstrain->mkdir(cdcstrain_name);
515 TDirectory* drms = (TDirectory*)
jfile->Get(
"rms");
516 if(drms==NULL) drms=
jfile->mkdir(
"rms");
517 char cdrms_name[32];
sprintf(cdrms_name,
"rms-f%d",ifactor);
518 cdrms=drms->mkdir(cdrms_name);
525 if(TMath::IsNaN(px->
mean()))
526 {cout <<
"cwb2G::DataConditioning - Error : found NaN in strain data !!!" << endl;
EXIT(1);}
534 if(TMath::IsNaN(px->
mean()))
535 {cout <<
"cwb2G::DataConditioning - Error : found NaN in MDC data !!!" << endl;
EXIT(1);}
538 int nshift =
int(factor*px->
rate());
540 int jstart = nshift<0 ? -nshift : 0;
541 int jstop = nshift<0 ? px->
size() : px->
size()-nshift;
542 for(
int j=jstart;
j<jstop;
j++) px->
data[
j+nshift] = pX.data[
j];
544 double tshift=nshift/px->
rate();
549 for(
int k=0;
k<(
int)
pD[i]->IWFP.size();
k++) {
554 for(
int k=0;
k<(
int)
pD[i]->TIME.size();
k++)
pD[i]->TIME[
k]+=tshift;
558 vector<TString>
ifos(nIFO);
585 rr.
add(*
hot[i],const_cast<char*>(
"target"));
607 double layer =
j+0.01;
621 sprintf(info,
"-IFO:%d-RMS:%g",i,
hot[i]->rms());
644 cout<<
"After "<<
ifo[
i]<<
" data conditioning"<<endl;
663 if(ifile!=NULL) ifile->Close();
671 vector<TString> delObjList;
691 char cdrms_name[32];
sprintf(cdrms_name,
"rms-f%d",ifactor);
692 char cdcstrain_name[32];
sprintf(cdcstrain_name,
"cstrain-f%d",ifactor);
695 TFile*
ifile =
new TFile(fName);
697 {cout <<
"cwb2G::DataConditioning - Error opening root file : " << fName << endl;
EXIT(1);}
701 {cout <<
"cwb::DataConditioning - Error opening root file : " <<
jname << endl;
EXIT(1);}
704 TDirectory* jcdrms = NULL;
705 TDirectory* jcdcstrain = NULL;
707 TDirectory* dcstrain = (TDirectory*)
jfile->Get(
"cstrain");
708 if(dcstrain==NULL) dcstrain=
jfile->mkdir(
"cstrain");
709 jcdcstrain=dcstrain->mkdir(cdcstrain_name);
711 TDirectory* drms = (TDirectory*)
jfile->Get(
"rms");
712 if(drms==NULL) drms=
jfile->mkdir(
"rms");
713 jcdrms=drms->mkdir(cdrms_name);
732 {cout <<
"cwb2G::DataConditioning - Error : cstrain not present, job terminated!!!" << endl;
EXIT(1);}
734 {
jfile->cd();jcdcstrain->cd();pws->Write(
ifo[
i]);}
744 {cout <<
"cwb2G::DataConditioning - Error : strain rms not present, job terminated!!!" << endl;
EXIT(1);}
746 {
jfile->cd();jcdcstrain->cd();jcdrms->cd();pws->Write(
ifo[
i]);}
791 char sfactor[8];
sprintf(sfactor,
"%d",ifactor);
804 std::vector<SSeries<double> >
vSS;
805 for(
int n=0;
n<
nIFO;
n++) vSS.push_back(SS);
811 {cout <<
"cwb2G::Coherence - Error : file " <<
jname <<
" not found" << endl;
EXIT(1);}
814 TDirectory* sdir = (TDirectory*)
jfile->Get(
"csparse");
815 if(sdir==NULL) sdir =
jfile->mkdir(
"csparse");
818 char sfactor[8];
sprintf(sfactor,
"%d",ifactor);
824 int upN =
rateANA/1024;
if(upN<1) upN=1;
830 int layers = level>0 ? 1<<level : 0;
832 cout <<
"level : " << level <<
"\t rate(hz) : " << rate
833 <<
"\t layers : " << layers <<
"\t df(hz) : " <<
rateANA/2./double(1<<level)
834 <<
"\t dt(ms) : " << 1000./rate << endl;
852 sprintf(cmd,
"gILEVEL = %lu;",level);
EXPORT(
size_t,gILEVEL,cmd);
864 cout<<
"thresholds in units of noise variance: Eo="<<Eo<<
" Emax="<<Eo*2<<endl;
866 sprintf(info,
"-RES:%d-THR:%g",
i,Eo);
870 cout<<
"live time in zero lag: "<<TL<<endl;
876 sprintf(cmd,
"gILEVEL = %lu;",level);
EXPORT(
size_t,gILEVEL,cmd);
885 vSS[
n].SetMap(&WS[n]);
886 vSS[
n].SetHalo(
mTau);
889 vSS[1].SetMap(&WS[1]);
890 vSS[1].SetHalo(
mTau);
896 if(
cfg.
simulation) {cout<<
"ifactor|clusters|pixels ";cout.flush();}
897 else {cout<<
"lag|clusters|pixels "; cout.flush();}
898 int csize_tot=0;
int psize_tot=0;
905 wc.
cpf(*(pwc),
false);
912 pwc->
write(
jfile,
"coherence",
"clusters",0,cycle);
914 cout<<cycle<<
"|"<<pwc->
csize()<<
"|"<<pwc->
size()<<
" ";cout.flush();
915 csize_tot+=pwc->
csize(); psize_tot+=pwc->
size();
924 TList*
fList = gDirectory->GetList();
926 TIter nextobj(fList);
927 while ((obj = (TObject *) nextobj())) {
928 if(
TString(obj->GetName()).Contains(
"clusters"))
delete obj;
931 for(
int n=0;
n<
nIFO;
n++) {vSS[
n].UpdateSparseTable();}
939 vSS[
n].Write(wsname,TObject::kWriteDelete);
942 watchCoh.Stop(); cout<<endl;
943 PrintElapsedTime(watchCoh.RealTime(),
"Coherence Elapsed Time for this level : ");
944 cout<<endl; watchCoh.Reset();
952 char sfactor[8];
sprintf(sfactor,
"%d",ifactor);
988 {cout <<
"cwb2G::SuperCluster - Error opening : " <<
jname << endl;
EXIT(1);}
991 if(
froot==NULL) {cout <<
"cwb2G::SuperCluster - Error opening : " <<
froot->GetPath() << endl;
EXIT(1);}
992 TDirectory* wdir = (TDirectory*)
froot->Get(
"histogram");
993 if(wdir==NULL) wdir =
froot->mkdir(
"histogram");
1013 cout <<
"Loading sparse TF map ... " << endl;
1024 cout <<
"cwb2G::SuperCluster : sparse map " << swname
1025 <<
" not exist in job file" << endl;
EXIT(1);
1028 pD[
n]->
vSS.push_back(SS);
1036 if(ifile!=NULL)
sprintf(cmd,
"gIFILE = (void*)%p;",ifile);
1037 else sprintf(cmd,
"gIFILE = NULL;");
1038 EXPORT(
void*,gIFILE,cmd);
1056 if(ifile!=NULL) wc.
read(ifile,
"coherence",
"clusters",0,cycle);
1057 else wc.
read(
jfile,
"coherence",
"clusters",0,cycle);
1060 for(
int i=nRES-1;
i>=0;
i--)
1063 for(
int i=nRES-1;
i>=0;
i--)
1066 cout<<
"-----------------------------------------------------"<<endl;
1067 char cycle_name[32];
1069 else sprintf(cycle_name,
" lag=%d",cycle);
1070 cout<<
"-> Processing " <<cycle_name<<
" ..."<<endl;
1071 cout<<
" --------------------------------------------------"<<endl;
1072 cout<<
" coher clusters|pixels : " 1073 <<setfill(
' ')<<setw(6)<<wc.
csize()<<
"|"<<wc.
size()<<endl;
1079 cout<<
" super clusters|pixels : " 1080 <<setfill(
' ')<<setw(6)<<wc.
esize(0)<<
"|"<<wc.
psize(0)<<endl;
1085 cout<<
" defrag clusters|pixels : " 1086 <<setfill(
' ')<<setw(6)<<wc.
esize(0)<<
"|"<<wc.
psize(0)<<
"\n";
1091 pwc->
cpf(wc,
false);
1102 double pfrac = ptot>0 ? double(psel)/double(ptot) : 0.;
1104 if(count<10000)
break;
1106 cout<<
" subnet clusters|pixels : " 1107 <<setfill(
' ')<<setw(6)<<
NET.
events()<<
"|"<<pwc->
psize(-1)<<
"\n";
1112 cout<<
" defrag clusters|pixels : " 1113 <<setfill(
' ')<<setw(6)<<
NET.
events()<<
"|"<<pwc->
psize(-1)<<
"\n";
1117 nnn += pwc->
psize(-1);
1121 pwc->
write(
jfile,
"supercluster",
"clusters",0,cycle);
1122 pwc->
write(
jfile,
"supercluster",
"clusters",-1,cycle);
1126 cout<<endl;cout.flush();
1130 cout<<endl<<
"Supercluster done"<<endl;
1131 if(mmm) cout<<
"total clusters|pixels|frac : " 1132 <<setfill(
' ')<<setw(6)<<nevt<<
"|"<<nnn<<
"|"<<nnn/double(mmm)<<
"\n";
1133 else cout<<
"total clusters : "<<nevt<<
"\n";
1134 cout<<endl;cout.flush();
1137 if(mmm)
sprintf(info,
"-NEVT:%d-PSIZE:%d-FRAC:%g",nevt,nnn,nnn/
double(mmm));
1138 else sprintf(info,
"-NEVT:%d",nevt);
1175 vector<TString> delObjList;
1178 delObjList.push_back(
"coherence");
1182 char cdrms_name[32];
sprintf(cdrms_name,
"rms/rms-f%d",ifactor);
1183 char cdcstrain_name[32];
sprintf(cdcstrain_name,
"cstrain/cstrain-f%d",ifactor);
1184 delObjList.push_back(cdrms_name);
1185 delObjList.push_back(cdcstrain_name);
1227 EXPORT(
float*,gSLAGSHIFT,TString::Format(
"gSLAGSHIFT = (float*)%p;",
slagShift).Data());
1232 jfile = (xname==
jname) ?
new TFile(xname,
"UPDATE") :
new TFile(xname);
1234 {cout <<
"cwb2G::Likelihood - Error : file " << xname.Data() <<
" not found" << endl;
EXIT(1);}
1244 ifile =
new TFile(
jname,
"UPDATE");
1245 if(ifile==NULL||!ifile->IsOpen()) {
1246 cout <<
"cwb2G::Likelihood - Error : file " <<
jname <<
" not found" << endl;
EXIT(1); }
1253 if(xname!=
jname) ifile->Close();
1261 cout <<
"Loading sparse TF map ... " << endl;
1270 cout <<
"cwb2G::Likelihood : sparse map " << swname
1271 <<
" not exist in job file" << endl;
EXIT(1);
1274 pD[
n]->
vSS.push_back(SS);
1281 {netburst->
dopen(outDump,const_cast<char*>(
"a"),
true); netburst->
dclose();}
1288 EXPORT(
int,gLRETRY,
"gLRETRY = 0;")
1290 EXPORT(
int,gRET,
"gRET = 0;")
1294 int gRET=0;
IMPORT(
int,gRET)
if(gRET==-1)
continue;
1296 int gLRETRY=0;
IMPORT(
int,gLRETRY) LRETRY=gLRETRY;
1304 vector<int> clist = pwc->
read(
jfile,
"supercluster",
"clusters",0,cycle);
1308 char cycle_name[32];
1310 else sprintf(cycle_name,
" lag=%d",cycle);
1311 cout<<
"-------------------------------------------------------"<<endl;
1312 cout<<
"-> Processing " << clist.size() <<
" clusters in" << cycle_name<<endl;
1313 cout<<
" ----------------------------------------------------"<<endl;cout.flush();
1319 int nselected_core_pixels = 0;
1320 int nrejected_weak_pixels = 0;
1321 int nrejected_loud_pixels = 0;
1322 for(
int k=0;
k<(
int)clist.size();
k++) {
1325 pwc->
read(
jfile,
"supercluster",
"clusters",nmax,cycle,0,clist[
k]);
1331 int id = size_t(cid.
data[cid.
size()-1]+0.1);
1339 int selected_core_pixels = 0;
1350 if(
cfg.
dump) netburst->
dopen(outDump,const_cast<char*>(
"a"),
false);
1354 int rejected_weak_pixels = 0;
1355 int rejected_loud_pixels = 0;
1360 cout<<
" cluster-id|pixels: "<<setfill(
' ')<<setw(5)<<clist[
k]<<
"|"<<pwc->
size()-npixels;
1361 if(detected) cout <<
"\t -> SELECTED !!!" << endl;
1362 else cout <<
"\t <- rejected " << endl;
1369 ifile =
new TFile(
jname,
"UPDATE");
1370 if(ifile==NULL||!ifile->IsOpen()) {
1371 cout <<
"cwb2G::Likelihood - Error : file " <<
jname <<
" not found" << endl;
EXIT(1); }
1373 pwc->
write(ifile,
"likelihood",
"clusters",0,cycle);
1374 pwc->
write(ifile,
"likelihood",
"clusters",-1,cycle,0,k+1);
1375 if(detected) cout<<
"saved"<<endl;cout.flush();
1377 if(xname!=
jname) ifile->Close();
1380 if(detected) nevents++;
1381 if(gLRETRY==0) npixels=pwc->
size();
1382 if(!detected && gLRETRY==LRETRY) {
1384 EXPORT(
int,gLRETRY,TString::Format(
"gLRETRY = %d;",gLRETRY).Data())
1394 ifile =
new TFile(
jname,
"UPDATE");
1395 if(ifile==NULL||!ifile->IsOpen()) {
1396 cout <<
"cwb2G::Likelihood - Error : file " <<
jname <<
" not found" << endl;
EXIT(1); }
1401 cout<<
"dump ced to job file ..." <<endl;
1402 TDirectory* cdced = NULL;
1403 cdced = (TDirectory*)ifile->Get(
"ced");
1404 if(cdced == NULL) cdced = ifile->mkdir(
"ced");
1407 cout<<
"dump ced to disk ..." <<endl;
1429 if(
ced->
Write(ofactor,fullCED)) ceddir = 1;
1431 if(xname!=
jname) ifile->Close();
1439 clist = pwc->
read(
jfile,
"supercluster",
"clusters",0,cycle);
1440 if(pwc->
cList.size()) {
1442 for(
int p=0;p<
npix;p++) pwc->
pList.pop_back();
1448 }
else gLRETRY=LRETRY;
1449 EXPORT(
int,gLRETRY,TString::Format(
"gLRETRY = %d;",gLRETRY).Data())
1463 cout<<endl<<endl;cout.flush();
1468 sprintf(info,
"-NEVT:%d",NEVENTS);
1475 vector<TString> delObjList;
1499 TDirectory* sdir = (TDirectory*)jfile->Get(tdir);
1500 if(sdir==NULL) sdir = jfile->mkdir(tdir);
1507 std::vector<SSeries<double> >
vSS;
1508 for(
int n=0;
n<
nIFO;
n++) vSS.push_back(ss);
1516 vSS[
n].SetMap(
pTF[n]);
1517 vSS[
n].SetHalo(
mTau);
1521 vSS[1].SetMap(
pTF[1]);
1522 vSS[1].SetHalo(
mTau);
1532 wc.
read(jfile,tname,
"clusters",0,cycle);
1535 if(tname==
"coherence") {
1536 wc.
read(jfile,tname,
"clusters",-2,cycle,-vSS[0].wrate());
1538 wc.
read(jfile,tname,
"clusters",-1,cycle,vSS[0].wrate());
1548 for(
int n=0;
n<
nIFO;
n++) {vSS[
n].UpdateSparseTable();}
1554 for(
int n=0;
n<
nIFO;
n++) {ncore+=vSS[
n].GetSparseSize();ncluster+=vSS[
n].GetSparseSize(0);}
1555 ccluster = 2*(vSS[0].GetHaloSlice()+vSS[0].GetHaloSlice(
true))+1;
1556 ccluster*= 2*vSS[0].GetHaloLayer()+1;
1577 vSS[
n].Write(wsname,TObject::kWriteDelete);
1582 cout <<
"----------------- FINAL SPARSE STATISTIC / DETECTOR -------------------" << endl;
1583 cout <<
"npix_core_tot|npix_cluster_tot|3*npix_cluster_tot|ccluster_tot|ratio : " << endl;
1584 cout << ncore_tot <<
" | " << ncluster_tot <<
" | " << 3*ncluster_tot
1585 <<
" | " << ccluster_tot <<
" | " << (ccluster_tot?double(3*ncluster_tot)/(ccluster_tot):0) << endl;
1605 cout <<
"cwb2G::FillSparseTFmap ..." << endl;
1612 pD[
n]->
vSS.push_back(SS);
1618 pD[1]->
vSS.push_back(SS);
1630 wc.
read(jfile,
"coherence",
"clusters",0,cycle);
1633 if(tname==
"coherence") {
1634 wc.
read(jfile,tname,
"clusters",-2,cycle,-
pD[0]->
vSS[
i].wrate());
1636 wc.
read(jfile,tname,
"clusters",-1,cycle,
pD[0]->
vSS[
i].wrate());
1668 {cout <<
"cwb2G::SaveWaveforms - NULL input root file : " <<
jname << endl;
EXIT(1);}
1670 {cout <<
"cwb2G::SaveWaveforms - save mode must be 1/2/3 " << endl;
EXIT(1);}
1673 TDirectory*
wf = (TDirectory*)jfile->Get(
"waveform");
1674 if(wf==NULL) wf=jfile->mkdir(
"waveform");
1675 char cdwf_name[32];
sprintf(cdwf_name,
"waveform-f%d",ifactor);
1676 TDirectory* cdwf = (TDirectory*)wf->Get(cdwf_name);
1677 if(cdwf==NULL) cdwf=wf->mkdir(cdwf_name);
1681 pD->Write(pD->
Name,TObject::kSingleKey);
1704 {cout <<
"cwb2G::LoadWaveforms - NULL input root file : " <<
iname << endl;
EXIT(1);}
1706 {cout <<
"cwb2G::LoadWaveforms - save mode must be 1/2/3 " << endl;
EXIT(1);}
1710 char cdwf_name[32];
sprintf(cdwf_name,
"waveform/waveform-f%d",ifactor);
1712 if(pd==NULL)
return;
char channelNamesMDC[NIFO_MAX][50]
monster wdmMRA
list of pixel pointers for MRA
void SetOptions(int simulation, double rho, double inRate, bool useSparse=false, char *gtype=const_cast< char *>("png"), int paletteId=0)
detector * getifo(size_t n)
param: detector index
virtual size_t defragment(double T, double F, TH2F *=NULL)
T - maximum time gap in seconds F - maximum frequency gap in Hz.
virtual void resize(unsigned int)
CWB::frame fr[2 *NIFO_MAX]
void PrintElapsedTime(int job_elapsed_time, TString info)
void LoadWaveforms(TFile *ifile, detector *pD, int ifactor, int wfSAVE=1)
size_t write(const char *file, int app=0)
#define REGRESSION_FILTER_LENGTH
void white(double dT=0, int wtype=1, double offset=0., double stride=0.)
what it does: see algorithm description in wseries.hh
static size_t GetProcInfo(bool mvirtual=true)
std::vector< vector_int > cRate
void WriteSparseTFmap(TFile *jfile, int ifactor, TString tdir, TString tname)
void setMatrix(double edge=0., double f=1.)
void SaveWaveforms(TFile *jfile, detector *pD, int ifactor, int wfSAVE=1)
void FileGarbageCollector(TString ifName, TString ofName="", vector< TString > delObjList=vector< TString >())
std::vector< netcluster > wc_List
size_t add(detector *)
param: detector structure return number of detectors in the network
bool singleDetector
used for the stage stuff
WSeries< double > * pTF[NIFO_MAX]
pointers to detectors
void setAntenna(detector *)
param: detector (use theta, phi index array)
std::vector< double > * getmdcTime()
size_t setsnr(wavearray< double > &, std::vector< double > *, std::vector< double > *, double=5., double=8.)
WDM< double > * pwdm[NRES_MAX]
virtual void rate(double r)
#define REGRESSION_MATRIX_FRACTION
bool Likelihood(int ifactor, char *ced_dir, netevent *output=NULL, TTree *net_tree=NULL, char *outDump=NULL)
long getNetworkPixels(int LAG, double Eo, double DD=1., TH1F *hist=NULL)
TString iname
stage benchmark
void CWB_Plugin(TFile *jfile, CWB::config *, network *, WSeries< double > *, TString, int)
COHERENCE.
detector * pD[NIFO_MAX]
noise variability
size_t setIndexMode(size_t=0)
void bandpass(wavearray< DataType_t > &ts, double flow, double fhigh, int n=-1)
TFile * jfile
output root file
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
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
void add(const wavearray< DataType_t > &, int=0, int=0, int=0)
std::vector< vector_float > sArea
std::vector< SSeries< double > > vSS[NIFO_MAX]
virtual double ReadData(double mdcShift, int ifactor)
void Coherence(int ifactor)
size_t setcore(bool core, int id=0)
virtual size_t supercluster(char atype, double S, bool core)
param: statistic: E - excess power, L - likelihood param: selection threshold T for likelihood cluste...
std::vector< vector_int > cList
size_t add(WSeries< double > &target, char *name, double fL=0., double fH=0.)
virtual void start(double s)
std::vector< double > mdcTime
void setTDFilter(int nCoeffs, int L=1)
long subNetCut(int lag, float subnet=0.6, float subcut=0.33, TH2F *hist=NULL)
network NET
pointers to WSeries
size_t cluster(int kt, int kf)
param: time gap in pixels return: number of reconstructed clusters
long likelihood2G(char mode, int lag, int ID, TH2F *hist=NULL)
std::vector< vector_float > p_Map
void dopen(const char *fname, char *mode, bool header=true)
cout<< "Selected Pixels : "<< nPix<< endl;wc.cluster(1, 1);SSeries< double > ss
fprintf(stdout,"start=%f duration=%f rate=%f\, x.start(), x.size()/x.rate(), x.rate())
virtual size_t size() const
long likelihoodWP(char mode, int lag, int ID, TH2F *hist=NULL)
void output2G(TTree *, network *, size_t, int, double)
double setVeto(double=5.)
param: time window around injections
void PrintAnalysisInfo(CWB_STAGE stage, TString comment, TString info, bool out=true, bool log=true)
#define IMPORT(TYPE, VAR)
double maxEnergy(wavearray< DataType_t > &ts, Wavelet &w, double=0, int=1, int=0, TH1F *=NULL)
void SuperCluster(int ifactor)
std::vector< netpixel > pList
wavearray< double > * getHoT()
param: no parameters
gwavearray< double > * gx
void setDelayIndex(double rate)
param: MDC log file
std::vector< std::string > mdcList
void Resample(const wavearray< DataType_t > &, double, int=6)
#define MIN_SKYRES_HEALPIX
void apply(double threshold=0., char c='a')
void FillSparseTFmap(TFile *jfile, int ifactor, TString tname)
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
printf("total live time: non-zero lags = %10.1f \, liveTot)
std::vector< vector_int > p_Ind
double dataShift[NIFO_MAX]
std::vector< float > cFreq
int Write(double factor, bool fullCED=true)
char channelNamesRaw[NIFO_MAX][50]
void SetChannelName(char *chName)
void PrintStageInfo(CWB_STAGE stage, TString comment, bool out=true, bool log=true, TString fname="")
void solve(double th, int nE=0, char c='s')
char jname[1024]
job file object
#define REGRESSION_SOLVE_EIGEN_THR
std::vector< wavearray< double > * > IWFP
size_t cpf(const netcluster &, bool=false, int=0)
double THRESHOLD(double bpp)
param: selected fraction of LTF pixels assuming Gaussian noise
void setDelay(const char *="L1")
size_t loadTDampSSE(network &net, char c, size_t BATCH=10000, size_t LOUD=0)
wavearray< double > getClean()
#define REGRESSION_APPLY_THR
#define EXPORT(TYPE, VAR, CMD)
std::vector< float > cTime
virtual double mean() const
WSeries< double > * getTFmap()
param: no parameters
unsigned int jobfOptions
history object
CWB::ced * ced
temporary time series
int SetSkyMask(network *net, CWB::config *cfg, char *options, char skycoord, double skyres=-1)
void readFrames(char *filename, char *channel, wavearray< double > &w)
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
netcluster * getwc(size_t n)
param: delay index
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
std::vector< clusterdata > cData
netevent * netburst
livetime object
virtual wavearray< double > select(char *, double)
void setMRAcatalog(char *fn)
double factors[FACTORS_MAX]
void DataConditioning(int ifactor)
wavearray< double > * hot[NIFO_MAX]
wavelet pointers: pwdm[0] - l_high, wdm[nRES-1] l_low
cout<< "live time after cat 2 : "<< detSegs_ctime<< endl;if(detSegs_ctime< segTHR) {cout<< "job segment live time after cat2 < "<< segTHR<< " sec, job terminated !!!"<< endl;exit(1);} double Tb=detSegs[0].start;double Te=detSegs[0].stop;double dT=Te-Tb;char file[512], tdf00[512], tdf90[512], buFFer[1024];int rnID=int(gRandom->Rndm(13) *1.e9);if(simulation) { i=NET.readMDClog(injectionList, double(long(Tb)) -mdcShift);printf("GPS: %16.6f saved, injections: %d\", double(long(Tb)), i);frTB[nIFO].shiftBurstMDCLog(NET.mdcList, ifos, mdcShift);for(int i=0;i< NET.mdcTime.size();i++) NET.mdcTime[i]+=mdcShift;vector< waveSegment > mdcSegs(NET.mdcTime.size());for(int k=0;k< NET.mdcTime.size();k++) {mdcSegs[k].start=NET.mdcTime[k]-gap;mdcSegs[k].stop=NET.mdcTime[k]+gap;} vector< waveSegment > mdcSegs_dq2=slagTB.mergeSegLists(detSegs_dq2, mdcSegs);double mdcSegs_ctime=slagTB.getTimeSegList(mdcSegs_dq2);cout<< "live time in zero lag after cat2+inj : "<< mdcSegs_ctime<< endl;if(mdcSegs_ctime==0) {cout<< "job segment with zero cat2+inj live time in zero lag, job terminated !!!"<< endl;exit(1);} } if(dump_infos_and_exit) exit(0);if(mask >0.) NET.setSkyMask(mask, skyMaskFile);for(i=0;i< nIFO;i++) { frTB[i].readFrames(FRF[i], channelNamesRaw[i], x);x.start(x.start()+dataShift[i]);x.start(x.start() -segLen *(segID[i]-segID[0]));if(singleDetector) TB.resampleToPowerOfTwo(x);sprintf(file,"%s/%s_%d_%s_%d_%d.dat", nodedir, ifo[i], int(Tb), data_label, runID, rnID);if(dump_sensitivity_and_exit) { sprintf(file,"%s/sensitivity_%s_%d_%s_job%d.txt", dump_dir, ifo[i], int(Tb), data_label, runID);cout<< endl<< "Dump Sensitivity : "<< file<< endl<< endl;TB.makeSpectrum(file, x);continue;} if(dcCal[i]>0.) x *=dcCal[i];if(fResample >0) { x.FFT(1);x.resize(fResample/x.rate() *x.size());x.FFT(-1);x.rate(fResample);} pTF[i]=pD[i]-> getTFmap()
std::vector< SSeries< double > > vSS
size_t read(const char *)
double ReadData(double mdcShift, int ifactor)
#define REGRESSION_SOLVE_EIGEN_NUM
#define REGRESSION_SOLVE_REGULATOR
virtual void resize(unsigned int)
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
void putLayer(wavearray< DataType_t > &, double n)
param: n - layer number
size_t setsim(WSeries< double > &, std::vector< double > *, double=5., double=8., bool saveWF=false)
void setSkyMaps(double, double=0., double=180., double=0., double=360.)
param: sky map granularity step, degrees param: theta begin, degrees param: theta end...
virtual WSeries< double > white(double, int, double=0., double=0.)
what it does: each wavelet layer is devided into k intervals.