5 TB.
checkFile(gSystem->Getenv(
"CWB_ROOTLOGON_FILE"));
6 TB.
checkFile(gSystem->Getenv(
"CWB_PARAMETERS_FILE"));
7 TB.
checkFile(gSystem->Getenv(
"CWB_UPARAMETERS_FILE"));
13 if(
search==
'b' ||
search==
'B') cout<<
"\n un-modeled search (single stream): "<<
SEARCH()<<endl;
14 if(
search==
'r' ||
search==
'R') cout<<
"\n un-modeled search (dual stream): "<<
SEARCH()<<endl;
23 cout <<
"------> cWB in Sigle Detector Mode !!!" << endl;
47 for(i=0; i<
nIFO; i++) cout<<
ifo[i]<<
" ";
48 cout<<
" detectors\n\n";
70 for(i=0; i<
nIFO; i++) NET.
add(pD[i]);
85 cout<<
"maximum time delay between detectors: "<<mTau<<endl;
86 cout<<
" maximum time delay difference: "<<dTau<<endl;
87 cout<<
" skymap search mode: "<<
mode<<endl;
88 cout<<
" skymap angular resolution: "<<
angle<<endl;
89 cout<<
" skymap size in polar angle: "<<
Theta1<<
", "<<
Theta2<<endl;
90 cout<<
" skymap size in azimuthal angle: "<<
Phi1<<
", "<<
Phi2<<endl;
91 cout<<
" netRHO and netCC: "<<
netRHO<<
", "<<
netCC<<endl;
92 cout<<
" regulator delta, local: "<<
delta<<
" "<<NET.
local<<endl<<endl;
100 gSystem->Exec(
"/bin/date");
101 gSystem->Exec(
"/bin/hostname");
106 int dump_infos_and_exit =
TString(gSystem->Getenv(
"CWB_DUMP_INFOS_AND_EXIT")).Atoi();
107 int dump_sensitivity_and_exit =
TString(gSystem->Getenv(
"CWB_DUMP_SENSITIVITY_AND_EXIT")).Atoi();
110 runID = srunID.Atoi();
112 cout<<
"job ID: "<<runID<<endl;
124 else sprintf(job_stage,
"PRODUCTION");
127 if(cwbBuffer!=NULL) {
137 for(
int i=0;i<gApplication->Argc();i++)
sprintf(cmd_line,
"%s %s",cmd_line,gApplication->Argv(i));
141 if(rootlogonBuffer!=NULL) {
148 for(
int i=0;i<gApplication->Argc()-1;i++) {
149 if(
TString(gApplication->Argv(i)).Contains(
".C")) {
150 char* parametersBuffer = histTB.
readFile(gApplication->Argv(i));
151 if(parametersBuffer!=NULL) {
152 if(
TString(gApplication->Argv(i))==cwb_parameters_name) {
158 delete [] parametersBuffer;
179 vector<TString>
ifos(nIFO);
183 cout <<
"Error : slagSize<=1 & lagSize==1 in simulation mode !!!" << endl;
exit(1);
188 cout << endl <<
"START SLAG Init ..." << endl << endl;
203 if(SLAG.
jobId!=runID) {cout <<
"jobID " << runID <<
" not found in the slag list !!!" << endl;
exit(1);}
206 for(n=0; n<
nIFO; n++) segID[n]=SLAG.
segId[n];
207 cout <<
"SuperLag=" << slagID <<
" jobID=" << jobID;
208 for(n=0; n<nIFO; n++) cout <<
" segID[" <<
ifo[
n] <<
"]=" << segID[
n];cout << endl;
219 cout <<
"Slag Segments" << endl;
221 cout << endl <<
"END SLAG Init ..." << endl << endl;
227 for(
int n=0;n<
nIFO;n++) segID[n]=0;
234 cout <<
"Standard Segments" << endl;
239 for(n=0; n<
nIFO; n++) slagShift[n]=(segID[n]-segID[0])*
segLen;
241 netburst.setSLags(slagShift);
245 for(
int i=0;i<
nIFO;i++) {
246 cout <<
"detSegs_dq1[" <<
ifo[
i] <<
"] Range : " << detSegs[
i].start <<
"-" << detSegs[
i].stop << endl;
248 if(detSegs.size()==0) {cout <<
"no segments found for this job, job terminated !!!" << endl;
exit(1);}
252 for(n=0; n<
nIFO; n++) {
253 nfrFiles[
n]=frTB[
n].frl2FrTree(
frFiles[n]);
254 cout <<
ifo[
n] <<
" -> nfrFiles : " << nfrFiles[
n] << endl;
260 cout <<
"MDC " <<
" -> nfrFiles : " << nfrFiles[
nIFO] << endl;
264 cout <<
"mdc_range : " << mdc_range.
start <<
" " << mdc_range.
stop << endl;
269 cout <<
"mdcShift : " << mdcShift << endl;
270 FRF[
nIFO] = frTB[
nIFO].getFrList(detSegs[0].start-mdcShift, detSegs[0].
stop-mdcShift,
segEdge);
280 detSegs_dq2.push_back(detSegs[0]);
282 for(
int i=0;i<detSegs_dq2.size();i++) {
283 cout <<
"detSegs_dq2[" << i <<
"] Range : " << detSegs_dq2[
i].start <<
"-" << detSegs_dq2[
i].stop << endl;
286 cout <<
"live time after cat 2 : " << detSegs_ctime << endl;
287 if(detSegs_ctime<
segTHR) {cout <<
"job segment live time after cat2 < " <<
segTHR <<
" sec, job terminated !!!" << endl;
exit(1);}
289 double Tb=detSegs[0].start;
290 double Te=detSegs[0].stop;
296 int rnID =
int(gRandom->Rndm(13)*1.e9);
300 printf(
"GPS: %16.6f saved, injections: %d\n",
double(
long(Tb)),i);
305 vector<waveSegment> mdcSegs(NET.
mdcTime.size());
307 vector<waveSegment> mdcSegs_dq2 = slagTB.
mergeSegLists(detSegs_dq2,mdcSegs);
309 cout <<
"live time in zero lag after cat2+inj : " << mdcSegs_ctime << endl;
310 if(mdcSegs_ctime==0) {cout <<
"job segment with zero cat2+inj live time in zero lag, job terminated !!!" << endl;
exit(1);}
313 if(dump_infos_and_exit)
exit(0);
317 for(i=0; i<
nIFO; i++) {
322 sprintf(file,
"%s/%s_%d_%s_%d_%d.dat",
324 if(dump_sensitivity_and_exit) {
326 cout << endl <<
"Dump Sensitivity : " << file << endl << endl;
341 fprintf(stdout,
"start=%f duration=%f rate=%f\n",
344 cout <<
"net.C - Error : data not synchronized" << endl;
345 cout <<
ifo[
i] <<
" " << x.
start() <<
" != " <<
ifo[0] <<
" " << pTF[0]->
start() << endl;
349 cout <<
"net.C - Error : data have different rates" << endl;
350 cout <<
ifo[
i] <<
" " << x.
rate() <<
" != " <<
ifo[0] <<
" " << pTF[0]->
rate() << endl;
357 sprintf(file,
"%s/mdc%s_%d_%s_%d_%d.dat",
369 fprintf(stdout,
"start=%f duration=%f rate=%f\n",
375 if(singleDetector)
break;
377 if(dump_sensitivity_and_exit)
exit(0);
400 for(
int iii=0; iii<
nfactor; iii++) {
404 cout<<
"factor="<<factor<<endl;
416 if(!gSystem->GetPathInfo(endFile,fstemp)) {
417 printf(
"The file %s already exists - skip\n",endFile);
420 if(!rf.IsZombie())
continue;
425 char prod_label[512];
426 sprintf(prod_label,
"%d_%d_%s_slag%d_lag%d_%d_job%d",
438 cout<<
"output file on the node: "<<outFile<<endl;
439 cout<<
"final output file name : "<<endFile<<endl;
440 cout<<
"temporary output file : "<<tmpFile<<endl;
448 for(i=0; i<
nIFO; i++) {
449 sprintf(file,
"%s/%s_%d_%s_%d_%d.dat",
455 sprintf(file,
"%s/mdc%s_%d_%s_%d_%d.dat",
465 pD[
i]->
white(60.,1,8.,20.);
474 cout<<
"After "<<
ifo[
i]<<
" data conditioning"<<endl;
477 if(singleDetector) {*pD[1]=*pD[0];
break;}
495 TFile *
froot =
new TFile(tmpFile,
"RECREATE");
496 TTree*
net_tree = netburst.setTree();
497 TTree* live_tree= live.
setTree();
500 TTree* mdc_tree =
mdc.setTree();
503 TTree* var_tree = wavevar.
setTree();
504 TTree* noise_tree = noiserms.
setTree();
513 Long_t xid,xsize,xflags,xmt;
514 int xestat = gSystem->GetPathInfo(outDump,&xid,&xsize,&xflags,&xmt);
516 sprintf(command,
"/bin/rm %s",outDump);
517 gSystem->Exec(command);
525 for(
int mlag=mlagOff;mlag<mlagSize;mlag+=
mlagStep) {
529 if(lagSize==0)
continue;
531 cout <<
"lagSize : " << lagSize <<
" lagOff : " << lagOff << endl;
536 cout<<
"lag step: "<<
lagStep<<endl;
537 cout<<
"number of time lags: "<<lags<<endl;
542 cout<<
"live time in zero lag: "<<TL<<endl<<endl;
543 if(TL <= 0.)
exit(1);
551 cout<<
"Start coherent search: "; gSystem->Exec(
"/bin/date");
566 cout<<
"pixel threshold in units of noise rms: "<<Ao<<endl;
567 cout<<
"2 OR threshold in units of noise var: "<<
x2or*Ao*Ao<<endl;
569 cout<<
"total pixels: "<<NET.
coherence(Ao)<<
" ";
571 n = size_t(2.*
Tgap*pD[0]->
getTFmap()->resolution(0)+0.1);
572 m = size_t(
Fgap/pD[0]->
getTFmap()->resolution(0)+0.0001);
574 cout<<
"clusters: "<<NET.
cluster(n,m)<<
" ";
577 for(j=0; j<NET.
nLag; j++) {
578 sprintf(file,
"%s/pix_%d_%s_%d_%d_%d.lag",
580 wc = *(NET.
getwc(j));
581 if(!append)
np.data[
j] = wc.
write(file,0);
582 else np.data[
j] += wc.
write(file,1);
598 for(j=0; j<lags; j++){
599 sprintf(file,
"%s/pix_%d_%s_%d_%d_%d.lag",
603 pwc = NET.
getwc(j); pwc->
cpf(wc,
true);
604 cout<<m<<
"|"<<pwc->
size()<<
" ";
605 sprintf(command,
"/bin/rm %s",file);
606 gSystem->Exec(command);
609 cout<<
"events in the buffer: "<<NET.
events()<<
"\n";
628 cout<<
"rejected weak pixels: "<<NET.
netcut(
netRHO,
'r',0,1)<<
"\n";
629 cout<<
"rejected loud pixels: "<<NET.
netcut(
netCC,
'c',0,1)<<
"\n";
630 cout<<
"events in the buffer: "<<NET.
events()<<
"\n";
633 cout<<
"dump ced into "<<out_CED<<
"\n";
636 CWB::ced ced(&NET,&netburst,out_CED);
639 bool fullCED = singleDetector ?
false :
true;
640 if(ced.Write(factor,fullCED)) ceddir = 1;
646 gSystem->Exec(
"/bin/date");
648 cout<<
"\nSearch done\n";
649 cout<<
"reconstructed events: "<<NET.
events()<<
"\n";
657 if(
dump) netburst.dopen(outDump,
"w");
660 live.
output(live_tree,&NET,slagShift);
663 netburst.output(net_tree,&NET,factor);
664 mdc.output(mdc_tree,&NET,factor);
667 netburst.output(net_tree,&NET);
668 for(i=0; i<
nIFO; i++) {
670 noiserms.
output(noise_tree,&pD[i].nRMS,i+1,R/2);
678 if(
dump) netburst.dclose();
684 sprintf(command,
"/bin/mv %s %s", tmpFile, outFile);
685 if(!
cedDump) gSystem->Exec(command);
686 sprintf(command,
"/bin/mv %s %s",outFile,endFile);
687 if(!
cedDump) gSystem->Exec(command);
688 sprintf(command,
"/bin/mv %s %s",outDump,endDump);
690 xestat = gSystem->GetPathInfo(end_CED,&xid,&xsize,&xflags,&xmt);
692 sprintf(command,
"/bin/mv %s/* %s/.",out_CED,end_CED);
694 sprintf(command,
"/bin/mv %s %s",out_CED,end_CED);
696 if(
cedDump && ceddir) gSystem->Exec(command);
702 for(i=0; i<
nIFO; i++) {
704 sprintf(command,
"/bin/rm %s",file);
705 gSystem->Exec(command);
709 sprintf(command,
"/bin/rm %s",file);
710 gSystem->Exec(command);
712 if(singleDetector)
break;
714 cout<<
"Stopping the job "<<runID<<endl;
715 gSystem->Exec(
"/bin/date");
725 printf(
"Job Elapsed Time - %02d:%02d:%02d (hh:mm:ss)\n",job_elapsed_hour,job_elapsed_min,job_elapsed_sec);
726 printf(
"Job Speed Factor - %2.1fX\n",job_speed_factor);
virtual size_t size() const
size_t write(const char *file, int app=0)
void white(double dT=0, int wtype=1, double offset=0., double stride=0.)
what it does: see algorithm description in wseries.hh
size_t add(detector *)
param: detector structure return number of detectors in the network
printf("total live time: non-zero lags = %10.1f \n", liveTot)
void set2or(double p)
param: threshold
size_t readMDClog(char *, double=0., int=11, int=12)
param: MDC log file param: approximate gps time
void setAntenna(detector *)
param: detector (use theta, phi index array)
std::vector< double > * getmdcTime()
fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(), x.size()/x.rate(), x.rate())
char channelNamesMDC[NIFO_MAX][50]
virtual void rate(double r)
WSeries< double > * pTF[nIFO]
char channelNamesRaw[NIFO_MAX][50]
virtual WSeries< float > variability(double=0., double=-1., double=-1.)
param: first - time window to calculate normalization constants second - low frequency boundary for c...
char * watversion(char c='s')
size_t setIndexMode(size_t=0)
void add(const wavearray< DataType_t > &, int=0, int=0, int=0)
virtual void ReadBinary(const char *, int=0)
char frFiles[NIFO_MAX+1][256]
size_t setSkyMask(double f, char *fname)
CWB::Toolbox frTB[nIFO+1]
virtual size_t supercluster(char atype, double S, bool core)
param: statistic: E - excess power, L - likelihood param: selection threshold T for likelihood cluste...
virtual void start(double s)
vector< waveSegment > detSegs_dq2
std::vector< double > mdcTime
long coherence(double, double=0., double=0.)
param: threshold on lognormal pixel energy (in units of noise rms) param: threshold on total pixel en...
size_t cluster(int kt, int kf)
param: time gap in pixels return: number of reconstructed clusters
void setRunID(size_t n)
param: run
void constraint(double d=1., double g=0.0001)
param: constraint parameter, p=0 - no constraint
void Forward(size_t k)
param: number of steps
double setVeto(double=5.)
param: time window around injections
void setDelayIndex(double rate)
param: MDC log file
std::vector< std::string > mdcList
double dataShift[NIFO_MAX]
double getDelay(const char *c="")
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
std::vector< waveSegment > segList
size_t cpf(const netcluster &, bool=false, int=0)
void output(TTree *waveTree, network *net, float *slag=NULL, vector< waveSegment > detSegs=DEFAULT_WAVESEGMENT)
void setDelay(const char *="L1")
void bandPass(double f1, double f2, double a=0.)
long likelihood(char='E', double=sqrt(2.), int=0, size_t=0, int=-1, bool=false)
void output(TTree *, WSeries< double > *, int=0, double=1.)
virtual void DumpBinary(const char *, int=0)
vector< waveSegment > detSegs
cout<<" network of ";for(i=0;i< nIFO;i++) cout<< ifo[i]<<" ";cout<<" detectors\n\n";Meyer< double > B(1024)
WSeries< double > * getTFmap()
param: no parameters
TString cwb_parameters_name
virtual void lprFilter(double, int=0, double=0., double=0.)
void AddLog(char *Stage, char *Log, TDatime *Time=NULL)
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
vector< TString > ifos(nIFO)
strcpy(RunLabel, RUN_LABEL)
Meyer< double > S(1024, 2)
netcluster * getwc(size_t n)
param: delay index
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
vector< waveSegment > cat2List
void setDelayFilters(detector *=NULL)
size_t netcut(double, char='L', size_t=0, int=1)
param: threshold param: minimum cluster size processed by the corrcut param: cluster type return: num...
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\n", 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()
size_t read(const char *)
virtual void resize(unsigned int)
int setTimeShifts(size_t=1, double=1., size_t=0, size_t=0, const char *=NULL, const char *="w", size_t *=NULL)
param number of time lags param time shift step in seconds param first lag ID param maximum lag ID pa...
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
size_t setsim(WSeries< double > &, std::vector< double > *, double=5., double=8., bool saveWF=false)
void AddHistory(char *Stage, char *Type, char *History, TDatime *Time=NULL)
void setSkyMaps(double, double=0., double=180., double=0., double=360.)
param: sky map granularity step, degrees param: theta begin, degrees param: theta end...
double threshold(double, double)
param: selected fraction of LTF pixels assuming Gaussian noise param: maximum time delay between dete...
vector< waveSegment > cat1List