20 #include "TTreeFormula.h" 23 #include <Riostream.h> 25 #define CCAST(PAR) const_cast<char*>(PAR) \ 34 #define MAX_THREADS 16 36 #define MAX_TREE_SIZE 100000000000LL 51 MergeParms mp = *((MergeParms*) ptr);
55 sprintf(ofName,
"%s/mergeList.T%d.txt",mp.odir.Data(),mp.threadID);
61 sprintf(cmd,
"root -l -b '%s/cwb_merge_thread.C(\"%s\",%d,\"%s\",\"%s\",%d,%d,%d)'",
62 gSystem->ExpandPathName(
"$CWB_MACROS"),
ofName, mp.simulation,
63 mp.odir.Data(), mp.label.Data(), mp.brms, mp.bvar, mp.bpsm);
67 gSystem->Exec(TString::Format(
"rm %s",ofName));
87 cout <<
"CWB::Toolbox::readSegments - Error Opening File : " << ifile << endl;
97 vector<waveSegment>
iseg;
100 in.getline(str,1024);
101 if(str[0] ==
'#')
continue;
102 in.seekg(fpos, ios::beg);
105 if (!in.good())
break;
107 in.seekg(fpos+1, ios::beg);
135 vector<waveSegment> vsegs;
136 int n = ilist.size();
137 if(n==0)
return vsegs;
143 for(
int i=0;
i<
n;
i++) {
145 cout <<
"CWB::Toolbox::unionSegments - Error : start must be <= stop - start = " 146 << ilist[
i].start <<
" stop = " << ilist[
i].stop << endl;
150 cout <<
"CWB::Toolbox::unionSegments - Error : start must be positive - start = " 151 << ilist[
i].start << endl;
158 double right_most = -1;
160 for (
int i = 0 ;
i <
n ;
i++) {
161 if (isegs[
i].
start > right_most) {
162 right_most = isegs[
i].
stop;
167 if (isegs[
i].
stop > right_most) {
168 right_most = isegs[
i].
stop;
175 for(
int i=0;
i<cnt+1;
i++)
if(osegs[
i].
stop>0) {seg=osegs[
i];seg.
index=vcnt++;vsegs.push_back(seg);}
194 int size = ilist.size();
199 start[
i] = ilist[
i].start;
200 stop[
i] = ilist[
i].stop;
205 Int_t *
id =
new Int_t[
size];
206 TMath::Sort(size,start,
id,
false);
209 if(start[
id[
i]]<=0) flag=
false;
210 if(start[
id[
i]]>stop[
id[
i]]) flag=
false;
211 if(start[
id[i]]<stop[
id[i-1]]) flag=
false;
214 cout <<
"CWB::Toolbox::invertSegments - Error in segment list (duplicated veto) : " << endl;
215 cout <<
id[i-1]+1 <<
" " << start[
id[i-1]] <<
" " << stop[
id[i-1]] <<
" flag " << flag << endl;
216 cout <<
id[
i]+1 <<
" " << start[
id[
i]] <<
" " << stop[
id[
i]] <<
" flag " << flag << endl;
223 vector<waveSegment>
olist;
233 SEG.
stop = stop[
id[
i]];
234 ctime+=stop[
id[
i]]-start[
id[
i]];
235 olist.push_back(SEG);
262 int size = ilist.size();
267 start[
i] = ilist[
i].start;
268 stop[
i] = ilist[
i].stop;
273 Int_t *
id =
new Int_t[
size];
274 TMath::Sort(size,start,
id,
false);
277 if(start[
id[
i]]<=0) flag=
false;
278 if(start[
id[
i]]>stop[
id[
i]]) flag=
false;
279 if(start[
id[i]]<stop[
id[i-1]]) flag=
false;
282 cout <<
"CWB::Toolbox::invertSegments - Error in segment list (duplicated veto) : " << endl;
283 cout <<
id[i-1]+1 <<
" " << start[
id[i-1]] <<
" " << stop[
id[i-1]] <<
" flag " << flag << endl;
284 cout <<
id[
i]+1 <<
" " << start[
id[
i]] <<
" " << stop[
id[
i]] <<
" flag " << flag << endl;
291 vector<waveSegment>
olist;
300 SEG.
stop = start[
id[0]];
301 olist.push_back(SEG);
302 for(
int i=0;
i<size-1;
i++) {
305 SEG.
stop = start[
id[
i+1]];
307 olist.push_back(SEG);
310 SEG.
start = stop[
id[size-1]];
311 SEG.
stop = 4000000000;
312 olist.push_back(SEG);
332 for (
int i = 0;
i <
n;
i++)
335 f = ((double)
i) / ((double) (n - 1));
336 window[
i] = 0.35875 -
343 for (
int i=0;
i<
n;
i++) norm += pow(window[
i],2);
345 for (
int i=0;i<
n;i++) window[i] /= sqrt(norm);
372 int n1=ilist1.size();
373 int n2=ilist2.size();
378 vector<waveSegment>
olist;
382 while (i<n1 && j<n2) {
383 if (ilist2[j].stop<=ilist1[i].start) j++;
384 else if (ilist1[i].stop <= ilist2[j].start) i++;
386 if (ilist2[j].start < ilist1[i].start) start=ilist1[
i].start;
387 else start = ilist2[
j].start;
388 if (ilist2[j].stop > ilist1[i].stop) stop=ilist1[
i].stop;
389 else stop=ilist2[
j].stop;
390 if (ilist2[j].stop >= ilist1[i].stop) i++;
397 olist.push_back(SEG);
421 vector<waveSegment>
olist;
427 if (!in.good()) {cout <<
"CWB::Toolbox::readSegList - Error Opening File : " << DQF.
file << endl;gSystem->Exit(1);}
433 in.getline(str,1024);
434 if (!in.good())
break;
435 if(str[0] !=
'#') size++;
438 in.clear(ios::goodbit);
439 in.seekg(0, ios::beg);
440 if (size==0) {cout <<
"CWB::Toolbox::readSegList - Error : File " << DQF.
file <<
" is empty" << endl;gSystem->Exit(1);}
448 in.getline(str,1024);
449 if(str[0] ==
'#')
continue;
450 in.seekg(fpos, ios::beg);
452 if(DQF.
c4) in >> dummy >> start[
N] >> stop[
N] >> dummy;
453 else in >> start[
N] >> stop[
N];
454 if (!in.good())
break;
456 in.seekg(fpos+1, ios::beg);
459 if(stop[N]<=start[N]) {
461 cout <<
"CWB::Toolbox::readSegList - Error Ranges : " << start[
N] <<
" " << stop[
N] << endl;
466 cout <<
"CWB::Toolbox::readSegList - Error Max Range : " << N <<
" " << size << endl;
474 Int_t *
id =
new Int_t[
N];
475 TMath::Sort(N,start,
id,
false);
476 for(
int i=1;
i<
N;
i++) {
478 if(start[
id[
i]]<=0) flag=
false;
479 if(start[
id[
i]]>stop[
id[
i]]) flag=
false;
480 if(start[
id[i]]<stop[
id[i-1]]) flag=
false;
483 cout <<
"CWB::Toolbox::readSegList - Error in segment list file (duplicated veto) : " << DQF.
file << endl;
484 cout <<
id[i-1]+1 <<
" " << start[
id[i-1]] <<
" " << stop[
id[i-1]] <<
" flag " << flag << endl;
485 cout <<
id[
i]+1 <<
" " << start[
id[
i]] <<
" " << stop[
id[
i]] <<
" flag " << flag << endl;
499 SEG.
stop = start[
id[0]];
500 olist.push_back(SEG);
501 for(
int i=0;
i<N-1;
i++) {
504 SEG.
stop = start[
id[
i+1]];
506 olist.push_back(SEG);
509 SEG.
start = stop[
id[N-1]];
510 SEG.
stop = 4000000000;
511 olist.push_back(SEG);
513 for(
int i=0;
i<
N;
i++) {
516 SEG.
stop = stop[
id[
i]];
517 ctime+=stop[
id[
i]]-start[
id[
i]];
518 olist.push_back(SEG);
543 vector<waveSegment>
olist;
548 if(DQF[
n].
cat<=dqcat) dqf[ndqf++]=DQF[
n];
551 cout <<
"CWB::Toolbox::readSegList - no CWB_CAT=" << dqcat <<
" files in the list" << endl;
555 vector<waveSegment>*
list =
new vector<waveSegment>[ndqf];
557 for(
int n=0;
n<ndqf;
n++) list[
n]=readSegList(dqf[
n]);
560 for(
int n=1;n<ndqf;n++) olist = mergeSegLists(list[n],olist);
562 for(
int n=0;n<ndqf;n++) list[n].
clear();
583 if((fP = fopen(fName.Data(),
"w")) == NULL) {
584 cout <<
"cannot open output file " << fName.Data() <<
". \n";
587 cout <<
"Write output file : " << fName.Data() << endl;
590 fprintf(fP,
"# seg start stop duration\n");
591 for(
int i=0;
i<(
int)list.size();
i++) {
593 fprintf(fP,
"%-d\t%-d\t%-d\t%-d\n",
597 int(list[
i].stop-list[
i].start));
636 vector<waveSegment>
olist;
637 olist=getJobList(ilist, segLen, segMLS, segEdge);
670 cout << endl <<
"CWB::Toolbox::getJobList - Error : segMLS must be <= segLen" << endl;
671 cout <<
"segMLS = " << segMLS <<
" - segLen = " << segLen << endl << endl;
678 int size = ilist.size();
680 vector<waveSegment>
olist;
687 start = fmod(ilist[
i].start,1)!=0 ?
int(ilist[
i].start+1) : ilist[
i].start;
688 stop = fmod(ilist[
i].stop ,1)!=0 ?
int(ilist[
i].stop) : ilist[
i].stop;
704 olist.push_back(SEG);
710 half=
int(remainder/2);
715 olist.push_back(SEG);
719 olist.push_back(SEG);
724 olist.push_back(SEG);
730 olist.push_back(SEG);
735 for(
int j=0;
j<n-1;
j++) {
739 olist.push_back(SEG);
741 remainder=stop-SEG.
stop;
742 half=
int(remainder/2);
747 olist.push_back(SEG);
751 olist.push_back(SEG);
756 olist.push_back(SEG);
760 printf(
"Toolbox::getJobList : lost livetime after building of the standard job list = %d sec\n",lostlivetime);
784 vector<waveSegment> jobList1=getJobList(cat1List, segLen, segMLS, segEdge);
787 vector<waveSegment> jobList2;
788 for(
int i=0;
i<jobList1.size();
i++) {
791 detSegs_dq2.push_back(jobList1[
i]);
792 detSegs_dq2 = mergeSegLists(detSegs_dq2,cat2List);
794 if(detSegs_ctime<segTHR) {
795 cout <<
"CWB::Toolbox::getJobList : job segment=" << i+1
796 <<
" live time after cat2 : " << detSegs_ctime << endl;
797 cout <<
" segTHR=" << segTHR
798 <<
" sec -> job removed !!!" << endl;
800 jobList2.push_back(jobList1[i]);
866 if(slagSize<1) slagSize=1;
868 if((
int)slagMax>slagSegs) {
869 cout <<
"CWB::Toolbox::makeSlagList : Error - slagMax must be < slagSegs" << endl;
874 cout <<
"CWB::Toolbox::makeSlagList : Error - slagSegs must be positive" << endl;
878 if((slagMax>0)&&(slagMin>slagMax)) {
879 cout <<
"CWB::Toolbox::getSlagList : Error - slagMin must be < slagMax" << endl;
883 if(slagSite!=NULL)
for(
int n=0;
n<(
int)nIFO;
n++) {
884 if(slagSite[
n] >= nIFO) {
885 cout <<
"CWB::Toolbox::getSlagList : Error slagSite - value out of range " << endl;
895 cout << endl <<
"CWB::Toolbox::getSlagList : read slags from file -> " << slagFile << endl;
900 cout <<
"CWB::Toolbox::getSlagList : Error Opening File : " << slagFile << endl;
910 in.getline(str,1024);
911 if(str[0] ==
'#')
continue;
912 in.seekg(fpos, ios::beg);
915 if (!in.good())
break;
917 for(
int j=0;
j<(
int)slagList.size();
j++) {
918 if(
id==slagList[
j].
jobId) {
919 cout <<
"CWB::Toolbox::getSlagList : duplicated slag id " 920 <<
" or wrong format in slag file -> " << slagFile << endl;
925 for(
int n=0;
n<(
int)nIFO;
n++) {
926 in >>
id; SLAG.
slagId.push_back(
id);
928 slagList.push_back(SLAG);
929 if (slagList.size()==
size)
break;
930 if (!in.good())
break;
932 in.seekg(fpos+1, ios::beg);
937 if(size>slagList.size()) {
938 cout <<
"CWB::Toolbox::getSlagList : Error - slagOff+slagSize " << size <<
" is > entries in slagFile " << slagList.size() << endl;
944 cout << endl <<
"CWB::Toolbox::getSlagList : built-in slags" << endl;
948 for(
int n=0;
n<(
int)nIFO;
n++) {
951 for(
int m=0;
m<(
int)ifo.size();
m++)
if((
int)slagSite[
n]==ifo[
m]) unique=
false;
952 if(unique) ifo.push_back(slagSite[
n]);
964 for(
int n=0;
n<(
int)nIFO;
n++) SLAG.
slagId.push_back(0);
965 slagList.push_back(SLAG);
967 for(
int m=(
int)slagMin;
m<=(
int)slagMax;
m++) {
970 int nifo = ifo.size();
972 getSlagList(slags,size,
m,nifo,
id);
974 gRandom->SetSeed(
m+1);
975 while(slags.size()>0) {
976 int k =
int(gRandom->Uniform(0,
int(slags.size())));
977 slags[
k].jobId=jobId++;
978 slagList.push_back(slags[k]);
979 slags.erase(slags.begin()+
k);
984 for(
int j=0;
j<(
int)slagList.size();
j++) {
985 slag SLAG = slagList[
j];
986 slagList[
j].slagId.resize(nIFO);
987 for(
int n=0;
n<(
int)nIFO;
n++) {
988 for(
int m=0;
m<(
int)ifo.size();
m++)
989 if((
int)slagSite[
n]==ifo[
m]) slagList[
j].slagId[
n] = SLAG.
slagId[
m];
995 if((slagSize>slagList.size()-
slagOff)) {
996 cout <<
"CWB::Toolbox::getSlagList : Error - the number of available slags = " 997 << slagList.size() <<
" are less than the requested slag (slagSize) = " << slagSize << endl;
1006 for(
int i=(
int)slagOff;
i<
int(slagOff+slagSize);
i++) {
1017 for(
int j=(
int)slagOff;
j<
int(slagOff+slagSize);
j++){
1021 for (
int n=0;
n<(
int)nIFO;
n++)
if((
k+slagList[
j].slagId[
n])>
slagSegs) check=
false;
1022 for (
int n=0; n<(
int)nIFO; n++)
if((
k+slagList[
j].slagId[n])<0) check=
false;
1028 SLAG.
slagId.push_back(++nseg);
1029 for(
int n=0; n<(
int)nIFO; n++) {
1030 SLAG.
slagId.push_back(slagList[
j].slagId[n]);
1031 SLAG.
segId.push_back(slagList[
j].slagId[n]+
k);
1033 slist.push_back(SLAG);
1056 if(slagSize && ((
int)slagList.size()==
slagSize))
return;
1057 for(
int j=-slagRank;
j<=slagRank;
j++) {
1060 for(
int n=0;
n<(
int)
id.
size();
n++)
if(
j==
id[
n]) unique=
false;
1061 if(!unique)
continue;
1065 for(
int n=0; n<(
int)
id.
size(); n++) m+=abs(
id[n]);
1068 SLAG.
jobId=slagList.size();
1069 SLAG.
slagId.push_back(0);
1070 for(
int n=0; n<(
int)
id.
size(); n++) SLAG.
slagId.push_back(
id[n]);
1071 slagList.push_back(SLAG);
1072 if(slagSize && ((
int)slagList.size()==
slagSize))
return;
1074 id.resize(
id.
size()-1);
1078 getSlagList(slagList, slagSize, slagRank, nifo-1,
id);
1079 id.resize(
id.
size()-1);
1104 if(((
int)slagList.size()>0)&&(slagFile!=
"")) {
1107 if((fP = fopen(slagFile.Data(),
"w")) == NULL) {
1108 cout <<
"CWB::Toolbox::makeSlagList : Error - cannot open file " << slagFile.Data() << endl;
1112 int nIFO = slagList[0].segId.size();
1117 fprintf(fP,
"# Super Lags List - %d jobs\n",(
int)slagList.size());
1119 fprintf(fP,
"# nIFO %13d\n",
int(nIFO));
1128 for(
int j=0;
j<(
int)slagList.size();
j++){
1130 if(slagList[
j].slagId[1]==1) {
1131 fprintf(fP,
"%14d", slagList[
j].slagId[0]);
1137 fprintf(fP,
"%14d", slagList[
j].slagId[0]);
1165 vector<waveSegment> _jobList(jobList.size());
1166 for(
int i=0;
i<(
int)jobList.size();
i++) _jobList[
i].
index=jobList[
i];
1167 return createDagFile(_jobList, condor_dir, label, jobFiles, stage, jobmin, jobmax);
1185 return createDagFile(jobList, condor_dir, label, jobFiles,
"CWB_STAGE_FULL", jobmin, jobmax);
1204 if(jobFiles.size()==0) {
1206 if(gSystem->Getenv(
"CWB_UPARAMETERS_FILE")!=NULL) {
1207 cwb_uparameters_file=
TString(gSystem->Getenv(
"CWB_UPARAMETERS_FILE"));
1208 if(cwb_uparameters_file!=
"")
checkFile(cwb_uparameters_file);
1209 jobFiles.resize(jobList.size());
1210 for(
int n=0;
n<(
int)jobFiles.size();
n++) jobFiles[
n]=cwb_uparameters_file;
1212 cout <<
"CWB::Toolbox::createDagFile : Error - CWB_UPARAMETERS_FILE env is not defined" << endl;
1216 if(jobFiles.size()!=jobList.size()) {
1217 cout <<
"CWB::Toolbox::createDagFile : Error - jobFiles size " << jobFiles.size() <<
1218 "is != jobList size " << jobList.size() << endl;
1223 int end = jobmax<1||jobmax>(
int)jobList.size() ? jobList.size() :
jobmax;
1226 sprintf(ofile,
"%s/%s.dag",condor_dir.Data(),label.Data());
1233 if(jobFiles[
n-1]==
"")
continue;
1234 jID = jobList[
n-1].index;
1236 sprintf(ostring,
"JOB A%i %s/%s.sub",jID,condor_dir.Data(),label.Data());
1237 out << ostring << endl;
1238 sprintf(ostring,
"VARS A%i PID=\"%i\" CWB_UFILE=\"%s\" CWB_STAGE=\"%s\"",
1239 jID,jID,jobFiles[
n-1].Data(),stage.Data());
1240 out << ostring << endl;
1241 if(gSystem->Getenv(
"_USE_OSG")!=NULL) {
1242 sprintf(ostring,
"SCRIPT POST A%i %s/cwb_net_osg_post.sh %i",jID,gSystem->Getenv(
"CWB_SCRIPTS"),jID);
1243 out << ostring << endl;
1245 sprintf(ostring,
"RETRY A%i 3000",jID);
1246 out << ostring << endl;
1268 return createDagFile(slagList, condor_dir, label, jobFiles,
"CWB_STAGE_FULL", jobmin, jobmax);
1287 if(jobFiles.size()==0) {
1289 if(gSystem->Getenv(
"CWB_UPARAMETERS_FILE")!=NULL) {
1290 cwb_uparameters_file=
TString(gSystem->Getenv(
"CWB_UPARAMETERS_FILE"));
1291 if(cwb_uparameters_file!=
"")
checkFile(cwb_uparameters_file);
1292 jobFiles.resize(slagList.size());
1293 for(
int n=0;
n<(
int)jobFiles.size();
n++) jobFiles[
n]=cwb_uparameters_file;
1295 cout <<
"CWB::Toolbox::createDagFile : Error - CWB_UPARAMETERS_FILE env is not defined" << endl;
1299 if(jobFiles.size()!=slagList.size()) {
1300 cout <<
"CWB::Toolbox::createDagFile : Error - jobFiles size " << jobFiles.size() <<
1301 "is != slagList size " << slagList.size() << endl;
1306 int end = jobmax<1||jobmax>(
int)slagList.size() ? slagList.size() :
jobmax;
1309 sprintf(ofile,
"%s/%s.dag",condor_dir.Data(),label.Data());
1316 if(jobFiles[
n-1]==
"")
continue;
1317 jID = slagList[
n-1].jobId;
1319 sprintf(ostring,
"JOB A%i %s/%s.sub",jID,condor_dir.Data(),label.Data());
1320 out << ostring << endl;
1321 sprintf(ostring,
"VARS A%i PID=\"%i\" CWB_UFILE=\"%s\" CWB_STAGE=\"%s\"",
1322 jID,jID,jobFiles[
n-1].Data(),stage.Data());
1323 out << ostring << endl;
1324 sprintf(ostring,
"RETRY A%i 3000",jID);
1325 out << ostring << endl;
1326 if(gSystem->Getenv(
"_USE_OSG")!=NULL) {
1327 sprintf(ostring,
"SCRIPT POST A%i %s/cwb_net_osg_post.sh",jID,gSystem->Getenv(
"CWB_SCRIPTS"));
1328 out << ostring << endl;
1352 sprintf(ofile,
"%s/%s.sub",condor_dir.Data(),label.Data());
1354 sprintf(ofile,
"%s/%s.sub.%s",condor_dir.Data(),label.Data(),ext.Data());
1357 if((fP = fopen(ofile,
"w")) == NULL) {
1358 cout <<
"CWB::Toolbox::createSubFile : Error - cannot open file " << ofile << endl;
1362 fprintf(fP,
"universe = vanilla\n");
1363 fprintf(fP,
"getenv = true\n");
1364 fprintf(fP,
"priority = $(PRI)\n");
1365 fprintf(fP,
"on_exit_hold = ( ExitCode != 0 )\n");
1366 fprintf(fP,
"request_memory = 3000\n");
1367 fprintf(fP,
"executable = cwb.sh\n");
1368 fprintf(fP,
"job_machine_attrs = Machine\n");
1369 fprintf(fP,
"job_machine_attrs_history_length = 5\n");
1370 fprintf(fP,
"requirements = target.machine =!= MachineAttrMachine1 && target.machine =!= MachineAttrMachine2 && target.machine =!= MachineAttrMachine3 && target.machine =!= MachineAttrMachine4 && target.machine =!= MachineAttrMachine5\n");
1371 fprintf(fP,
"environment = CWB_JOBID=$(PID);CWB_UFILE=$(CWB_UFILE);CWB_STAGE=$(CWB_STAGE)\n");
1372 if(condor_tag!=
"")
fprintf(fP,
"accounting_group = %s\n",condor_tag.Data());
1373 fprintf(fP,
"output = %s/$(PID)_%s_$(CWB_STAGE).out\n",out_dir.Data(),label.Data());
1374 fprintf(fP,
"error = %s/$(PID)_%s_$(CWB_STAGE).err\n",err_dir.Data(),label.Data());
1375 if(gSystem->Getenv(
"_USE_OSG")!=NULL) {
1378 fprintf(fP,
"log = %s/condor/%s.log\n",home_dir.Data(),label.Data());
1379 fprintf(fP,
"should_transfer_files = YES\n");
1380 fprintf(fP,
"when_to_transfer_output = ON_EXIT\n");
1381 fprintf(fP,
"transfer_input_files = %s/.rootrc, %s/%s.tgz\n",home_dir.Data(), condor_dir.Data(),label.Data());
1382 fprintf(fP,
"transfer_output_files = %s/output, %s/log\n", label.Data(), label.Data());
1383 fprintf(fP,
"request_memory = 2000\n");
1385 fprintf(fP,
"log = %s/%s.log\n",log_dir.Data(),label.Data());
1387 fprintf(fP,
"notification = never\n");
1407 sprintf(ifile,
"%s/%s.dag",condor_dir.Data(),label.Data());
1410 if(!in.good()) {cout <<
"Error Opening File : " << ifile << endl;gSystem->Exit(1);}
1416 in.getline(istring,1024);
1417 if (!in.good())
break;
1419 TObjString* stoken =(TObjString*)token->At(0);
1420 TString jobLabel = stoken->GetString();
1421 if(jobLabel.CompareTo(
"JOB")!=0)
continue;
1422 stoken =(TObjString*)token->At(1);
1424 jobs.push_back(jobID);
1425 if(token)
delete token;
1431 Int_t *
id =
new Int_t[jobs.size()];
1432 Int_t *jd =
new Int_t[jobs.size()];
1433 for(
int i=0;
i<(
int)jobs.size();
i++) jd[
i]=jobs[
i];
1434 TMath::Sort((
int)jobs.size(),jd,
id,
false);
1435 for(
int i=0;
i<(
int)jobs.size();
i++) jobList.push_back(jd[
id[
i]]);
1450 return getJobBenchmark(ifName, stageID, -1, -1, bench);
1484 vector<float> vbench(2);
1488 TFile *
ifile = TFile::Open(ifName);
1490 cout <<
"Failed to open " << ifName.Data() << endl;
1495 if(ihistory==NULL) {
1496 cout <<
"Error : history is not present!!!" << endl;
1501 char stageName[64];
sprintf(stageName,
"STG:%d-",stageID);
1504 char resName[64];
sprintf(resName,
"RES:%d-",resID);
1507 char factorName[64];
sprintf(factorName,
"FCT:%d-",factorID);
1510 for(
int i=0;
i<log_size;
i++) {
1513 if(!log.Contains(stageName))
continue;
1514 if(resID>=0 && !log.Contains(resName))
continue;
1515 if(factorID>=0 && !log.Contains(factorName))
continue;
1516 if(log.Contains(bench+
TString(
":"))) {
1517 TObjArray*
token = log.Tokenize(
'\n');
1518 for(
int j=0;
j<token->GetEntries();
j++) {
1521 if(line.Contains(bench+
TString(
":"))) {
1524 TObjArray* ltoken = line.Tokenize(
'-');
1525 for(
int k=0;
k<ltoken->GetEntries();
k++) {
1528 TObjArray* stoken = stat.Tokenize(
':');
1529 if(stoken->GetEntries()!=2)
continue;
1534 if(stat_name==bench) {
1535 float value = stat_value.Atof();
1536 if(value>vbench[1]) vbench[1]=
value;
1539 if(stat_name==
"JOB") {
1540 vbench[0] = stat_value.Atoi();
1554 if(vbench[0]==-1) vbench[0] = getJobId(ifName);
1578 UserGroup_t*
uinfo = gSystem->GetUserInfo();
1583 if(gSystem->Getenv(
"HOME_WAT")!=NULL) {
1584 cwb_home_wat=
TString(gSystem->Getenv(
"HOME_WAT"));
1586 if(cwb_home_wat==
"") {
1587 cout <<
"CWB::Toolbox::DAG2LSF : Error : HOME_WAT not defined !!!" << endl;
1593 if(gSystem->Getenv(
"LSF_QUEUE")!=NULL) {
1594 lsf_queue=
TString(gSystem->Getenv(
"LSF_QUEUE"));
1597 cout <<
"CWB::Toolbox::DAG2LSF : Error : LSF_QUEUE not defined !!!" << endl;
1609 cout <<
"CWB::Toolbox::DAG2LSF - Error Opening File : " << dagFile << endl;
1615 lsfFile.ReplaceAll(
".dag",
".lsf");
1619 cout <<
"CWB::Toolbox::DAG2LSF - Error Opening File : " << lsfFile << endl;
1623 out <<
"#!/bin/bash" << endl << endl;
1631 in.getline(istring,1024);
1632 if (!in.good())
break;
1633 if(!
TString(istring).BeginsWith(
"VARS"))
continue;
1635 for(
int j=2;
j<token->GetEntries();
j++) {
1637 if(item.BeginsWith(
"PID=")) {
1638 item.ReplaceAll(
"PID=",
"");
1639 item.ReplaceAll(
"\"",
"");
1643 if(item.BeginsWith(
"CWB_UFILE=")) {
1644 item.ReplaceAll(
"CWB_UFILE=",
"");
1645 item.ReplaceAll(
"\"",
"");
1649 if(item.BeginsWith(
"CWB_STAGE=")) {
1650 item.ReplaceAll(
"CWB_STAGE=",
"");
1651 item.ReplaceAll(
"\"",
"");
1658 char lsf_label[1024];
1659 if(CWB_STAGE==
"CWB_STAGE_FULL") {
1660 sprintf(lsf_label,
"%s",data_label);
1662 sprintf(lsf_label,
"%s_%s",data_label,CWB_STAGE.Data());
1664 sprintf(lsf_cmd,
"bsub -q %s -J A%d -g /%s/%s \ 1665 \\\n -f \"%s > %s/%d_%s.ufile\" \ 1666 \\\n -f \"%s/%s.tgz > %s/%d_%s.tgz\" \ 1667 \\\n -o %d_%s.out -f \"%s/%s/%d_%s.out < %d_%s.out\" \ 1668 \\\n -e %d_%s.err -f \"%s/%s/%d_%s.err < %d_%s.err\" \ 1669 \\\n -f \"%s/%s/%d_%s.tgz < %s/%d_%s.tgz\" \ 1670 \\\n -Ep \"rm %d_%s.out\" \ 1671 \\\n -Ep \"rm %d_%s.err\" \ 1672 \\\n -Ep \"rm %s/%d_%s.tgz\" \ 1673 \\\n \"/bin/bash %s/tools/cwb/scripts/cwb_net_lsf.sh %d %s %s %s %s %s\"",
1674 lsf_queue.Data(), PID, uname.Data(),
data_label,
1675 CWB_UFILE.Data(),
nodedir, PID, lsf_label,
1684 out << lsf_cmd << endl << endl;
1686 if(token)
delete token;
1693 return lsf_job_cnt>0 ? lsfFile :
"";
1710 sprintf(ifile,
"%s/merge_%s.M%d.lst",merge_dir.Data(),label.Data(),
version);
1711 vector<TString> jfname;
1712 return getMergeJobList(ifile,jfname);
1725 vector<TString> jfname;
1726 return getMergeJobList(ifname,jfname);
1742 in.open(ifname.Data());
1744 cout <<
"CWB::Toolbox::getMergeJobList : Error Opening File : " << ifname << endl;
1749 vector<TString> jfname;
1753 in.getline(istring,1024);
1754 if (!in.good())
break;
1756 TObjString* stoken =(TObjString*)token->At(token->GetEntries()-1);
1758 if(jobID==0)
continue;
1759 job.push_back(jobID);
1760 jfname.push_back(istring);
1761 if(token)
delete token;
1766 jobFileList.clear();
1768 Int_t *
id =
new Int_t[job.size()];
1769 Int_t *jd =
new Int_t[job.size()];
1770 for(
int i=0;
i<(
int)job.size();
i++) jd[
i]=job[
i];
1771 TMath::Sort((
int)job.size(),jd,
id,
false);
1772 for(
int i=0;
i<(
int)job.size();
i++) {
1773 jobList.push_back(jd[
id[
i]]);
1774 jobFileList.push_back(jfname[
id[i]]);
1804 if(ilist.size()==0) {cout <<
"CWB::Toolbox::getSlagJobList - Error ilist size=0" << endl;gSystem->Exit(1);}
1805 if(seglen<=0) {cout <<
"CWB::Toolbox::getSlagJobList - Error seglen<=0" << endl;gSystem->Exit(1);}
1808 vector<waveSegment>
jlist;
1810 int start=ilist[0].start;
1811 int stop=ilist[ilist.size()-1].stop;
1812 start=seglen*TMath::Nint(
double(start/seglen));
1813 stop=seglen*TMath::Nint(
double(stop/seglen));
1815 int njob=(stop-
start)/seglen;
1816 for(
int n=0;
n<njob;
n++) {
1817 SEG.
start=start+
n*seglen;
1819 jlist.push_back(SEG);
1836 for(
int i=0;
i<(
int)list.size();
i++) {
1854 for(
int j=0;
j<(
int)slagList.size();
j++)
if(slagList[
j].
jobId==jobid) {SLAG=slagList[
j];
break;}
1880 cout <<
"CWB::Toolbox::getSlagList : dqcat must be CWB_CAT1 or CWB_CAT2 !!!" << endl;
1886 int lsize=islagList.size();
1892 vector<waveSegment> dqList;
1896 vector<slag> oslagList;
1899 if(lsize>0) nIFO=islagList[0].segId.size();
1905 dq1List=readSegList(nDQF, DQF,
CWB_CAT1);
1907 jobList=getSlagJobList(dq1List, segLen);
1908 cout <<
"input list size " << lsize << endl;
1909 cout <<
"jobList size " << jobList.size() << endl;
1913 vector<waveSegment>* ilist =
new vector<waveSegment>[
nDQF];
1914 for(
int i=0;
i<
nDQF;
i++) ilist[
i]=readSegList(DQF[
i]);
1917 for(
int j=0;
j<lsize;
j++) {
1918 if(
j%1000==0)
printf(
"%6d - Rejected slags = %d/%lu\n",
j,nRejected,islagList.size());
1919 SLAG = islagList[
j];
1920 slagID = SLAG.slagId[0];
1922 for(
int n=0;
n<
nIFO;
n++) segID[
n]=SLAG.segId[
n];
1926 for(
int i=0;i<
nDQF;i++) DQF[i].
shift=0.;
1927 setSlagShifts(SLAG, ifos, segLen, nDQF, DQF);
1932 for(
int i=0;i<
nDQF;i++) {
1933 if(DQF[i].
cat<=dqcat) {
1938 if(first) {dqList=ilist[
i];first=
false;}
1939 else dqList = mergeSegLists(ilist[i],dqList);
1948 if(dqList.size()==0)
continue;
1959 for(
int i=0;i<
nDQF;i++) {
1965 if(first) {dq1List=ilist[
i];first=
false;}
1966 else dq1List = mergeSegLists(ilist[i],dq1List);
1974 segList.push_back(SEG);
1975 segList = mergeSegLists(dq1List,segList);
1976 SEG = getMaxSeg(segList);
1982 segList.push_back(SEG);
1983 segList = mergeSegLists(dqList,segList);
1986 MSEG = getMaxSeg(segList);
1991 segLength=getTimeSegList(segList);
1994 livetime1+=segLength;
1995 if(segLength<lenTHR) {
1998 livetime2+=segLength;
1999 oslagList.push_back(SLAG);
2002 printf(
"%6lu - Rejected slags = %d/%lu\n\n",islagList.size(),nRejected,islagList.size());
2003 printf(
"Slag livetime before dq cat%d is approximately %d sec \t= %.2f days\n",dqcat,livetime1,
double(livetime1)/86400.);
2004 printf(
"Slag livetime after dq cat%d is approximately %d sec \t= %.2f days\n",dqcat,livetime2,
double(livetime2)/86400.);
2008 for(
int i=0;i<
nDQF;i++) ilist[i].
clear();
2031 for(
int j=0;
j<(
int)ifos.size();
j++)
if(ifos[
j]==DQF[
i].
ifo) ifoID=
j;
2033 cout <<
"CWB::Toolbox::setSlagShifts - " << DQF[
i].
ifo <<
" is not in the contained into ifos vector" << endl;
2036 if(ifoID>=(
int)SLAG.
segId.size()) {
2037 cout <<
"CWB::Toolbox::setSlagShifts - " << DQF[
i].
ifo <<
" index " << ifoID <<
"is not correct" << endl;
2050 vector<waveSegment> dqList) {
2072 segList.push_back(SEG);
2073 segList = mergeSegLists(dqList,segList);
2074 MSEG = getMaxSeg(segList);
2076 if(segLength<segMLS+2*segEdge) {
2077 cout <<
"CWB::Toolbox::getSegList - Segment length too small : " << segLength << endl;
2085 segList.push_back(SEG);
2096 vector<double> dummy;
2097 return getLiveTime(jobList, dqList, dummy);
2121 if(shiftList.size()==0) shiftList.push_back(0.);
2122 int nshift = shiftList.size();
2125 vector<waveSegment>*
segList =
new vector<waveSegment>[nshift];
2126 vector<waveSegment> mergeList=mergeSegLists(jobList,dqList);
2128 for(
int i=0;
i<nshift;
i++) {
2129 double shift = shiftList[
i];
2130 if(shift==0) {segList[
i]=mergeList;
continue;}
2131 int jsize = jobList.size();
2132 for(
int j=0;
j<jsize;
j++) {
2133 double jstart = jobList[
j].start;
2134 double jstop = jobList[
j].stop;
2135 double jlen = jstop-jstart;
2136 int ksize = mergeList.size();
2137 for(
int k=0;
k<ksize;
k++) {
2138 double kstart = mergeList[
k].start;
2139 double kstop = mergeList[
k].stop;
2140 if((kstop<=jstart)||(kstart>=jstop))
continue;
2146 seg.
start=kstart;seg.
stop=kstop;segList[
i].push_back(seg);
2149 seg.
start=kstart-jlen;seg.
stop=kstop-jlen;segList[
i].push_back(seg);
2151 if((kstart<jstop)&&(kstop>jstop)) {
2152 seg.
start=kstart;seg.
stop=jstop;segList[
i].push_back(seg);
2153 seg.
start=jstart;seg.
stop=kstop-jlen;segList[
i].push_back(seg);
2160 for(
int i=0;
i<nshift;
i++) segList[
i]=unionSegments(segList[
i]);
2163 for(
int i=1;i<nshift;i++) segList[0]=mergeSegLists(segList[0],segList[i]);
2165 double livetime=getTimeSegList(segList[0]);
2166 for(
int i=0;i<nshift;i++) segList[i].
clear();
2184 for(
int i=0;
i<(
int)mdcList.size();
i++) {
2185 char mdc_string[1024]=
"";
2188 sprintf(mdc_string,
"%s",(((TObjString*)token->At(0))->
GetString()).Data());
2189 for(
int j=1;
j<9;
j++)
2190 sprintf(mdc_string,
"%s %s",mdc_string,(((TObjString*)token->At(
j))->
GetString()).Data());
2191 double frTime = (((TObjString*)token->At(9))->
GetString()).Atof();
2192 sprintf(mdc_string,
"%s %10.6f",mdc_string,frTime+mdc_shift);
2194 double mdcTime = (((TObjString*)token->At(10))->
GetString()).Atof();
2195 sprintf(mdc_string,
"%s %10.6f",mdc_string,mdcTime+mdc_shift);
2197 for(
int j=11;
j<15;
j++)
2198 sprintf(mdc_string,
"%s %s",mdc_string,(((TObjString*)token->At(
j))->
GetString()).Data());
2199 for(
int j=15;
j<token->GetEntries();
j++) {
2201 sprintf(mdc_string,
"%s %s",mdc_string,stoken.Data());
2202 for(
int n=0;
n<(
int)ifos.size();
n++) {
2203 if(ifos[
n]==stoken) {
2204 double ifoTime = (((TObjString*)token->At(
j+1))->
GetString()).Atof();
2205 sprintf(mdc_string,
"%s %10.6f",mdc_string,ifoTime+mdc_shift);
2211 mdcList[
i]=mdc_string;
2213 if(token)
delete token;
2232 char sval[32];
sprintf(sval,
"%e",val);
2233 return SetMDCLog(log, pos, sval);
2249 char log_string[1024]=
"";
2252 for(
int n=0;
n<token->GetEntries();
n++) {
2253 if(
n==pos)
sprintf(log_string,
"%s %s",log_string, val.Data());
2254 else sprintf(log_string,
"%s %s",log_string, (((TObjString*)token->At(
n))->
GetString()).Data());
2256 if(token)
delete token;
2271 int size = token->GetEntries();
2287 for(
int n=0;
n<token->GetEntries();
n++) {
2288 if(
n==pos)
return (((TObjString*)token->At(
n))->
GetString()).Data();
2322 if(mlength<=0)
return 0.;
2343 mergeCWBTrees(nthreads, fileList, simulation, odir, label, brms, bvar, bpsm);
2359 cout <<
"CWB::Toolbox::mergeCWBTrees : Error - nthreads must be <= : " <<
MAX_THREADS << endl;
2368 vector<TString> waveList;
2369 vector<TString> mdcList;
2370 vector<TString> liveList;
2371 vector<TString> rmsList;
2372 vector<TString> varList;
2376 if(fileList.size()<=nthreads) {
2378 lSize[0]=fileList.size();
2380 if(fileList.size()%nthreads!=0) {
2381 for(
int i=0;
i<nthreads;
i++) lSize[
i]=
int(fileList.size()/nthreads);
2382 lSize[0]+=fileList.size()%nthreads;
2384 for(
int i=0;
i<nthreads;
i++) lSize[
i]=
int(fileList.size()/nthreads);
2390 for(
int i=0;
i<nthreads;
i++) {
2392 tLabel[
i]=TString::Format(
"%s.T%d",label.Data(),
i);
2394 waveList.push_back(TString::Format(
"%s/wave_%s.root",odir.Data(),tLabel[
i].Data()));
2396 mdcList.push_back(TString::Format(
"%s/mdc_%s.root",odir.Data(),tLabel[
i].Data()));
2398 liveList.push_back(TString::Format(
"%s/live_%s.root",odir.Data(),tLabel[
i].Data()));
2399 if(brms) rmsList.push_back(TString::Format(
"%s/rms_%s.root",odir.Data(),tLabel[
i].Data()));
2400 if(bvar) varList.push_back(TString::Format(
"%s/var_%s.root",odir.Data(),tLabel[
i].Data()));
2405 for(
int i=0;
i<nthreads;
i++) {
2406 mergeParms[
i].threadID =
i;
2407 mergeParms[
i].fileList = fList[
i];
2409 mergeParms[
i].odir =
odir;
2410 mergeParms[
i].label = tLabel[
i];
2411 mergeParms[
i].brms = brms;
2412 mergeParms[
i].bvar = bvar;
2413 mergeParms[
i].bpsm = bpsm;
2417 for(
int i=0;
i<nthreads;
i++) {
2419 thread[
i] =
new TThread(name,
MergeHandle, (
void*) &mergeParms[
i]);
2421 printf(
"Starting Thread : %s\n",name);
2425 for(
int i=0;
i<nthreads;
i++) thread[
i]->Join();
2430 fName=TString::Format(
"wave_%s",label.Data());
2431 if(waveList.size()) {
2432 mergeTrees(waveList,
"waveburst", odir, fName,
true);
2433 for(
int i=0;
i<waveList.size();
i++) gSystem->Exec(
"rm "+waveList[
i]);
2435 fName=TString::Format(
"mdc_%s",label.Data());
2436 if(mdcList.size()) {
2437 mergeTrees(mdcList,
"mdc", odir, fName,
true);
2438 for(
int i=0;
i<waveList.size();
i++) gSystem->Exec(
"rm "+mdcList[
i]);
2440 fName=TString::Format(
"live_%s",label.Data());
2441 if(liveList.size()) {
2442 mergeTrees(liveList,
"liveTime", odir, fName,
true);
2443 for(
int i=0;
i<liveList.size();
i++) gSystem->Exec(
"rm "+liveList[
i]);
2445 fName=TString::Format(
"rms_%s",label.Data());
2446 if(rmsList.size()) {
2447 mergeTrees(rmsList,
"noise", odir, fName,
true);
2448 for(
int i=0;
i<rmsList.size();
i++) gSystem->Exec(
"rm "+rmsList[
i]);
2450 fName=TString::Format(
"var_%s",label.Data());
2451 if(varList.size()) {
2452 mergeTrees(varList,
"variability", odir, fName,
true);
2453 for(
int i=0;
i<varList.size();
i++) gSystem->Exec(
"rm "+varList[
i]);
2457 for(
int i=0;
i<nthreads;
i++)
delete thread[
i];
2476 mergeCWBTrees(fileList, simulation, odir, label, brms, bvar, bpsm);
2511 if (simulation==1) c =
's';
2514 cout <<
"CWB::Toolbox::mergeCWBTrees - Start merging " 2515 << fileList.size() <<
" files (simulation) ..." << endl;
2517 cout <<
"CWB::Toolbox::mergeCWBTrees - Start merging " 2518 << fileList.size() <<
" files (production) ..." << endl;
2520 TChain wav(
"waveburst");
2521 if(!bpsm) wav.SetBranchStatus(
"Psm",
false);
2524 TChain rms(
"noise");
2525 TChain var(
"variability");
2526 TChain liv(
"liveTime");
2533 cout <<
"CWB::Toolbox::mergeCWBTrees - Add file to chain in progress ..." << endl;
2534 for(
int n=0;
n<(
int)fileList.size();
n++) {
2535 sprintf(f,
"%s",fileList[
n].Data());
2570 if (++nfiles%1000==0) cout <<
"CWB::Toolbox::mergeCWBTrees - " << nfiles
2571 <<
"/" << fileList.size() <<
" files" << endl;
2572 if(c==
's') mdc.Add(f);
2573 else { liv.Add(f);
if(brms) rms.Add(f);
if(bvar) var.Add(f);}
2576 sprintf(s,
"%s/wave_%s.root",odir.Data(),label.Data());
2577 cout <<
"CWB::Toolbox::mergeCWBTrees - Merging " << s <<
" in progress ..." << endl;
2588 sprintf(work_dir,
"%s",gSystem->WorkingDirectory());
2594 TFile fwav(s,
"UPDATE");
2595 history->Write(
"history");
2600 sprintf(s,
"%s/mdc_%s.root",odir.Data(),label.Data());
2601 cout <<
"CWB::Toolbox::mergeCWBTrees - Merging " << s <<
" in progress ..." << endl;
2605 TFile fmdc(s,
"UPDATE");
2611 sprintf(s,
"%s/live_%s.root",odir.Data(),label.Data());
2612 cout <<
"CWB::Toolbox::mergeCWBTrees - Merging " << s <<
" in progress ..." << endl;
2616 TFile fliv(s,
"UPDATE");
2621 sprintf(s,
"%s/rms_%s.root",odir.Data(),label.Data());
2622 cout <<
"CWB::Toolbox::mergeCWBTrees - Merging " << s <<
" in progress ..." << endl;
2625 TFile frms(s,
"UPDATE");
2631 sprintf(s,
"%s/var_%s.root",odir.Data(),label.Data());
2632 cout <<
"CWB::Toolbox::mergeCWBTrees - Merging " << s <<
" in progress ..." << endl;
2635 TFile fvar(s,
"UPDATE");
2641 if(history!=NULL)
delete history;
2643 cout <<
"CWB::Toolbox::mergeCWBTrees - Merged files : " << nfiles<<endl;
2668 cout <<
"CWB::Toolbox::mergeTrees - Start merging " 2669 << fileList.size() <<
" files (simulation) ..." << endl;
2671 TChain
tree(treeName);
2678 cout <<
"CWB::Toolbox::mergeTrees - Add file to chain in progress ..." << endl;
2679 for(
int n=0;
n<(
int)fileList.size();
n++) {
2680 sprintf(f,
"%s",fileList[
n].Data());
2686 sprintf(cmd,
"mv %s %s.zombie",f,f);
2694 int estat = gSystem->GetPathInfo(txtfile,&
id,&size,&flags,&mt);
2697 cmd2.ReplaceAll(
".root",
".txt");
2698 gSystem->Exec(cmd2.Data());
2699 cout<<
"Zombie file!!! Renamed root and txt files as :"<<f<<
".zombie"<<endl;
2712 if (++nfiles%1000==0) cout <<
"CWB::Toolbox::mergeTrees - " << nfiles
2713 <<
"/" << fileList.size() <<
" files" << endl;
2716 if(ofName.EndsWith(
".root")) {
2717 sprintf(s,
"%s/%s",odir.Data(),ofName.Data());
2719 sprintf(s,
"%s/%s.root",odir.Data(),ofName.Data());
2721 cout <<
"CWB::Toolbox::mergeTrees - Merging " << s <<
" in progress ..." << endl;
2726 TFile
froot(s,
"UPDATE");
2727 history->Write(
"history");
2730 if(history!=NULL)
delete history;
2732 cout <<
"CWB::Toolbox::mergeTrees - Merged files : " << nfiles<<endl;
2733 cout <<
"CWB::Toolbox::mergeTrees - Zombie files : " << ZombieCnt<<endl;
2749 in.open(ifName.Data(),
ios::in);
2751 cout <<
"CWB::Toolbox::readFileList - Error Opening File : " << ifName.Data() << endl;
2759 in.getline(istring,8192);
2760 if (!in.good())
break;
2761 fileList.push_back(istring);
2777 char ofile_name[1024];
2778 if(ofName[0]!=
'/') {
2779 sprintf(ofile_name,
"%s/%s",gSystem->WorkingDirectory(),ofName.Data());
2781 sprintf(ofile_name,
"%s",ofName.Data());
2783 cout << ofile_name << endl;
2786 if (!out.good()) {cout <<
"CWB::Toolbox::dumpFileList - Error Opening File : " << ofName.Data() << endl;gSystem->Exit(1);}
2788 for (
int i=0;
i<(
int)fileList.size();
i++) {
2789 out << fileList[
i] << endl;
2801 vector<waveSegment> dqList) {
2816 if (jobId<1) {cout <<
"CWB::Toolbox::getSegList - Error : jobId must be > 0" << endl;gSystem->Exit(1);}
2818 vector<waveSegment>
jobList = getJobList(dqList, segLen, segMLS, segEdge);
2820 if (jobId>(
int)jobList.size()) {
2821 cout <<
"CWB::Toolbox::getSegList - Error : jobId is > max jobs " << jobList.size() << endl;
2825 vector<waveSegment>
detSegs(nIFO);
2826 for(
int i=0;
i<
nIFO;
i++) detSegs[
i] = jobList[jobId-1];
2842 char* fileBuffer=NULL;
2849 in.open(ifName.Data(),
ios::in);
2851 cout <<
"CWB::Toolbox::readFile - Error Opening File : " << ifName.Data() << endl;
2855 fileBuffer =
new char[2*size+1];
2856 bzero(fileBuffer,(2*size+1)*
sizeof(
char));
2861 in.getline(istring,8192);
2862 if (!in.good())
break;
2864 int len = strlen(istring);
2866 strncpy(fileBuffer+fileLength,istring,len+1);
2867 fileLength += len+1;
2892 cout <<
"CWB::Toolbox::setCuts - cuts is empty : nothing to do" << endl;
2896 TString ifname = idir+
"/"+ifName;
2899 cout<<
"Opening File : " << ifname.Data() << endl;
2900 TFile
ifile(ifname);
2901 if (ifile.IsZombie()) {
2902 cout <<
"CWB::Toolbox::setCuts - Error opening file " << ifname.Data() << endl;
2905 TTree *
itree = (TTree*)ifile.Get(trname.Data());
2907 cout <<
"CWB::Toolbox::setCuts - tree " << trname
2908 <<
" is not present in file " << ifname.Data() << endl;
2911 Int_t
isize = (Int_t)itree->GetEntries();
2912 cout <<
"tree size : " << isize << endl;
2914 cout <<
"CWB::Toolbox::setCuts - tree " << trname
2915 <<
" is empty in file " << ifname.Data() << endl;
2921 int err =
formula.Compile(cuts.Data());
2923 cout <<
"CWB::Toolbox::setCuts - wrong input cuts " << cuts << endl;
2928 itree->SetEstimate(isize);
2929 itree->Draw(
"Entry$",cuts.Data(),
"goff",
isize);
2930 int nentries = itree->GetSelectedRows();
2931 cout <<
"CWB::Toolbox::setCuts - selected entries : " 2932 << nentries <<
"/" << itree->GetEntries() << endl;
2935 TString ofname = odir+
"/"+ifName;
2936 ofname.ReplaceAll(
".root",
TString(
".C_")+olabel+
".root");
2937 cout <<
"CWB::Toolbox::setCuts - Create cuts file : " << ofname << endl;
2939 if(!overwrite) gSystem->Exit(0);
2942 TFile
ofile(ofname,
"RECREATE");
2943 TTree*
otree = itree->CopyTree(cuts.Data());
2944 otree->SetDirectory(&ofile);
2954 for(
int i=0;
i<stageList->GetSize();
i++) {
2955 TObjString* stageObjString = (TObjString*)stageList->At(
i);
2956 TString stageName = stageObjString->GetString();
2957 char* stage =
const_cast<char*
>(stageName.Data());
2958 if(stageName==
"CUTS") pcuts=history->
GetHistory(stage,const_cast<char*>(
"PARAMETERS"));
2966 sprintf(work_dir,
"%s",gSystem->WorkingDirectory());
2968 TString _cuts = (pcuts!=
"") ? pcuts+
" "+cuts : cuts;
2970 char logmsg[2048];
sprintf(logmsg,
"Apply Selection Cuts : %s",cuts.Data());
2979 TFile
ofile(ofname,
"UPDATE");
2980 history->Write(
"history");
3010 if(irho!=0 && irho!=1) {
3011 cout <<
"CWB::Toolbox::setIFAR - bad irho value : irho must be 0/1" << endl;
3016 cout <<
"CWB::Toolbox::setIFAR - sels is empty : nothing to do" << endl;
3022 if(!exist) gSystem->Exit(0);
3024 vector<double> xrho,xfar;
3025 int xsize = GetStepFunction(farFile, xrho, xfar);
3027 vector<double> trho,tfar;
3028 for(
int i=0;
i<xsize;
i++) {
3030 trho.push_back(xrho[
i]);
3031 tfar.push_back(xfar[i]);
3034 xrho=trho; xfar=tfar;
3036 std::vector<double>::iterator it;
3037 it = xrho.begin(); xrho.insert(it,-DBL_MAX);
3038 it = xfar.begin(); xfar.insert(it,xfar[0]);
3040 xrho.push_back(DBL_MAX/2);
3041 xfar.push_back(xfar[xfar.size()-1]);
3046 for(
int i=1;
i<xsize;
i++) {
3047 double dx = (xrho[
i]-xrho[
i-1])/100000.;
3048 gr.SetPoint(cnt++,xrho[
i]-dx,xfar[
i-1]);
3049 gr.SetPoint(cnt++,xrho[
i]+dx,xfar[
i]);
3052 TString ifname = idir+
"/"+ifName;
3055 cout<<
"Opening File : " << ifname.Data() << endl;
3056 TFile
ifile(ifname);
3057 if (ifile.IsZombie()) {
3058 cout <<
"CWB::Toolbox::setIFAR - Error opening file " << ifname.Data() << endl;
3061 TTree *
itree = (TTree*)ifile.Get(trname.Data());
3063 cout <<
"CWB::Toolbox::setIFAR - tree " << trname
3064 <<
" is not present in file " << ifname.Data() << endl;
3067 Int_t
isize = (Int_t)itree->GetEntries();
3068 itree->SetEstimate(isize);
3069 cout <<
"tree size : " << isize << endl;
3071 cout <<
"CWB::Toolbox::setIFAR - tree " << trname
3072 <<
" is empty in file " << ifname.Data() << endl;
3077 TTreeFormula fsel(
"sels", sels.Data(),
itree);
3078 int err = fsel.Compile(sels.Data());
3080 cout <<
"CWB::Toolbox::setIFAR - wrong input sels " << sels << endl;
3085 TString ofname = odir+
"/"+ifName;
3086 ofname.ReplaceAll(
".root",
TString(
".S_")+olabel+
".root");
3087 cout <<
"CWB::Toolbox::setIFAR - Create sels file : " << ofname << endl;
3089 if(!overwrite) gSystem->Exit(0);
3092 TFile
ofile(ofname,
"RECREATE");
3093 TTree *
otree = (TTree*)itree->CloneTree(0);
3098 bool ifar_exists=
false;
3099 bool bin_exists=
false;
3100 TIter
next(itree->GetListOfBranches());
3101 while ((branch=(TBranch*)
next())) {
3102 if(
TString(
"ifar").CompareTo(branch->GetName())==0) ifar_exists=
true;
3103 if(
TString(
"bin").CompareTo(branch->GetName())==0) bin_exists=
true;
3107 if (ifar_exists) itree->SetBranchAddress(
"ifar",&ifar);
3108 else otree->Branch(
"ifar",&ifar,
"ifar/F");
3109 string* bin =
new string();
3110 if (bin_exists) itree->SetBranchAddress(
"bin",&bin);
3111 else otree->Branch(
"bin",bin);
3114 itree->Draw(
"Entry$",sels.Data(),
"goff",
isize);
3115 double*
entry = itree->GetV1();
3116 int nsel = itree->GetSelectedRows();
3117 cout <<
"CWB::Toolbox::setIFAR - selected entries : " 3118 << nsel <<
"/" << itree->GetEntries() << endl;
3121 itree->SetBranchAddress(
"rho",rho);
3124 itree->GetEntry(
int(entry[
i]));
3125 double far = gr.Eval(rho[irho]);
3126 double xifar = far>0 ? 1./far : -1;
3134 *bin = olabel.Data();
3139 itree->Draw(
"Entry$",TString::Format(
"!(%s)",sels.Data()).Data(),
"goff",
isize);
3140 int _nsel = itree->GetSelectedRows();
3141 double* _entry = itree->GetV1();
3142 cout <<
"CWB::Toolbox::setIFAR - not selected entries : " 3143 << _nsel <<
"/" << itree->GetEntries() << endl;
3145 for(
int i=0;
i<_nsel;
i++) {
3146 itree->GetEntry(
int(_entry[
i]));
3147 if(!ifar_exists) ifar=0;
3148 if(!bin_exists) *bin=
"";
3162 for(
int i=0;
i<stageList->GetSize();
i++) {
3163 TObjString* stageObjString = (TObjString*)stageList->At(
i);
3164 TString stageName = stageObjString->GetString();
3165 char* stage =
const_cast<char*
>(stageName.Data());
3166 if(stageName==
"IFAR") psels=history->
GetHistory(stage,
CCAST(
"PARAMETERS"));
3167 if(stageName==
"IFAR") pfile=history->
GetHistory(stage,
CCAST(
"FARFILE"));
3176 sprintf(work_dir,
"%s",gSystem->WorkingDirectory());
3178 TString _sels = (psels!=
"") ? psels+
"\n"+sels : sels;
3180 char logmsg[2048];
sprintf(logmsg,
"Created ifar on the bin set (%s) specified by the selection : %s",olabel.Data(),sels.Data());
3182 char* buffer = readFile(farFile);
3184 char* farBuffer=NULL;
3186 farBuffer =
new char[strlen(pfile)+strlen(buffer)+farFile.Sizeof()+10];
3187 sprintf(farBuffer,
"%s\n%s\n%s\n",pfile,farFile.Data(),buffer);
3189 farBuffer =
new char[strlen(buffer)+farFile.Sizeof()+10];
3190 sprintf(farBuffer,
"%s\n%s\n",farFile.Data(),buffer);
3194 delete [] farBuffer;
3203 TFile
ofile(ofname,
"UPDATE");
3204 history->Write(
"history");
3227 cout <<
"CWB::Toolbox::setChunk - bad chunk chunk value : chunk must be >0" << endl;
3231 TString ifname = idir+
"/"+ifName;
3234 cout<<
"Opening File : " << ifname.Data() << endl;
3235 TFile
ifile(ifname);
3236 if (ifile.IsZombie()) {
3237 cout <<
"CWB::Toolbox::setChunk - Error opening file " << ifname.Data() << endl;
3240 TTree *
itree = (TTree*)ifile.Get(trname.Data());
3242 cout <<
"CWB::Toolbox::setChunk - tree " << trname
3243 <<
" is not present in file " << ifname.Data() << endl;
3246 Int_t
isize = (Int_t)itree->GetEntries();
3247 itree->SetEstimate(isize);
3248 cout <<
"tree size : " << isize << endl;
3250 cout <<
"CWB::Toolbox::setChunk - tree " << trname
3251 <<
" is empty in file " << ifname.Data() << endl;
3256 TString ofname = odir+
"/"+ifName;
3258 ofname.ReplaceAll(
".root",
TString(
".K_")+schunk+
".root");
3259 cout <<
"CWB::Toolbox::setChunk - Create sels file : " << ofname << endl;
3261 if(!overwrite) gSystem->Exit(0);
3264 TFile
ofile(ofname,
"RECREATE");
3265 TTree *
otree = (TTree*)itree->CloneTree(0);
3270 bool chunk_exists=
false;
3271 TIter
next(itree->GetListOfBranches());
3272 while ((branch=(TBranch*)
next())) {
3273 if(
TString(
"chunk").CompareTo(branch->GetName())==0) chunk_exists=
true;
3276 if (chunk_exists) itree->SetBranchAddress(
"chunk",&chunk);
3277 else otree->Branch(
"chunk",&chunk,
"chunk/I");
3280 itree->Draw(
"Entry$",
"",
"goff",isize);
3281 double*
entry = itree->GetV1();
3282 int nsel = itree->GetSelectedRows();
3283 cout <<
"CWB::Toolbox::setChunk - selected entries : " 3284 << nsel <<
"/" << itree->GetEntries() << endl;
3287 itree->GetEntry(
int(entry[
i]));
3300 for(
int i=0;
i<stageList->GetSize();
i++) {
3301 TObjString* stageObjString = (TObjString*)stageList->At(
i);
3302 TString stageName = stageObjString->GetString();
3303 char* stage =
const_cast<char*
>(stageName.Data());
3304 if(stageName==
"CHUNK") psels=history->
GetHistory(stage,
CCAST(
"PARAMETERS"));
3305 if(stageName==
"CHUNK") pfile=history->
GetHistory(stage,
CCAST(
"FARFILE"));
3313 sprintf(work_dir,
"%s",gSystem->WorkingDirectory());
3315 char logmsg[2048];
sprintf(logmsg,
"Created chunk leaf : %d",chunk);
3324 TFile
ofile(ofname,
"UPDATE");
3325 history->Write(
"history");
3349 cout <<
"CWB::Toolbox::setUniqueEvents - Create unique events file : " << ofwave << endl;
3351 if(!overwrite) gSystem->Exit(0);
3353 TFile* iwroot = TFile::Open(ifwave);
3354 if(iwroot==NULL||!iwroot->IsOpen()) {
3355 cout <<
"CWB::Toolbox::setUniqueEvents - Error : input file " << ifwave <<
" not found" << endl;
3359 TTree* iwtree = (TTree*) iwroot->Get(
"waveburst");
3361 cout <<
"CWB::Toolbox::setUniqueEvents - Error : tree waveburst not found in file " << ifwave << endl;
3365 cout <<
"CWB::Toolbox::setUniqueEvents : number of input events : " << iwtree->GetEntries() << endl;
3368 TFile* owroot =
new TFile(ofwave,
"RECREATE");
3369 if(owroot->IsZombie()) {
3370 cout <<
"CWB::Toolbox::setUniqueEvents - Error : output file " << ofwave <<
" not opened" << endl;
3374 TTree *owtree = (TTree*)iwtree->CloneTree(0);
3382 bool replaceLeaf=
false;
3383 TIter
next(iwtree->GetListOfBranches());
3384 while ((branch=(TBranch*)
next())) {
3385 if (leaf_name.CompareTo(branch->GetName())==0) {
3386 cout <<
"fragment leaf [" << leaf_name <<
"] already applied" << endl;
3390 cout <<
"CWB::Toolbox::setUniqueEvents : Do you want to overwrite the previous leaf ? (y/n) ";
3392 }
while ((strcmp(answer,
"y")!=0)&&(strcmp(answer,
"n")!=0));
3393 if (strcmp(answer,
"y")==0) {
3402 if (replaceLeaf) owtree->SetBranchAddress(
"fsize",&fsize);
3403 else owtree->Branch(
"fsize",&fsize,
"fsize/I");
3406 sprintf(sel,
"Entry$:rho[%i]:time[%i]:factor",pp_irho,nIFO);
3407 int iwsize=iwtree->GetEntries();
3408 iwtree->Draw(sel,
"",
"goff",iwsize);
3409 double*
entry = iwtree->GetV1();
3410 double*
rho = iwtree->GetV2();
3411 double* time = iwtree->GetV3();
3412 double*
factor = iwtree->GetV4();
3415 std::map <float, float> mfactors;
3416 for(
int i=0;
i<iwsize;
i++) {
3417 iwtree->GetEntry(
i);
3418 mfactors[factor[
i]]=factor[
i];
3421 std::map<float, float>::iterator iter;
3422 for (iter=mfactors.begin(); iter!=mfactors.end(); iter++) {
3423 factors.push_back(iter->second);
3429 sprintf(cuts,
"abs(factor-%g)/%g<0.0001",factors[
i],factors[i]);
3430 iwtree->Draw(sel,cuts,
"goff",iwsize);
3431 int size=iwtree->GetSelectedRows();
3434 iwtree->GetEntry(entry[0]);
3439 TMath::Sort(size,time,index,kFALSE);
3441 float rhomax=rho[
j];
3448 if(abs(time1-time2)>0.001 ||
k==size-1) {
3449 iwtree->GetEntry(entry[imax]);
3456 if(rho[j]>rhomax) {rhomax=rho[
j];imax=
j;}
3463 cout <<
" factor id = " << i <<
"\t reconstructed events : " 3464 << iwtree->GetSelectedRows() <<
"\t unique events : " 3465 << unique_evt <<
"\tfactor = " << factors[
i] << endl;
3478 for(
int i=0;
i<stageList->GetSize();
i++) {
3479 TObjString* stageObjString = (TObjString*)stageList->At(
i);
3480 TString stageName = stageObjString->GetString();
3481 char* stage =
const_cast<char*
>(stageName.Data());
3482 if(stageName==
"CUTS") pcuts=history->
GetHistory(stage,const_cast<char*>(
"PARAMETERS"));
3490 sprintf(work_dir,
"%s",gSystem->WorkingDirectory());
3492 TString cuts = (pcuts!=
"") ? pcuts+
" (U)" :
"(U)";
3494 char logmsg[2048];
sprintf(logmsg,
"Apply Selection Cuts : %s",
"(U)");
3503 TFile
ofile(ofwave,
"UPDATE");
3504 history->Write(
"history");
3572 #define YEAR (24.*3600.*365.) 3573 #define TRIALS 2 // number of trials used in the combine algorithm 3574 #define TIME_UNCERTAINTY 0.1 // (sec) used to select entries belonging to the same event 3579 cout <<
"CWB::Toolbox::CombineCBC - Create combine wave output file : " << ofwave << endl;
3581 if(!overwrite) gSystem->Exit(1);
3585 cout <<
"CWB::Toolbox::CombineCBC - Error : ifarthr must be >=0" << endl;
3590 if(msearch!=
"IMBHB" && msearch!=
"BBH") {
3591 cout <<
"CWB::Toolbox::CombineCBC - Error : msearch not defined, valid values are: IMBHB,BBH" << endl;
3594 int msearch_id = (msearch==
"IMBHB") ? 0 : 1;
3596 bool simulation = ifmdc.size()>0 ? true :
false;
3601 if(ifwave.size()!=2) {cout <<
"CWB::Toolbox::CombineCBC - Error in input wave files declaration" << endl;
exit(1);}
3602 for(
int n=0;
n<2;
n++) {
3603 iwroot[
n] = TFile::Open(ifwave[
n]);
3604 if(iwroot==NULL||!iwroot[n]->IsOpen()) {
3605 cout <<
"CWB::Toolbox::CombineCBC - Error : input wave file " << ifwave[
n] <<
" not found" << endl;
3609 iwtree[
n] = (TTree*) iwroot[n]->Get(
"waveburst");
3610 if(iwtree[n]==NULL) {
3611 cout <<
"CWB::Toolbox::CombineCBC - Error : tree waveburst not found in file " << ifwave[
n] << endl;
3615 cout <<
"CWB::Toolbox::CombineCBC : number of input wave events : " << iwtree[
n]->GetEntries()
3616 <<
" ( " << ((n==0)?
"IMBHB":
"BBH") <<
" )" << endl;
3623 if(simulation && ifmdc.size()!=2) {cout <<
"CWB::Toolbox::CombineCBC - Error in input mdc files declaration" << endl;
exit(1);}
3625 for(
int n=0;
n<2;
n++) {
3626 imroot[
n] = TFile::Open(ifmdc[
n]);
3627 if(imroot==NULL||!imroot[n]->IsOpen()) {
3628 cout <<
"CWB::Toolbox::CombineCBC - Error : input mdc file " << ifmdc[
n] <<
" not found" << endl;
3632 imtree[
n] = (TTree*) imroot[n]->Get(
"mdc");
3633 if(imtree[n]==NULL) {
3634 cout <<
"CWB::Toolbox::CombineCBC - Error : tree mdc not found in file " << ifmdc[
n] << endl;
3638 cout <<
"CWB::Toolbox::CombineCBC : number of input mdc events : " << imtree[
n]->GetEntries()
3639 <<
" ( " << ((n==0)?
"IMBHB":
"BBH") <<
" )" << endl;
3649 for(
int n=0;
n<2;
n++) {
3651 tfwave[
n].ReplaceAll(
".root",TString::Format(
"_tmp%d.root",
n));
3652 twroot[
n] =
new TFile(tfwave[
n],
"RECREATE");
3653 if(twroot[n]->IsZombie()) {
3654 cout <<
"CWB::Toolbox::CombineCBC - Error : temporary output wave file " << tfwave[
n] <<
" not opened" << endl;
3667 for(
int n=0;
n<2;
n++) {
3669 tfmdc[
n].ReplaceAll(
".root",TString::Format(
"_tmp%d.root",
n));
3670 tmroot[
n] =
new TFile(tfmdc[
n],
"RECREATE");
3671 if(tmroot[n]->IsZombie()) {
3672 cout <<
"CWB::Toolbox::CombineCBC - Error : temporary output mdc file " << tfmdc[
n] <<
" not opened" << endl;
3686 for(
int n=0;
n<2;
n++) {
3687 bool check_ifar=
false;
3689 TIter
next(iwtree[
n]->GetListOfBranches());
3690 while ((branch=(TBranch*)
next())) {
3691 if(
TString(branch->GetName())==
"search") {
3692 cout <<
"CWB::Toolbox::CombineCBC - Error : 'search' leaf already exist"; gSystem->Exit(1);
3694 if(
TString(branch->GetName())==
"combine") {
3695 cout <<
"CWB::Toolbox::CombineCBC - Error : 'combine' leaf already exist"; gSystem->Exit(1);
3697 if(
TString(branch->GetName())==
"ifar") check_ifar=
true;
3701 cout <<
"CWB::Toolbox::CombineCBC - Error : ifar is not defined in wave root file " << ifwave[
n] << endl;
3706 string*
SEARCH =
new string();
3710 for(
int n=0;
n<2;
n++) {
3711 twtree[
n]->Branch(
"search",SEARCH);
3712 twtree[
n]->Branch(
"combine",&COMBINE,
"combine/I");
3713 twtree[
n]->SetBranchAddress(
"ifar",&IFAR);
3718 vector<double> vGPS;
3719 vector<string> vSEARCH;
3720 vector<float> vIFAR;
3721 vector<int> vCOMBINE;
3725 std::map <float, float> mfactors;
3727 for(
int n=0;
n<2;
n++) {
3728 int size=imtree[
n]->GetEntries();
3729 imtree[
n]->Draw(
"factor",
"",
"goff",size);
3730 double*
factor = imtree[
n]->GetV1();
3732 mfactors[factor[
i]]=factor[
i];
3736 std::map<float, float>::iterator iter;
3737 for (iter=mfactors.begin(); iter!=mfactors.end(); iter++) {
3738 factors.push_back(iter->second);
3740 int nfactor= simulation ? factors.size() : 1;
3746 char wcuts[1024]=
"";
3747 char mcuts[1024]=
"";
3749 sprintf(wcuts,
"(abs(factor-%g)/%g<0.0001)",factors[
i],factors[i]);
3750 sprintf(mcuts,
"(abs(factor-%g)/%g<0.0001)",factors[i],factors[i]);
3752 sprintf(wcuts,
"lag[%d]==0&&slag[%d]==0",nIFO,nIFO);
3754 if(ifarthr>=0) {
char cuts[1024];
sprintf(cuts,
"%s&&(ifar>%e)",wcuts,ifarthr*
YEAR);
strcpy(wcuts,cuts);}
3755 cout << endl <<
" Factor " <<
i+1 <<
"\twave cuts -> " << wcuts << endl;
3756 if(simulation) cout <<
" Factor " <<
i+1 <<
"\tmdc cuts -> " << mcuts << endl;
3764 int cwsize[2]={0,0};
3767 sprintf(wsel,
"Entry$:time[%i]:ifar:frequency[0]",nIFO);
3769 sprintf(wsel,
"Entry$:time[0]:ifar:frequency[0]");
3771 for(
int n=0;
n<2;
n++) {
3772 iwsize[
n]=iwtree[
n]->GetEntries();
3773 iwtree[
n]->Draw(wsel,wcuts,
"goff",iwsize[
n]);
3774 iwentry[
n] = iwtree[
n]->GetV1();
3775 iwtime[
n] = iwtree[
n]->GetV2();
3776 iwifar[
n] = iwtree[
n]->GetV3();
3777 iwfreq[
n] = iwtree[
n]->GetV4();
3778 cwsize[
n] = iwtree[
n]->GetSelectedRows();
3785 int cmsize[2]={0,0};
3787 sprintf(msel,
"Entry$:time[0]");
3788 if(simulation) cout <<
" Factor " <<
i+1 <<
"\tmdc selections-> " << msel << endl << endl;
3790 for(
int n=0;
n<2;
n++) {
3791 imsize[
n]=imtree[
n]->GetEntries();
3792 imtree[
n]->Draw(msel,mcuts,
"goff",imsize[
n]);
3793 imentry[
n] = imtree[
n]->GetV1();
3794 imtime[
n] = imtree[
n]->GetV2();
3795 cmsize[
n] = imtree[
n]->GetSelectedRows();
3800 int size = cwsize[0]+cwsize[1]+cmsize[0]+cmsize[1];
3801 if(size==0) {cout <<
"CWB::Toolbox::CombineCBC - Error : no entries in wave/mdc tree" << endl;
return 1;}
3805 int* ttype =
new int[
size];
3806 int* stype =
new int[
size];
3807 double* time =
new double[
size];
3808 double* ifar =
new double[
size];
3810 string*
log =
new string[
size];
3811 string* ilog =
new string();
3814 for(
int n=0;
n<2;
n++) {
3816 for(
int i=0;
i<cwsize[
n];
i++) {
3817 entry[
k] = iwentry[
n][
i];
3820 time[
k] = iwtime[
n][
i];
3821 ifar[
k] = iwifar[
n][
i];
3822 freq[
k] = iwfreq[
n][
i];
3826 if(!simulation)
continue;
3828 imtree[
n]->SetBranchAddress(
"log",&ilog);
3829 for(
int i=0;
i<cmsize[
n];
i++) {
3830 entry[
k] = imentry[
n][
i];
3833 time[
k] = imtime[
n][
i];
3836 imtree[
n]->GetEntry(entry[k]);
3843 cout << endl <<
"List of uncombined events found in " << msearch <<
" search" << endl << endl;
3844 for(
int k=0;k<
size;k++) {
3845 cout << k+1 <<
"\tsearch -> " << (stype[
k]==0 ?
"IMBHB" :
"BBH")
3846 <<
"\tGPS -> " << std::setprecision(14) << time[
k]
3847 <<
"\tFREQ -> " << std::setprecision(4) << freq[
k]
3848 <<
"\tIFAR -> " << std::setprecision(6) << ifar[
k]/
YEAR <<
" (years)" << endl;
3854 TMath::Sort(size,time,index,kFALSE);
3860 float eifar[2]={0,0};
3861 float efreq[2]={0,0};
3862 string elog[2]={
"",
""};
3863 int ewidx[2]={-1,-1};
3864 int emidx[2]={-1,-1};
3869 int nX[8]={0,0,0,0,0,0,0,0};
3873 double ptime=time[
j];
3877 eifar[stype[
j]]=ifar[
j];
3878 efreq[stype[
j]]=freq[
j];
3882 elog[stype[
j]]=log[
j];
3888 int kstart = size==1 ? 0 : 1;
3890 for(
int k=kstart; k<kstop; k++) {
3905 cout <<
"CWB::Toolbox::CombineCBC - Error : ";
3906 if(simulation) cout <<
"number of reconstructed events per injection is > 2" << endl;
3907 else cout <<
"number of reconstructed events within " <<
TIME_UNCERTAINTY <<
" (sec) is > 2" << endl;
3911 if(esize[1]>ifmdc.size()) {
3912 cout <<
"CWB::Toolbox::CombineCBC - Error : number of injections per event is > " << ifmdc.size() << endl;
3917 if(esize[1]==ifmdc.size()) {
3920 if(elog[0].compare(elog[1])!=0) {
3921 cout << endl <<
"CWB::Toolbox::CombineCBC - Error : mdc log string are different" << endl << endl;
3922 cout <<
"IMBHB MDC log string: " << endl << elog[0].c_str() << endl << endl;
3923 cout <<
"BBH MDC log string: " << endl << elog[1].c_str() << endl << endl;
3927 if(simulation) nMDC++;
3931 iSEARCH=stype[ewidx[ewidx[0]!=-1?0:1]];
3932 if(iSEARCH==0 && efreq[0]<=fthr) {COMBINE=5;IFAR=eifar[0];}
3933 if(iSEARCH==1 && efreq[1]<=fthr) {COMBINE=6;IFAR=eifar[1]/
TRIALS;}
3934 if(iSEARCH==0 && efreq[0] >fthr) {COMBINE=7;IFAR=eifar[0]/
TRIALS;}
3935 if(iSEARCH==1 && efreq[1] >fthr) {COMBINE=8;IFAR=eifar[1];}
3939 if(efreq[0]<=fthr && efreq[1]>fthr) COMBINE=3;
3940 if(efreq[1]<=fthr && efreq[0]>fthr) COMBINE=4;
3941 iSEARCH = eifar[0]>eifar[1] ? 0 : 1;
3942 IFAR = eifar[0]>eifar[1] ? eifar[0]/
TRIALS : eifar[1]/
TRIALS;
3943 if(efreq[0]<=fthr && efreq[1]<=fthr) {COMBINE=1;IFAR=eifar[0];iSEARCH=0;}
3944 if(efreq[0] >fthr && efreq[1] >fthr) {COMBINE=2;IFAR=eifar[1];iSEARCH=1;}
3951 imtree[iSEARCH]->GetEntry(entry[emidx[iSEARCH]]);
3952 tmtree[iSEARCH]->Fill();
3957 *SEARCH = iSEARCH==0 ?
"IMBHB" :
"BBH";
3959 iwtree[iSEARCH]->GetEntry(entry[ewidx[iSEARCH]]);
3960 twtree[iSEARCH]->Fill();
3962 vGPS.push_back(ptime);
3963 vSEARCH.push_back(*SEARCH);
3964 vIFAR.push_back(IFAR);
3965 vCOMBINE.push_back(COMBINE);
3971 esize[0]= 0;esize[1]= 0;
3972 ewidx[0]=-1;ewidx[1]=-1;
3973 emidx[0]=-1;emidx[1]=-1;
3979 eifar[stype[
j]]=ifar[
j];
3980 efreq[stype[
j]]=freq[
j];
3984 elog[stype[
j]]=log[
j];
3992 cout <<
"----------------------------------------------------------------------------------------------" << endl;
3993 cout <<
" factor id = " <<
i+1 << endl;
3994 cout <<
" " <<
" -> injected events IMBHB/BBH : " << cmsize[0] <<
"/" << cmsize[1] << endl;
3995 cout <<
" " <<
" -> reconstructed events IMBHB/BBH : " << cwsize[0] <<
"/" << cwsize[1] << endl;
3996 cout <<
" " <<
" -> combined reconstructed/injected events : " << nWAVE <<
"/" << nMDC << endl;
3997 cout <<
"----------------------------------------------------------------------------------------------" << endl;
3998 cout <<
" Case 1 -> TRIALS=1 IMBHB band (IMBHB BBH), BBH band ( ) = " << nX[0] << endl;
3999 cout <<
" Case 2 -> TRIALS=1 IMBHB band ( ), BBH band (IMBHB BBH) = " << nX[1] << endl;
4000 cout <<
" Case 3 -> TRIALS=2 IMBHB band (IMBHB ), BBH band ( BBH) = " << nX[2] << endl;
4001 cout <<
" Case 4 -> TRIALS=2 IMBHB band ( BBH), BBH band (IMBHB ) = " << nX[3] << endl;
4002 cout <<
" Case 5 -> TRIALS=1 IMBHB band (IMBHB ), BBH band ( ) = " << nX[4] << endl;
4003 cout <<
" Case 6 -> TRIALS=2 IMBHB band ( BBH), BBH band ( ) = " << nX[5] << endl;
4004 cout <<
" Case 7 -> TRIALS=2 IMBHB band ( ), BBH band (IMBHB ) = " << nX[6] << endl;
4005 cout <<
" Case 8 -> TRIALS=1 IMBHB band ( ), BBH band ( BBH) = " << nX[7] << endl;
4006 cout <<
"----------------------------------------------------------------------------------------------" << endl;
4022 cout <<
" Total number of REC/INJ = " << nREC <<
"/" << nINJ << endl;
4028 for(
int n=0;
n<2;
n++) {
4033 for(
int n=0;
n<2;
n++) {
4040 TList *owlist =
new TList;
4041 for(
int n=0;
n<2;
n++) {
4042 twroot[
n] =
new TFile(tfwave[
n],
"READ");
4043 twtree[
n] = (TTree *)twroot[n]->Get(
"waveburst");
4044 owlist->Add(twtree[n]);
4046 TFile* owroot =
new TFile(ofwave,
"RECREATE");
4047 if(owroot->IsZombie()) {
4048 cout <<
"CWB::Toolbox::CombineCBC - Error : output wave file " << ofwave <<
" not opened" << endl;
4053 TTree* owtree = TTree::MergeTrees(owlist);
4054 owtree->SetName(
"waveburst");
4057 for(
int n=0;
n<2;
n++) twroot[
n]->
Close();
4060 TList *omlist =
new TList;
4061 for(
int n=0;
n<2;
n++) {
4062 tmroot[
n] =
new TFile(tfmdc[
n],
"READ");
4063 tmtree[
n] = (TTree *)tmroot[n]->Get(
"mdc");
4064 omlist->Add(tmtree[n]);
4066 TFile* omroot =
new TFile(ofmdc,
"RECREATE");
4067 if(omroot->IsZombie()) {
4068 cout <<
"CWB::Toolbox::CombineCBC - Error : output mdc file " << ofmdc <<
" not opened" << endl;
4072 TTree* omtree = TTree::MergeTrees(omlist);
4074 omtree->SetName(
"mdc");
4077 for(
int n=0;
n<2;
n++) twroot[
n]->
Close();
4080 for(
int n=0;
n<2;
n++) {
4081 char cmd[1024];
sprintf(cmd,
"rm %s",tfwave[
n].Data());
4085 for(
int n=0;
n<2;
n++) {
4086 char cmd[1024];
sprintf(cmd,
"rm %s",tfmdc[
n].Data());
4092 iwroot[msearch_id]->cd();
4093 int isize = iwtree[msearch_id]->GetEntries();
4095 TFile* owroot =
new TFile(ofwave,
"RECREATE");
4096 if(owroot->IsZombie()) {
4097 cout <<
"CWB::Toolbox::CombineCBC - Error : output wave file " << ofwave <<
" not opened" << endl;
4100 TTree* owtree = (TTree*)iwtree[msearch_id]->
CloneTree(0);
4105 string* oSEARCH =
new string();
4107 double* iTIME =
new double[2*
nIFO];
4108 float* iFREQ =
new float[
nIFO];
4109 float* iLAG =
new float[nIFO+1];
4110 float* iSLAG =
new float[nIFO+1];
4112 owtree->SetBranchAddress(
"ifar",&oIFAR);
4113 owtree->Branch(
"combine",&oCOMBINE,
"combine/I");
4114 owtree->Branch(
"search",oSEARCH);
4115 iwtree[msearch_id]->SetBranchAddress(
"ifar",&iIFAR);
4116 iwtree[msearch_id]->SetBranchAddress(
"time",iTIME);
4117 iwtree[msearch_id]->SetBranchAddress(
"frequency",iFREQ);
4118 iwtree[msearch_id]->SetBranchAddress(
"lag",iLAG);
4119 iwtree[msearch_id]->SetBranchAddress(
"slag",iSLAG);
4121 cout <<
"List of combined events found in " << msearch <<
" search" << endl;
4125 iwtree[msearch_id]->GetEntry(
i);
4126 if(iLAG[nIFO]==0 && iSLAG[nIFO]==0) {
4127 for(
int j=0;
j<vGPS.size();
j++) {
4129 *oSEARCH = vSEARCH[
j];
4130 oCOMBINE = vCOMBINE[
j];
4132 cout << k++ <<
"\tcombine -> " << oCOMBINE <<
"\tsearch -> " << oSEARCH->c_str()
4133 <<
"\tGPS -> " << std::setprecision(14) << vGPS[
j]
4134 <<
"\tFREQ -> " << std::setprecision(4) << iFREQ[0]
4135 <<
"\toIFAR -> " << std::setprecision(6) << oIFAR/
YEAR <<
" (years)" << endl;
4164 for(
int i=0;
i<stageList->GetSize();
i++) {
4165 TObjString* stageObjString = (TObjString*)stageList->At(
i);
4166 TString stageName = stageObjString->GetString();
4167 char* stage =
const_cast<char*
>(stageName.Data());
4168 if(stageName==
"CUTS") pcuts=history->
GetHistory(stage,const_cast<char*>(
"PARAMETERS"));
4176 sprintf(work_dir,
"%s",gSystem->WorkingDirectory());
4178 TString cuts = (pcuts!=
"") ? pcuts+
" ("+infos+
")" :
"("+infos+
")";
4181 sprintf(iparms,
"IMBHB_ifwave=%s, BBH_ifwave=%s, fthr=%0.2f",ifwave[0].Data(),ifwave[1].Data(),fthr);
4183 char logmsg[2048];
sprintf(logmsg,
"Combine IMBHB and BBH searches : %s",
"(J)");
4188 for(
int n=0;
n<2;
n++) iwroot[
n]->
Close();
4192 TFile owfile(ofwave,
"UPDATE");
4193 history->Write(
"history");
4196 TFile omfile(ofmdc,
"UPDATE");
4197 history->Write(
"history");
4229 int nIFO = ifos.size();
4235 vector<TString> wifos;
4240 for(
int i=0;
i<(
int)wifos.size();
i++)
if(wifos[
i].CompareTo(VDQF[
j].
ifo)==0) bifo=
false;
4241 for(
int i=0;
i<
nIFO;
i++)
if(ifos[
i].CompareTo(VDQF[
j].ifo)==0) bifo=
false;
4242 if(bifo) wifos.push_back(VDQF[
j].ifo);
4244 for(
int i=0;
i<(
int)wifos.size();
i++) cout << wifos[
i].Data() <<
" ";
4252 int nNET=TMath::Factorial(5);
4253 vector<TString>* xifos =
new vector<TString>[nNET];
4258 for(
int i=0;
i<
nIFO;
i++)
if(ifos[
i].CompareTo(VDQF[
j].
ifo)==0) bifo=
true;
4263 for(
int k=0;
k<(
int)wifos.size();
k++) {
4264 if(wifos[
k].CompareTo(VDQF[
j].ifo)==0) {
4266 strcpy(XDQF[nXDQF].ifo,ifos[
i].Data());
4269 xifos[
nXDQF].push_back(wifos[
k]);
4281 cout <<
"CWB::Toolbox::setVeto - No veto found : setVeto terminated" << endl;
4286 for(
int i=0;
i<(
int)xifos[
j].
size();
i++) cout << xifos[
j][
i].Data() <<
" ";
4292 vector<waveSegment>
jobList = getJobList(cat1List, segLen, segMLS, segEdge);
4296 char ifo_label[10][1024];
4297 char veto_label[10][1024];
4302 for(
int m=0;
m<
nset;
m++)
if(veto_name.CompareTo(veto_label[
m])==0) bnew=
false;
4303 if(bnew)
strcpy(veto_label[nset++],veto_name.Data());
4308 for(
int m=0;m<
nset;m++) {
4309 if(veto_name.CompareTo(veto_label[m])==0) {
4310 if(!
TString(ifo_label[m]).Contains(XDQF[
n].
ifo)) {
4311 sprintf(ifo_label[m],
"%s%s",ifo_label[m],XDQF[
n].ifo);
4316 char veto_file_label[1024]=
"";
4317 for(
int m=0;m<
nset;m++) {
4318 if(strlen(veto_file_label)==0) {
4319 sprintf(veto_file_label,
".V_%s%s",veto_label[m],ifo_label[m]);
4321 sprintf(veto_file_label,
"%s_%s%s",veto_file_label,veto_label[m],ifo_label[m]);
4326 sprintf(veto_file_label,
"%s.root",veto_file_label);
4327 cout <<
"veto file label : " << veto_file_label << endl;
4329 TString efname = odir+
"/"+ifName;
4330 efname.ReplaceAll(
".root",veto_file_label);
4332 if(!overwrite) gSystem->Exit(0);
4334 TString ifname = idir+
"/"+ifName;
4335 TString ofname = odir+
"/"+ifName;
4336 ofname.ReplaceAll(
".root",
".root.tmp.1");
4343 vector<waveSegment> vlist;
4347 for(
int i=0;
i<
nIFO;
i++)
if(ifos[
i].CompareTo(XDQF[
n].
ifo)==0) iIFO=
i;
4348 if (iIFO==-1)
continue;
4361 for(
int k=0;
k<
nIFO;
k++)
if(ifos[
k].CompareTo(VDQF[
j].ifo)==0) RDQF[nRDQF++]=VDQF[
j];
4365 vlist=readSegList(nRDQF, RDQF,
CWB_CAT2);
4373 for(
int k=0;
k<(
int)xifos[
n].
size();
k++)
if(xifos[
n][
k].CompareTo(VDQF[
j].ifo)==0) RDQF[nRDQF++]=VDQF[
j];
4376 vlist=readSegList(nRDQF, RDQF,
CWB_CAT2);
4383 vlist=readSegList(XDQF[
n]);
4385 double vlist_time = getTimeSegList(vlist);
4397 cout << veto_name <<
" list duration " <<
int(vlist_time) <<
" list size " << vlist.size() << endl;
4402 cout<<
"Opening BKG File : " << ifname.Data() << endl;
4403 TFile
ifile(ifname);
4404 if (ifile.IsZombie()) {
4405 cout <<
"CWB::Toolbox::setVeto - Error opening file " << ifname.Data() << endl;
4408 TTree *
itree = (TTree*)ifile.Get(
"waveburst");
4410 cout <<
"CWB::Toolbox::setVeto - tree waveburst is not present in file " << ifname.Data() << endl;
4413 Int_t
isize = (Int_t)itree->GetEntries();
4414 cout <<
"tree size : " << isize << endl;
4416 cout <<
"CWB::Toolbox::setVeto - tree waveburst is empty in file " << ifname.Data() << endl;
4419 itree->SetEstimate(isize);
4420 char selection[1024];
sprintf(selection,
"time[%d]",iIFO);
4421 itree->Draw(selection,
"",
"goff",isize);
4422 double* time = itree->GetV1();
4425 if(TMath::IsNaN(time[
i])) {
4427 cout <<
"CWB::Toolbox::setVeto - tree waveburst file " << ifname.Data() <<
" contains time NaN in ifo=" << ifos[iIFO] <<
" time=" << time[
i] << endl;
4432 Int_t *
id =
new Int_t[
isize];
4433 bool *bveto =
new bool[
isize];
4434 TMath::Sort((
int)isize,time,
id,
false);
4437 SEG.
start=time[
id[isize-1]];
4438 SEG.
stop=time[
id[isize-1]];
4439 vlist.push_back(SEG);
4440 int vsize = vlist.size();
4458 bool replaceVeto=
false;
4459 TIter
next(itree->GetListOfBranches());
4460 while ((branch=(TBranch*)
next())) {
4461 if (veto_name.CompareTo(branch->GetName())==0) {
4462 cout <<
"Veto [" << veto_name <<
"] already applied" << endl;
4466 cout <<
"Do you want to overwrite the previous spurious ? (y/n) ";
4468 }
while ((strcmp(answer,
"y")!=0)&&(strcmp(answer,
"n")!=0));
4469 if (strcmp(answer,
"y")==0) {
4482 TFile* ftmp =
new TFile(ofname,
"RECREATE");
4483 if (ftmp->IsZombie()) {
4484 cout <<
"CWB::Toolbox::setVeto - Error opening file " << ofname.Data() << endl;
4487 TTree *trtmp = (TTree*)itree->CloneTree(0);
4489 TString tVeto = veto_name;tVeto+=
"/b";
4491 trtmp->SetBranchAddress(veto_name.Data(),&bVeto);
4493 trtmp->Branch(veto_name.Data(),&bVeto,tVeto.Data());
4502 cout <<
"Start applying flag to time["<<iIFO<<
"]..." << endl;
4508 int ipc = (
int)((
double)(isize+1)/10.);
if(ipc==0) ipc=1;
4510 double ttime=time[
id[
h]];
4511 for (
int h=0;h<
isize;h++) bveto[h]=0;
4512 for (
int j=0;
j<vsize;
j++) {
4513 while ((ttime<=vlist[
j].
stop) && (h<
isize)) {
4514 if ((ttime>vlist[
j].
start)&&(ttime<=vlist[
j].stop)) {bveto[
id[
h]]=1;count++;}
else bveto[
id[
h]]=0;
4515 if(++h<isize) ttime=time[
id[
h]];
4516 if (h%ipc==0) {cout <<
pc;
if (pc<100) cout <<
" - ";pc+=10;cout.flush();}
4519 cout << pc << endl << endl;
4521 for (
int h=0;h<
isize;h++) {
4530 cout <<
"Writing new ofile "<< ofname.Data() << endl;
4531 Int_t osize = (Int_t)trtmp->GetEntries();
4532 cout <<
"osize : " << osize << endl;
4533 cout <<
"Flagged events: " << count <<
" Percentage: "<< (double)count/osize << endl;
4542 TString ecwb_pparameters_name =
TString(gSystem->Getenv(
"CWB_PPARAMETERS_FILE"));
4543 for(
int i=0;
i<gApplication->Argc()-1;
i++) {
4544 if(
TString(gApplication->Argv(
i)).Contains(
".C")) {
4545 char* parametersBuffer = readFile(gApplication->Argv(
i));
4546 if(parametersBuffer!=NULL) {
4547 if(
TString(gApplication->Argv(
i))==ecwb_pparameters_name) {
4553 delete [] parametersBuffer;
4559 sprintf(work_dir,
"%s",gSystem->WorkingDirectory());
4562 char veto_log[1024];
4563 sprintf(veto_log,
"%s Flagged events: %d/%d - Percentage : %f ",
4564 veto_name.Data(), count, osize, (double)count/osize );
4573 if(
n%2) ofname.ReplaceAll(
".root.tmp.2",
".root.tmp.1");
4574 else ofname.ReplaceAll(
".root.tmp.1",
".root.tmp.2");
4578 for(
int i=0;
i<nNET;
i++) xifos[
i].
clear();
4585 lfname.ReplaceAll(
".root",veto_file_label);
4586 lfname.ReplaceAll(
"wave_",
"live_");
4587 TString ilfname = idir+
"/"+ifName;
4588 ilfname.ReplaceAll(
"wave_",
"live_");
4590 imfname.ReplaceAll(
"live_",
"mdc_");
4592 mfname.ReplaceAll(
"live_",
"mdc_");
4594 ilstfname.ReplaceAll(
"wave_",
"merge_");
4595 ilstfname.ReplaceAll(
".root",
".lst");
4597 lstfname.ReplaceAll(
"live_",
"merge_");
4598 lstfname.ReplaceAll(
".root",
".lst");
4600 cout <<
" rfile : " << rfname.Data() << endl;
4601 cout <<
" ifile : " << ifname.Data() << endl;
4602 cout <<
" ofile : " << ofname.Data() << endl;
4603 cout <<
" lfile : " << lfname.Data() << endl;
4604 cout <<
" mfile : " << mfname.Data() << endl;
4605 cout <<
" ilstfile : " << ilstfname.Data() << endl;
4606 cout <<
" lstfile : " << lstfname.Data() << endl;
4609 sprintf(cmd,
"mv %s %s",ifname.Data(),ofname.Data());
4610 cout << cmd << endl;
4614 estat = gSystem->GetPathInfo(ilfname,&
id,&size,&flags,&mt);
4616 sprintf(cmd,
"cd %s;ln -sf ../%s %s",odir.Data(),ilfname.Data(),lfname.Data());
4617 cout << cmd << endl;
4620 estat = gSystem->GetPathInfo(imfname,&
id,&size,&flags,&mt);
4622 sprintf(cmd,
"cd %s;ln -sf ../%s %s",odir.Data(),imfname.Data(),mfname.Data());
4623 cout << cmd << endl;
4626 estat = gSystem->GetPathInfo(ilstfname,&
id,&size,&flags,&mt);
4628 sprintf(cmd,
"cd %s;ln -sf ../%s %s",odir.Data(),ilstfname.Data(),lstfname.Data());
4629 cout << cmd << endl;
4632 sprintf(cmd,
"rm %s",rfname.Data());
4633 cout << cmd << endl;
4638 TFile fhist(ofname,
"UPDATE");
4639 history->Write(
"history");
4665 return (estat!=0) ? false :
true;
4688 if ((estat!=0)&&(!question)) {
4689 cout <<
"CWB::Toolbox::checkFile - Error - File/Dir \"" << fName.Data() <<
"\" not exist" << endl;
4696 if ((estat==0)&&(question)) {
4700 cout <<
"File/Dir " << fName.Data() <<
" already exist" << endl;
4701 if(message.Sizeof()>1) cout << message.Data() << endl;
4702 cout <<
"Do you want to overwrite it ? (y/n) ";
4704 cout << endl << endl;
4705 }
while ((strcmp(answer,
"y")!=0)&&(strcmp(answer,
"n")!=0));
4706 if (strcmp(answer,
"n")==0) overwrite=
false;
4729 if((estat==0)&&(question==
true)) {
4734 cout <<
"dir \"" << dir.Data() <<
"\" already exist" << endl;
4735 cout <<
"Do you want to remove the files & recreate dir ? (y/n) ";
4738 }
while ((strcmp(answer,
"y")!=0)&&(strcmp(answer,
"n")!=0));
4739 if (strcmp(answer,
"y")==0) {
4740 if(
remove)
sprintf(cmd,
"rm %s/*",dir.Data());
4741 cout << cmd << endl;
4743 sprintf(cmd,
"mkdir -p %s",dir.Data());
4744 cout << cmd << endl;
4747 }
else if((estat==0)&&(question==
false)) {
4749 sprintf(cmd,
"rm %s/*",dir.Data());
4750 cout << cmd << endl;
4753 sprintf(cmd,
"mkdir -p %s",dir.Data());
4754 cout << cmd << endl;
4757 sprintf(cmd,
"mkdir -p %s",dir.Data());
4758 cout << cmd << endl;
4782 if((estat==0)&&(question==
true)) {
4787 sprintf(cmd,
"rm -rf %s",dir.Data());
4788 cout << cmd << endl;
4789 cout <<
"Do you want to remove the dir ? (y/n) ";
4792 }
while ((strcmp(answer,
"y")!=0)&&(strcmp(answer,
"n")!=0));
4793 if (strcmp(answer,
"y")==0) {
4797 }
else if((estat==0)&&(question==
false)) {
4798 sprintf(cmd,
"rm -rf %s",dir.Data());
4799 cout << cmd << endl;
4821 cout << question <<
" (y/n) ";
4824 }
while ((strcmp(answer,
"y")!=0)&&(strcmp(answer,
"n")!=0));
4825 return (strcmp(answer,
"y")==0) ? true :
false;
4841 bool extension = ofname==
"" ? true :
false;
4842 if(ofname==
"") ofname = macro.GetTitle();
4843 if(!ofname.EndsWith(
".C"))
4844 {cout <<
"CWB::Toolbox::addCWBFlags - Error : macro : " << ofname <<
" must have extension .C" << endl;gSystem->Exit(1);}
4845 if(extension) ofname.ReplaceAll(
".C",
"_CWBFlags.C");
4850 if(!out.good()) {cout <<
"CWB::Toolbox::addCWBFlags - Error : Opening File : " << ofname << endl;gSystem->Exit(1);}
4853 out <<
"// -------------------------------------------------------------------------" << endl;
4854 out <<
"// Compiler cWB Flags" << endl;
4855 out <<
"// -------------------------------------------------------------------------" << endl;
4856 if(gSystem->Getenv(
"_USE_HEALPIX")) {
4857 out <<
"#ifndef _USE_HEALPIX" << endl;
4858 out <<
"#define _USE_HEALPIX" << endl;
4859 out <<
"#endif" << endl;
4861 if(gSystem->Getenv(
"_USE_EBBH")) {
4862 out <<
"#ifndef _USE_EBBH" << endl;
4863 out <<
"#define _USE_EBBH" << endl;
4864 out <<
"#endif" << endl;
4866 if(gSystem->Getenv(
"_USE_LAL")) {
4867 out <<
"#ifndef _USE_LAL" << endl;
4868 out <<
"#define _USE_LAL" << endl;
4869 out <<
"#endif" << endl;
4871 out <<
"// -------------------------------------------------------------------------" << endl;
4875 TList* fLines = macro.GetListOfLines();
4878 while ((obj = (TObjString*)
next())) {
4880 out << line.Data() << endl;
4903 TIter
next(liv.GetListOfBranches());
4904 while ((branch=(TBranch*)
next())) {
4905 if (
TString(
"slag").CompareTo(branch->GetName())==0) slagFound=
true;
4909 cout <<
"CWB::Toolbox::getZeroLiveTime : Error - live tree do not contains slag leaf" << endl;
4913 long int ntot = liv.GetEntries();
4914 liv.Draw(
"live",TString::Format(
"lag[%d]==0&&slag[%d]==0",nIFO,nIFO).Data(),
"goff",ntot);
4915 long int nsel = liv.GetSelectedRows();
4916 double*
live = liv.GetV1();
4919 for(
long int i=0;
i<
nsel;
i++) liveZero += live[
i];
4929 int lag_number,
int slag_number,
int dummy) {
4969 std::map <double, std::map <double, int> > QMap;
4974 TIter
next(liv.GetListOfBranches());
4975 while ((branch=(TBranch*)
next())) {
4976 if (
TString(
"slag").CompareTo(branch->GetName())==0) slagFound=
true;
4980 cout <<
"CWB::Toolbox::getLiveTime : Error - live tree do not contains slag leaf" << endl;
4984 int run_max = liv.GetMaximum(
"run");
4987 long int ntot = liv.GetEntries();
4990 TLeaf* leaf = liv.GetLeaf(
"lag");
4992 cout <<
"CWB::Toolbox::getLiveTime : Error - lag leaf is not present" << endl;
4995 branch = leaf->GetBranch();
4997 for (
long int i=0;
i<ntot; ++
i) {
4998 long int entryNumber = liv.GetEntryNumber(
i);
4999 if (entryNumber < 0)
break;
5000 branch->GetEntry(entryNumber);
5001 double val = leaf->GetValue(nIFO);
5002 if (val>lag_max) lag_max = (
int)val;
5005 liv.SetBranchAddress(
"slag",xslag);
5006 liv.SetBranchAddress(
"lag",xlag);
5007 liv.SetBranchAddress(
"run",&xrun);
5008 liv.SetBranchAddress(
"live",&xlive);
5009 liv.SetBranchStatus(
"*",
false);
5010 liv.SetBranchStatus(
"live",
true);
5011 liv.SetBranchStatus(
"lag",
true);
5012 liv.SetBranchStatus(
"slag",
true);
5013 liv.SetBranchStatus(
"run",
true);
5016 int ipc = double(ntot)/10.;
if(ipc==0) ipc=1;
5017 for(
long int i=0;
i<ntot;
i++) {
5020 if(
i%ipc==0) {
if(ntot>100) {cout << pc<<
"%";
if (pc<100) cout <<
" - ";pc+=10;cout.flush();}}
5022 if((lag_number>=0)&&(slag_number>=0)) {
5023 Live = (xlag[
nIFO]==lag_number&&xslag[
nIFO]==slag_number) ? xlive : 0.;
5025 if((lag_number>=0)&&(slag_number==-1)) {
5026 Live = ((xlag[
nIFO]==lag_number)&&!((xlag[nIFO]==0)&&(xslag[
nIFO]==0))) ? xlive : 0.;
5028 if((lag_number==-1)&&(slag_number>=0)) {
5029 Live = ((xslag[
nIFO]==slag_number)&&!((xlag[nIFO]==0)&&(xslag[
nIFO]==0))) ? xlive : 0.;
5031 if((lag_number==-1)&&(slag_number==-1)) {
5032 Live = !((xlag[
nIFO]==0)&&(xslag[nIFO]==0)) ? xlive : 0.;
5041 if(QMap.find(xlag[nIFO])!=QMap.end()) {
5042 if(QMap[xlag[nIFO]].find(xslag[nIFO])!=QMap[xlag[
nIFO]].end()) {
5049 if((lag_number>=0)&&(slag_number>=0)) {
5050 save = (xlag[
nIFO]==lag_number&&xslag[
nIFO]==slag_number) ? save :
false;
5052 if((lag_number>=0)&&(slag_number==-1)) {
5053 save = ((xlag[
nIFO]==lag_number)&&!((xlag[nIFO]==0)&&(xslag[
nIFO]==0))) ? save :
false;
5055 if((lag_number==-1)&&(slag_number>=0)) {
5056 save = ((xslag[
nIFO]==slag_number)&&!((xlag[nIFO]==0)&&(xslag[
nIFO]==0))) ? save :
false;
5058 if((lag_number==-1)&&(slag_number==-1)) {
5059 save = !((xlag[
nIFO]==0)&&(xslag[nIFO]==0)) ? save :
false;
5063 for (
int j=0; j<
nIFO; j++) Qslag[j].
append(xslag[j]);
5070 dlag=xslag[
nIFO]*lag_max+xlag[
nIFO];
5071 Qdlag.
append(sqrt(dlag));
5074 if(pc==100) cout << pc<<
"%";;
5075 cout << endl << endl;
5078 double* dlag =
new double[Qdlag.
size()];
5080 Int_t *
id =
new Int_t[Qdlag.
size()];
5081 TMath::Sort((
int)Qdlag.
size(),dlag,
id,
false);
5086 for (
int j=0;
j<nIFO+1;
j++) {
5094 for (
int j=0;
j<nIFO+1;
j++) {
5121 TString wdir = gSystem->WorkingDirectory();
5122 TSystemDirectory gdir(
"", dir_name);
5123 TList *dfiles = NULL;
5125 void *dir = gSystem->OpenDirectory(dir_name);
5127 const char *file = 0;
5130 while ((file = gSystem->GetDirEntry(dir))) {
5131 dfiles->Add(
new TSystemFile(file, dir_name));
5133 gSystem->FreeDirectory(dir);
5136 dfiles = gdir.GetListOfFiles();
5139 cout <<
"CWB::Toolbox::getFileListFromDir : Error - dir not found !!!" << endl;
5142 TIter dnext(dfiles);
5148 while ((dfile = (TSystemFile*)dnext())) {
5149 fname = dfile->GetName();
5150 sprintf(path,
"%s/%s",dir_name.Data(),fname.Data());
5152 if ((endString!=
"")&&!fname.EndsWith(endString)) fsave=
false;
5153 if ((beginString!=
"")&&!fname.BeginsWith(beginString)) fsave=
false;
5154 if ((containString!=
"")&&!fname.Contains(containString)) fsave=
false;
5155 if(fsave) fileList.push_back(path);
5158 gSystem->ChangeDirectory(wdir);
5168 std::map<int,TString>
5184 std::map<int, TString> fileMap;
5185 for(
int i=0;
i<fileList.size();
i++) {
5186 fileMap[getJobId(fileList[
i])]=fileList[
i];
5209 double LT=0.,LAMBDA;
5210 double normalization=0;
5211 double enormalization=0;
5218 double ChiSquare=0.0;
5219 double dChiSquare=0.0;
5220 int minCountsPerBin=5;
5222 double Tlag0 = Tlag.
data[0];
5226 int min_fa_per_lag=1000000000;
5227 int max_fa_per_lag=0;
5229 if(Wlag[nIFO].data[
j]>0.) {
5230 if(
fabs(Tlag.
data[
j]-Tlag0)/Tlag0>0.1)
continue;
5232 if(fa_per_lag<min_fa_per_lag) min_fa_per_lag=fa_per_lag;
5233 if(fa_per_lag>max_fa_per_lag) max_fa_per_lag=fa_per_lag;
5237 int nbin = max_fa_per_lag-min_fa_per_lag;
5240 TH1I *h1 =
new TH1I(
"h1",
"h1",nbin,min_fa_per_lag,max_fa_per_lag);
5243 if(Wlag[nIFO].data[
j]>0.) {
5244 if(
fabs(Tlag.
data[
j]-Tlag0)/Tlag0>0.1)
continue;
5250 LAMBDA= LT>0 ? NC/LT : 0;
5252 cout <<
"Total Number of bkg coinc.: " << NC <<
" total Live Time: " 5253 << LT <<
" nb = " << LAMBDA << endl;
5257 TCanvas *
canvas =
new TCanvas(
"BackgroundPoissonFit",
"Background Poisson Fit",32,55,750,502);
5258 canvas->Range(-0.75,-1.23433,6.75,8.17182);
5259 canvas->SetBorderSize(1);
5260 canvas->SetFillColor(0);
5263 canvas->SetFrameFillColor(0);
5266 TF1 poisson(
"PoissonIFunction",
PoissonIFunction,min_fa_per_lag,max_fa_per_lag,2);
5267 poisson.FixParameter(0,h1->GetEntries());
5268 poisson.SetParameter(1,h1->GetMean());
5270 TH1D* hfit = (TH1D*)h1->Clone(
"hfit");
5273 for (
int k=1;
k<=hfit->GetNbinsX();
k++) {
5275 if (poisson.Eval(n)>minCountsPerBin) {
5277 dChiSquare=pow((hfit->GetBinContent(
k)-poisson.Eval(n)),2)/poisson.Eval(n);
5278 ChiSquare+=dChiSquare;
5279 cout <<
"NCoinc = " << n <<
" Found = " << hfit->GetBinContent(
k) <<
" Expected = " 5280 << poisson.Eval(n) <<
" Chi2_bin" << n <<
"= " << dChiSquare<<endl;
5285 cout <<
"CWB::Toolbox::doPoissonPlot - Warning !!! - Poisson NDF <=0 " << endl;
5290 double myPLevel=TMath::Prob(ChiSquare,myNDF);
5291 cout<<
"myChiSquare :"<<ChiSquare<<endl;
5292 cout <<
"NDF : " << myNDF << endl;
5293 cout<<
"myPlevel :"<<myPLevel<<endl;
5294 int XMIN=kMaxInt,
XMAX=kMinInt;
5296 while((hfit->GetBinContent(g)==0)&&(g<hfit->GetNbinsX())){XMIN=
g;g++;}
5298 while((hfit->GetBinContent(hfit->GetNbinsX()-
g)==0)&&(g<hfit->GetNbinsX())){XMAX=hfit->GetNbinsX()-
g;g++;}
5301 cout<<XMIN<<
" "<<XMAX<<endl;
5302 hfit->Fit(
"PoissonIFunction",
"R");
5303 TF1* fpois = hfit->GetFunction(
"PoissonIFunction");
5304 fpois->SetFillColor(kGreen);
5305 fpois->SetFillStyle(3002);
5306 fpois->SetLineColor(kGreen);
5307 fpois->SetLineStyle(1);
5308 fpois->SetLineWidth(1);
5309 fpois->SetLineColor(kGreen);
5310 hfit->SetTitle(
"Poisson Fit (Black : Data - Green : Fit)");
5312 hfit->SetLineWidth(4);
5313 hfit->GetXaxis()->SetTitle(
"#False Alarms/lag");
5314 hfit->GetYaxis()->SetTitle(
"#Counts");
5315 hfit->GetXaxis()->SetTitleSize(0.04);
5316 hfit->GetYaxis()->SetTitleSize(0.05);
5317 hfit->GetXaxis()->SetLabelSize(0.04);
5318 hfit->GetYaxis()->SetLabelSize(0.04);
5319 hfit->GetXaxis()->SetTitleOffset(0.97);
5320 hfit->GetYaxis()->SetTitleOffset(1.0);
5321 hfit->GetXaxis()->SetRange(XMIN,XMAX);
5322 gStyle->SetOptFit(0);
5323 gStyle->SetOptStat(0);
5328 chi2 = fpois->GetChisquare();
5329 ndf = fpois->GetNDF()-1;
5330 plevel = TMath::Prob(chi2,ndf);
5331 normalization = fpois->GetParameter(0);
5332 enormalization = fpois->GetParError(0);
5333 mean = fpois->GetParameter(1);
5334 emean = fpois->GetParError(1);
5336 cout <<
"Mean : "<<mean<<endl;
5337 cout <<
"Norm : "<<normalization<<endl;
5338 cout <<
"ChiSquare : " << chi2 << endl;
5339 cout <<
"NDF : " << ndf << endl;
5340 cout <<
"Plevel : " << plevel << endl;
5342 TLegend *legend =
new TLegend(0.597,0.660,0.961,0.921,NULL,
"brNDC");
5343 legend->SetLineColor(1);
5344 legend->SetTextSize(0.033);
5345 legend->SetBorderSize(1);
5346 legend->SetLineStyle(1);
5347 legend->SetLineWidth(1);
5348 legend->SetFillColor(10);
5349 legend->SetFillStyle(1001);
5351 legend->SetTextSize(0.03);
5353 sprintf(entry,
"# lags = %d ",Nlags);
5354 legend->AddEntry(
"",entry,
"");
5357 legend->AddEntry(
"*",entry,
"");
5358 sprintf(entry,
"Total Live time = %d s",
int(LT));
5359 legend->AddEntry(
"*",entry,
"");
5360 sprintf(entry,
"Mean = %.3f +/-%.3f FA/lag",mean,emean);
5361 legend->AddEntry(
"*",entry,
"");
5362 sprintf(entry,
"Chi2 (counts>%d) = %.3f",minCountsPerBin,ChiSquare);
5363 legend->AddEntry(
"*",entry,
"");
5364 sprintf(entry,
"NDF = %d",myNDF);
5365 legend->AddEntry(
"*",entry,
"");
5366 sprintf(entry,
"p-level = %.3f",myPLevel);
5367 legend->AddEntry(
"*",entry,
"");
5372 sprintf(fname,
"%s/Lag_PoissonFit.png",odir.Data());
5373 canvas->Print(fname);
5388 const int nCWB = 35;
5412 "CWB_ROOTLOGON_FILE",
5413 "CWB_PARAMETERS_FILE",
5414 "CWB_UPARAMETERS_FILE",
5415 "CWB_PPARAMETERS_FILE",
5416 "CWB_UPPARAMETERS_FILE",
5421 "CWB_HTML_BODY_PROD",
5422 "CWB_HTML_BODY_SIM",
5429 for(
int i=0;
i<nCWB;
i++) {
5430 if(gSystem->Getenv(ecwb_name[
i])==NULL) {
5434 ecwb_value[
i]=
TString(gSystem->Getenv(ecwb_name[i]));
5440 for(
int i=0;
i<nCWB;
i++) ecwb_csize+=ecwb_name[
i].Sizeof();
5441 for(
int i=0;
i<nCWB;
i++) ecwb_csize+=ecwb_value[
i].Sizeof();
5443 char* iBuffer =
new char[2*ecwb_csize];
5444 bzero(iBuffer,(2*ecwb_csize)*
sizeof(
char));
5448 for(
int i=0;
i<nCWB;
i++) {
5449 len = ecwb_name[
i].Sizeof();
5450 strncpy(iBuffer+iLength,ecwb_name[
i].Data(),len);
5452 iBuffer[iLength]=
'=';
5454 len = ecwb_value[
i].Sizeof();
5455 strncpy(iBuffer+iLength,ecwb_value[
i].Data(),len);
5457 iBuffer[iLength]=0x0a;
5478 TIter
next(itree->GetListOfBranches());
5479 while ((branch=(TBranch*)
next())) {
5480 if (
TString(leaf.Data()).CompareTo(branch->GetName())==0) bleaf=
true;
5503 cout <<
"CWB::Toolbox::makeSpectrum : Error - chuncklen<=0" << endl;
5508 cout <<
"CWB::Toolbox::makeSpectrum : Error - scratchlen<0" << endl;
5512 int scratchsize = scratchlen*x.
rate();
5513 int blocksize = chuncklen*x.
rate();
5514 scratchsize-=scratchsize%2;
5515 blocksize-=blocksize%2;
5516 double df=(double)x.
rate()/(double)(blocksize);
5518 if(x.
size()-2*scratchsize<blocksize) {
5519 cout <<
"CWB::Toolbox::makeSpectrum : Error - data size not enough to produce PSD" << endl;
5523 int loops = (x.
size()-2*scratchsize)/blocksize;
5524 cout <<
"Rate: " << x.
rate() << endl;
5526 double* window =
new double[blocksize];
5527 blackmanharris(window, blocksize);
5535 for (
int n=0;
n<loops;
n++) {
5539 for (
int i=0;
i<blocksize;
i++) y.
data[
i]=x.
data[
i+scratchsize+shift];
5543 for (
int i=0;
i<blocksize;
i+=2) psd[
i/2]+=pow(y.
data[
i],2)+pow(y.
data[
i+1],2);
5546 for (
int i=0;
i<blocksize/2;
i++) psd[
i]=sqrt(psd[
i]/(
double)loops);
5548 for (
int i=0;
i<blocksize/2;
i++) psd[
i]*=sqrt(2/df);
5550 for (
int i=0;
i<blocksize/2;
i++) psd[
i]*=sqrt(1/df);
5576 makeSpectrum(psd, x, chuncklen, scratchlen, oneside);
5581 if (!out.good()) {cout <<
"CWB::Toolbox::makeSpectrum - Error : Opening File : " << ofname.Data() << endl;gSystem->Exit(1);}
5582 for (
int i=0;
i<psd.
size();
i++) {
5584 out << freq <<
" " << psd[
i] << endl;
5610 #define OTIME_LENGHT 616 // minimum time lenght simulation 5611 #define TIME_SCRATCH 184 // time scratch used by the FFT 5612 #define LOW_CUT_FREQ 2.0 // output noise is 0 for freq<LOW_CUT_FREQ 5613 #define FREQ_RES 0.125 // input PSD is resampled with dF=FREQ_RES 5614 #define SRATE 16384. // output noise is produced @ rate=SRATE 5619 double ilenght = u.
size()/u.
rate();
5624 in.open(fName.Data(),
ios::in);
5626 cout <<
"CWB::Toolbox::getSimNoise - Error Opening File : " << fName.Data() << endl;
5633 in.getline(str,1024);
5634 if (!in.good())
break;
5635 if(str[0] !=
'#') size++;
5638 in.clear(ios::goodbit);
5639 in.seekg(0, ios::beg);
5647 if (!in.good())
break;
5648 if(ish.
data[cnt]<=0)
5649 {cout <<
"CWB::Toolbox::getSimNoise - input sensitivity file : " << fName.Data()
5650 <<
" contains zero at frequency : " << ifr.
data[
cnt] <<
" Hz " << endl;gSystem->Exit(1);}
5660 convertSampleRate(ifr,ish,ofr,osh);
5673 for (
int j=0;
j<time_factor;
j++) {
5678 y*=1./sqrt(time_factor);
5687 double am = y.
data[
i];
5700 int scratch=z.
size()/z.
rate()-otime_lenght;
5701 cout <<
"CWB::Toolbox::getSimNoise - scratch : " << scratch <<
" osize : " << z.
size()/z.
rate() << endl;
5702 if (scratch<0) {cout <<
"Error : bad data length -> " << z.
size()/z.
rate() << endl;gSystem->Exit(1);}
5710 random.SetSeed(seed+run);
5712 random.SetSeed(seed+run+1);
5717 cout <<
"CWB::Toolbox::getSimNoise - Error : whitened data size " << u.
size()
5718 <<
" is greater than temporary datat size " << x.
size() << endl;
5723 for(
int i=0;
i<u.
size();
i++) x[
i+iS]=u[
i];
5739 double df=u.rate()/u.size();
5740 int ipad =
fabs(seed)/
df;
5741 if(ipad>u.size()/2-1) ipad=u.size()/2-1;
5742 double fpad = sqrt(pow(u.data[2*ipad],2)+pow(u.data[2*ipad+1],2))*
run;
5743 for(
int i=0;
i<(
int)u.size()/2;
i++) {
5745 if(frequency<=
fabs(seed)) {
5746 u.data[2*
i] = gRandom->Gaus(0,fpad);
5747 u.data[2*i+1] = gRandom->Gaus(0,fpad);
5753 df=u.rate()/u.size();
5760 scratch=(osize-otime_lenght*u.rate())/2;
5761 for (
int i=0;
i<otime_lenght*u.rate();
i++) u.data[
i]=u.data[scratch+
i];
5763 u.resize(otime_lenght*u.rate());
5782 in.open(fName.Data(),
ios::in);
5784 {cout <<
"CWB::Toolbox::ReadDetectorPSD - Error Opening File : " 5785 << fName.Data() << endl;gSystem->Exit(1);}
5787 cout <<
"CWB::Toolbox::ReadDetectorPSD - Read File : " << fName.Data() << endl;
5792 in.getline(str,1024);
5793 if (!in.good())
break;
5794 if(str[0] !=
'#') size++;
5797 in.clear(ios::goodbit);
5798 in.seekg(0, ios::beg);
5806 if (!in.good())
break;
5807 if(ish.
data[cnt]<=0)
5808 {cout <<
"CWB::Toolbox::ReadDetectorPSD - input sensitivity file : " << fName.Data()
5809 <<
" contains zero at frequency : " << ifr.
data[
cnt] <<
" Hz " << endl;gSystem->Exit(1);}
5815 size=
int(fWidth/dFreq);
5820 convertSampleRate(ifr,ish,ofr,osh);
5822 osh.
rate(size*dFreq);
5839 if(iw.
size()<=0) {cout <<
"CWB::Toolbox::convertSampleRate : Error - input size <=0" << endl;gSystem->Exit(1);}
5840 if(iw.
rate()<=0) {cout <<
"CWB::Toolbox::convertSampleRate : Error - input rate <=0" << endl;gSystem->Exit(1);}
5841 if(ow.
size()<=0) {cout <<
"CWB::Toolbox::convertSampleRate : Error - output size <=0" << endl;gSystem->Exit(1);}
5842 if(ow.
rate()<=0) {cout <<
"CWB::Toolbox::convertSampleRate : Error - output rate <=0" << endl;gSystem->Exit(1);}
5846 double idx=1./iw.
rate();
5847 for (
int i=0;
i<(
int)ix.size();
i++) ix[
i]=
i*idx;
5851 double odx=1./ow.
rate();
5852 for (
int i=0;
i<(
int)ox.size();
i++) ox[
i]=
i*odx;
5855 TGraph* grin =
new TGraph(ix.size(), ix.data, iw.
data);
5856 TGraphSmooth* gs =
new TGraphSmooth(
"normal");
5857 TGraph* grout = gs->Approx(grin,
"linear", ox.size(), ox.data);
5860 for (
int i=0;
i<(
int)ox.size();
i++) {
5861 grout->GetPoint(
i,ox[
i],ow[i]);
5886 if(ix.
size()<=0) {cout <<
"CWB::Toolbox::convertSampleRate : Error - input size <=0" << endl;gSystem->Exit(1);}
5887 if(ox.
size()<=0) {cout <<
"CWB::Toolbox::convertSampleRate : Error - output size <=0" << endl;gSystem->Exit(1);}
5890 TGraph* grin =
new TGraph(ix.
size(), ix.
data, iy.
data);
5891 TGraphSmooth* gs =
new TGraphSmooth(
"normal");
5892 TGraph* grout = gs->Approx(grin,
"linear", ox.
size(), ox.
data);
5896 grout->GetPoint(
i,ox[
i],oy[i]);
5918 gRandom->SetSeed(1);
5923 TString efname = odir+
"/"+ifName;
5924 efname.ReplaceAll(
".root",
".N.root");
5926 if(!overwrite) gSystem->Exit(0);
5928 TString ifname = idir+
"/"+ifName;
5929 TString ofname = odir+
"/"+ifName;
5930 ofname.ReplaceAll(
".root",
".root.tmp.1");
5937 cout<<
"Opening BKG File : " << ifname.Data() << endl;
5938 TFile
ifile(ifname);
5939 if (ifile.IsZombie()) {
5940 cout <<
"CWB::Toolbox::setMultiplicity - Error opening file " << ifname.Data() << endl;
5943 TTree *
itree = (TTree*)ifile.Get(
"waveburst");
5945 cout <<
"CWB::Toolbox::setMultiplicity - tree waveburst is not present in file " << ifname.Data() << endl;
5948 Int_t tsize = (Int_t)itree->GetEntries();
5949 cout <<
"tree size : " << tsize << endl;
5951 cout <<
"CWB::Toolbox::setMultiplicity - tree waveburst is empty in file " << ifname.Data() << endl;
5954 itree->SetEstimate(tsize);
5955 char selection[1024];
sprintf(selection,
"time[%d]:Entry$",
n);
5958 int isize = itree->Draw(selection,cut,
"goff",tsize);
5959 double* time = itree->GetV1();
5960 double*
entry = itree->GetV2();
5962 Int_t *
id =
new Int_t[
isize];
5963 TMath::Sort((
int)isize,time,
id,
false);
5977 TIter
next(itree->GetListOfBranches());
5978 while((branch=(TBranch*)
next())) {
5979 if(
TString(
"Msize").CompareTo(branch->GetName())==0) {
5980 cout <<
"Multiplicity already applied" << endl;
5984 cout <<
"Do you want to overwrite the previous values ? (y/n) ";
5986 }
while((strcmp(answer,
"y")!=0)&&(strcmp(answer,
"n")!=0));
5987 if(strcmp(answer,
"n")==0) {
5999 int* Msize =
new int[
nIFO];
6000 int* Mid =
new int[
nIFO];
6001 UChar_t* Mm =
new UChar_t[
nIFO];
6003 char cMsize[32];
sprintf(cMsize,
"Msize[%1d]/I",nIFO);
6004 char cMid[32];
sprintf(cMid,
"Mid[%1d]/I",nIFO);
6005 char cMm[32];
sprintf(cMm,
"Mm[%1d]/b",nIFO);
6007 TFile* ftmp =
new TFile(ofname,
"RECREATE");
6008 if (ftmp->IsZombie()) {
6009 cout <<
"CWB::Toolbox::setMultiplicity - Error opening file " << ofname.Data() << endl;
6012 TTree *trtmp = (TTree*)itree->CloneTree(0);
6015 trtmp->SetBranchAddress(
"Msize",Msize);
6016 trtmp->SetBranchAddress(
"Mid",Mid);
6017 trtmp->SetBranchAddress(
"Mm",Mm);
6019 trtmp->Branch(
"Msize",Msize,cMsize);
6020 trtmp->Branch(
"Mid",Mid,cMid);
6021 trtmp->Branch(
"Mm",Mm,cMm);
6030 cout <<
"Start applying flag to time["<<
n<<
"]..." << endl;
6035 int ipc = (
int)((
double)(isize+1)/10.);
6036 int *oMsize =
new int[tsize];
6037 int *oMid =
new int[tsize];
6038 bool *oMm =
new bool[tsize];
6039 for(
int j=0;
j<tsize;
j++) {
6047 int ientry=
int(entry[
id[0]]);
6048 oMsize[ientry]=xMsize; oMid[ientry]=xMid; oMm[ientry]=
false;
6049 mlist.push_back(
id[0]);
6051 ientry=
int(entry[
id[
j]]);
6054 if((time[
id[j]]-time[
id[j-1]])<Tgap) {
6056 mlist.push_back(
id[j]);
6058 int xMm=mlist[
int(gRandom->Uniform(0,mlist.size()))];
6059 oMm[
int(entry[xMm])]=
true;
6060 for(
int m=0;
m<(
int)mlist.size();
m++) oMsize[
int(entry[mlist[
m]])]=xMsize;
6063 mlist.push_back(
id[j]);
6067 if(ipc!=0)
if(j%ipc==0) {cout <<
pc;
if (pc<100) cout <<
" - ";pc+=10;cout.flush();}
6069 cout << pc << endl << endl;
6071 int* iMsize =
new int[
nIFO];
6072 int* iMid =
new int[
nIFO];
6073 UChar_t* iMm =
new UChar_t[
nIFO];
6075 itree->SetBranchAddress(
"Msize",iMsize);
6076 itree->SetBranchAddress(
"Mid",iMid);
6077 itree->SetBranchAddress(
"Mm",iMm);
6080 for (
int j=0;
j<tsize;
j++) {
6090 if(Mm[
n]) nevt_master++;
6112 cout <<
"Writing new ofile "<< ofname.Data() << endl;
6113 Int_t osize = (Int_t)trtmp->GetEntries();
6114 cout <<
"osize : " << osize << endl;
6115 cout <<
"Master events: " << nevt_master <<
" Percentage: "<< 100.*double(nevt_master)/double(tsize) << endl;
6117 if((history!=NULL)&&(
n==nIFO-1)) {
6124 sprintf(work_dir,
"%s",gSystem->WorkingDirectory());
6126 char hTgap[1024];
sprintf(hTgap,
"Tgap = %f",Tgap);
6135 if(
n%2) ofname.ReplaceAll(
".root.tmp.2",
".root.tmp.1");
6136 else ofname.ReplaceAll(
".root.tmp.1",
".root.tmp.2");
6144 lfname.ReplaceAll(
".root",
".N.root");
6145 lfname.ReplaceAll(
"wave_",
"live_");
6146 TString ilfname = idir+
"/"+ifName;
6147 ilfname.ReplaceAll(
"wave_",
"live_");
6151 mfname.ReplaceAll(
".root",
".N.root");
6152 mfname.ReplaceAll(
"wave_",
"merge_");
6153 mfname.ReplaceAll(
".root",
".lst");
6154 TString imfname = idir+
"/"+ifName;
6155 imfname.ReplaceAll(
"wave_",
"merge_");
6156 imfname.ReplaceAll(
".root",
".lst");
6158 cout <<
" rfile " << rfname.Data() << endl;
6159 cout <<
" ifile " << ifname.Data() << endl;
6160 cout <<
" ofile " << ofname.Data() << endl;
6161 cout <<
" lfile " << lfname.Data() << endl;
6162 cout <<
" mfile " << mfname.Data() << endl;
6165 sprintf(cmd,
"mv %s %s",ifname.Data(),ofname.Data());
6166 cout << cmd << endl;
6169 sprintf(cmd,
"cd %s;ln -s ../%s %s",odir.Data(),ilfname.Data(),lfname.Data());
6170 cout << cmd << endl;
6173 sprintf(cmd,
"cd %s;ln -s ../%s %s",odir.Data(),imfname.Data(),mfname.Data());
6174 cout << cmd << endl;
6176 sprintf(cmd,
"rm %s",rfname.Data());
6177 cout << cmd << endl;
6182 TFile fhist(ofname,
"UPDATE");
6209 gRandom->SetSeed(1);
6212 float* iWslag =
new float[nIFO+1];
6213 W.
fChain->SetBranchAddress(
"slag",iWslag);
6218 for(
int i=0;
i<(
int)Tlag.
size();
i++) liveTot+=Tlag[
i];
6230 for(
int j=0;
j<
n;
j++) {
6231 if(slag[
j]==islag) {save=
false;
j=
n;}
6233 if(save) slag.
append(islag);
6235 int ilag=
int(Wlag[nIFO].data[i]);
6238 for(
int j=0;
j<
n;
j++) {
6239 if(lag[
j]==ilag) {save=
false;
j=
n;}
6241 if(save) lag.
append(ilag);
6248 int* slag2id =
new int[slag.
max()+1];
6249 int* lag2id =
new int[lag.
max()+1];
6250 for(
int i=0;
i<(
int)slag.
size();
i++) slag2id[slag[
i]]=
i;
6254 int** lag2live = (
int**)malloc((slag.
size())*
sizeof(
int*));
6255 for (
int i=0;
i<(
int)slag.
size();
i++) lag2live[
i] =
new int[lag.
size()];
6259 int ilag=
int(Wlag[nIFO].data[i]);
6260 lag2live[slag2id[
islag]][lag2id[
ilag]]=Tlag[
i];
6264 bool** lagRej = (
bool**)malloc((slag.
size())*
sizeof(
bool*));
6265 for (
int i=0;
i<(
int)slag.
size();
i++) lagRej[
i] =
new bool[lag.
size()];
6268 int ntrg = wav.GetEntries();
6270 int *
id =
new int[
ntrg];
6271 int *wslag =
new int[
ntrg];
6272 int *wlag =
new int[
ntrg];
6274 double** time = (
double**)malloc(nIFO*
sizeof(
double*));
6275 for (
int i=0;
i<
nIFO;
i++) time[
i] =
new double[ntrg];
6276 bool *trgMaster =
new bool[
ntrg];
6277 for (
int i=0;
i<
ntrg;
i++) trgMaster[
i] =
false;
6278 bool *wrej =
new bool[
ntrg];
6279 for (
int i=0;
i<
ntrg;
i++) wrej[
i] =
false;
6283 if(!Wsel[
i])
continue;
6287 for(
int j=0;
j<
nIFO;
j++)
if(!TMath::IsNaN(W.
time[
j])) time[
j][isel]=W.
time[
j];
else skip=
true;
6290 wslag[isel]=iWslag[
nIFO];
6300 vector<double> livlist;
6303 if(isel==0)
continue;
6304 TMath::Sort((
int)isel,time[
i],
id,
false);
6305 int islag = slag2id[wslag[
id[0]]];
6306 int ilag = lag2id[wlag[
id[0]]];
6308 livlist.push_back(live);
6309 idlist.push_back(
id[0]);
6310 for (
int j=1;
j<isel;
j++) {
6311 if(wrej[
id[
j]])
continue;
6312 if((time[i][
id[j]]-time[i][
id[j-1]])>Tgap) {
6320 int M=
int(gRandom->Uniform(0,(
int)idlist.size()));
6321 trgMaster[idlist[
M]]=
true;
6322 for(
int m=0;
m<(
int)idlist.size();
m++) {
6325 int islag = slag2id[wslag[idlist[
m]]];
6326 int ilag = lag2id[wlag[idlist[
m]]];
6327 if(!trgMaster[idlist[
m]]) {
6328 if(!lagRej[islag][ilag]) rejLive+=livlist[
m];
6330 wrej[idlist[
m]]=
true;
6331 Wsel[entry[idlist[
m]]]=2;
6338 int islag = slag2id[wslag[
id[
j]]];
6339 int ilag = lag2id[wlag[
id[
j]]];
6341 livlist.push_back(live);
6342 idlist.push_back(
id[j]);
6351 int M=
int(gRandom->Uniform(0,(
int)idlist.size()));
6352 trgMaster[idlist[
M]]=
true;
6353 for(
int m=0;
m<(
int)idlist.size();
m++) {
6356 int islag = slag2id[wslag[idlist[
m]]];
6357 int ilag = lag2id[wlag[idlist[
m]]];
6358 if(!trgMaster[idlist[
m]]) {
6359 if(!lagRej[islag][ilag]) rejLive+=livlist[
m];
6361 wrej[idlist[
m]]=
true;
6362 Wsel[entry[idlist[
m]]]=2;
6372 for (
int j=0;
j<isel;
j++) {
6373 if(wrej[
j])
continue;
6374 int islag = slag2id[wslag[
j]];
6375 int ilag = lag2id[wlag[
j]];
6376 if(lagRej[islag][ilag]) {
6409 for (
int i=0;
i<
nIFO;
i++)
delete [] time[
i];
6413 for (
int i=0;
i<(
int)slag.
size();
i++)
delete [] lag2live[
i];
6417 delete [] trgMaster;
6421 double rate = (isel-rejTrg)/(liveTot-rejLive);
6422 double erate = sqrt(isel-rejTrg)/(liveTot-rejLive);
6430 vector<double>& ex, vector<double>& ey) {
6432 GetStepFunction(
"", fName, 0, x, y, ex, ey);
6438 x.erase(x.begin()+size-1);
6439 y.erase(y.begin()+size-1);
6440 ex.erase(ex.begin()+size-1);
6441 ey.erase(ey.begin()+size-1);
6446 ex.erase(ex.begin());
6447 ey.erase(ey.begin());
6455 vector<double>&
x, vector<double>&
y, vector<double>& ex, vector<double>& ey) {
6484 if((option!=
"XMIN")&&(option!=
"XMAX")&&
6485 (option!=
"YMIN")&&(option!=
"YMAX")&&
6486 (option!=
"Y")&&(option!=
"EX")&&(option!=
"EY")) {
6487 cout<<
"GetLinearInterpolation : Error : option --value bad value -> " 6488 <<
"must be y/xmin/xmax/ymin/ymax"<<endl;
6491 if(option.BeginsWith(
"E")) {ferror=
true;option.Remove(0,1);}
6494 double xmin=DBL_MAX;
6495 double ymin=DBL_MAX;
6496 double xmax=-DBL_MAX;
6497 double ymax=-DBL_MAX;
6499 double exmin=DBL_MAX;
6500 double eymin=DBL_MAX;
6501 double exmax=-DBL_MAX;
6502 double eymax=-DBL_MAX;
6505 in.open(fName.Data());
6506 if (!in.good()) {cout <<
"GetGraphValue : Error Opening File : " << fName << endl;
exit(1);}
6511 in.getline(line,1024);
6512 if (in.eof())
break;
6513 std::stringstream linestream(line);
6515 if(linestream >> ix >> iy >> iex >> iey) {
6518 if(iex<exmin) exmin=iex;
6519 if(iey<eymin) eymin=iey;
6520 if(iex>exmax) exmax=iex;
6521 if(iey>eymax) eymax=iey;
6523 linestream.str(line);
6525 if(!(linestream >> ix >> iy)) {
6526 cout <<
"GetGraphValue : Wrong line format : must be 'x y' " << endl;
6529 ex.push_back(0); ey.push_back(0);
6530 exmin=0;eymin=0;exmax=0;eymax=0;
6533 if(!(linestream >> ix >> iy)) {
6534 cout <<
"GetGraphValue : Wrong line format : must be 'x y' " << endl;
6537 ex.push_back(0); ey.push_back(0);
6538 exmin=0;eymin=0;exmax=0;eymax=0;
6542 if(ix<xmin) xmin=ix;
6543 if(iy<ymin) ymin=iy;
6544 if(ix>xmax) xmax=ix;
6545 if(iy>ymax) ymax=iy;
6549 if(option==
"XMIN")
return ferror ? exmin : xmin;
6550 if(option==
"XMAX")
return ferror ? exmax : xmax;
6551 if(option==
"YMIN")
return ferror ? eymin : ymin;
6552 if(option==
"YMAX")
return ferror ? eymax : ymax;
6554 if(x.size()==0) cout <<
"GetGraphValue : Error - input file empty" << endl;
6557 int size = x.size();
6561 TMath::Sort(size,X.
data,I.
data,
false);
6565 std::vector<double>::iterator it;
6566 it = x.begin(); x.insert(it,-DBL_MAX);
6567 it = y.begin(); y.insert(it,y[I[0]]);
6568 it = ex.begin(); ex.insert(it,ex[I[0]]);
6569 it = ey.begin(); ey.insert(it,ey[I[0]]);
6576 TMath::Sort(size,X.
data,I.data,
false);
6578 x.push_back(x[I[size-1]]); y.push_back(-DBL_MAX);
6579 ex.push_back(ex[I[size-1]]); ey.push_back(ey[I[size-1]]);
6580 I.resize(size+1); I[
size]=I[size-1]; size++;
6583 if(V<x[I[0]])
return y[I[0]];
6584 if(V>=x[I[size-1]])
return y[I[size-1]];
6585 for(
int i=1;i<size;i++) if((V>=x[I[
i-1]]) && (V<x[I[
i]]))
return y[I[i-1]];
6588 if(V<x[I[0]])
return ex[I[0]];
6589 if(V>=x[I[size-1]])
return ex[I[size-1]];
6590 for(
int i=1;i<size;i++) if((V>=x[I[
i-1]]) && (V<x[I[
i]]))
return ex[I[i-1]];
6591 }
else if(option==
"Y") {
6592 if(V<x[I[0]])
return ey[I[0]];
6593 if(V>=x[I[size-1]])
return ey[I[size-1]];
6594 for(
int i=1;i<size;i++) if((V>=x[I[
i-1]]) && (V<x[I[
i]]))
return ey[I[i-1]];
6608 int R=1;
while (R < 2*(
int)w.
rate()) R*=2;
6610 if(R==2*w.
rate())
return;
6626 fprintf(stdout,
"--------------------------------\n");
6628 fprintf(stdout,
"--------------------------------\n");
6640 gSystem->Exec(
"date");
6642 FILE *
f = fopen(Form(
"/proc/%d/statm", gSystem->GetPid()),
"r");
6645 sscanf(s.Data(),
"%ld %ld", &total, &rss);
6646 cout << str.Data() <<
" virtual : " << total * 4 / 1024 <<
" (mb) rss : " << rss * 4 / 1024 <<
" (mb)" << endl;
6662 gRandom->SetSeed(0);
6663 int rnID =
int(gRandom->Rndm(13)*1.e9);
6664 UserGroup_t*
uinfo = gSystem->GetUserInfo();
6667 gSystem->Exec(
TString(
"mkdir -p /dev/shm/")+uname);
6669 char fName[1024]=
"";
6671 sprintf(fName,
"/dev/shm/%s/%s_%d.%s",uname.Data(),label.Data(),rnID,extension.Data());
6686 std::vector<std::string> vec(ilist.size());
6687 for(
int i=0;
i<ilist.size();
i++) vec[
i]=ilist[
i].Data();
6688 std::sort( vec.begin(), vec.end() ) ;
6689 vector<TString>
slist(ilist.size());
6690 for(
int i=0;
i<vec.size();
i++)
slist[
i]=vec[
i];
6703 if(!file.EndsWith(
"."+fext)) {
6704 cout <<
"CWB::Toolbox::getJobId : input file " 6705 << file.Data() <<
" must terminate with ." << fext << endl;
6709 file.ReplaceAll(
"."+fext,
"");
6710 TObjArray*
token = file.Tokenize(
'_');
6711 TObjString* sjobId = (TObjString*)token->At(token->GetEntries()-1);
6712 TString check = sjobId->GetString(); check.ReplaceAll(
"job",
"");
6713 if(!check.IsDigit()) {
6714 cout <<
"CWB::Toolbox::getJobId : input file " << file.Data()
6715 <<
" must terminate with _job#id." << fext << endl;
6716 cout <<
"#id is not a digit" << endl;
6720 if(token)
delete token;
6734 token = param.Tokenize(
TString(
" "));
6735 if((token->GetEntries())>0) {
6736 TObjString* otoken = (TObjString*)token->At(0);
6737 TString stoken = otoken->GetString();
6738 if(((token->GetEntries())>1) || (!stoken.BeginsWith(
"--"))) {
6739 cout <<
"CWB::Toolbox::getParameter : Bad param format \"" << param.Data() <<
"\"" << endl;
6740 cout <<
"Correct format is : --par1" << endl;
6744 if(token)
delete token;
6746 token = options.Tokenize(
TString(
" "));
6748 if((token->GetEntries())%2==1) {
6749 cout <<
"CWB::Toolbox::getParameter : Bad options format \"" << options.Data() <<
"\"" << endl;
6750 cout <<
"Correct format is : --par1 val1 --par2 val2 ..." << endl;
6753 for(
int i=0;
i<token->GetEntries();
i+=2) {
6754 TObjString* otoken = (TObjString*)token->At(
i);
6755 TString stoken = otoken->GetString();
6756 if(!stoken.BeginsWith(
"--")) {
6757 cout <<
"CWB::Toolbox::getParameter : Bad options format \"" << stoken.Data() <<
"\"" << endl;
6758 cout <<
"Correct format is : --par1 val1 --par2 val2 ..." << endl;
6762 for(
int i=0;
i<token->GetEntries();
i++) {
6763 TObjString* otoken = (TObjString*)token->At(
i);
6764 TString stoken = otoken->GetString();
6767 if(i<token->GetEntries()-1) {
6768 otoken = (TObjString*)token->At(
i+1);
6769 return otoken->GetString();
6773 if(token)
delete token;
6785 int MAXSIZE = 0xFFF;
6786 char proclnk[0xFFF];
6787 char filename[0xFFF];
6794 sprintf(proclnk,
"/proc/self/fd/%d", fno);
6795 r = readlink(proclnk, filename, MAXSIZE);
6796 if (r < 0)
return "";
6810 int MAXSIZE = 0xFFF;
6811 char proclnk[0xFFF];
6812 char filename[0xFFF];
6815 if (symlink != NULL)
6817 sprintf(proclnk,
"%s", symlink);
6818 r = readlink(symlink, filename, MAXSIZE);
6819 if (r < 0)
return "";
6833 vector<std::string> ifileList;
6834 vector<std::string> ipathList;
6835 vector<std::string> ofileList;
6836 vector<std::string> opathList;
6839 in.open(ifile.Data());
6841 cout <<
"CWB::Toolbox::getUniqueFileList : Error Opening Input File : " << ifile.Data() << endl;
6848 in.getline(istring,1024);
6849 if (!in.good())
break;
6852 TObjString* stoken =(TObjString*)token->At(token->GetEntries()-1);
6855 ipathList.push_back(istring);
6856 ifileList.push_back(fName.Data());
6861 for(
int i=0;
i<(
int)ifileList.size();
i++) {
6863 for(
int j=0;
j<(
int)ofileList.size();
j++) {
6867 ofileList.push_back(ifileList[
i]);
6868 opathList.push_back(ipathList[i]);
6878 cout <<
"CWB::Toolbox::getUniqueFileList : Error Opening Output File : " << ofile.Data() << endl;
6882 for(
int i=0;
i<(
int)ofileList.size();
i++) {
6883 out << opathList[
i] << endl;
6904 double temp=y[2*
i+1];
6927 y[
i]=sqrt(y[
i]*y[
i]+h[
i]*h[
i]);
6942 double twopi = TMath::TwoPi();
6951 p[
i]=TMath::ATan2(y[
i],h[i]);
6955 for(
int i=1;
i<(
int)p.size()/2-1;
i++) p[2*
i]=(p[2*
i-1]+p[2*
i+1])/2;
6959 y[
i] = (p[
i]-p[
i-1])/(twopi*dt);
6960 if(y[
i]>x.
rate()/2) y[
i]=0;
6963 if(y.
size()>1) y[0]=y[1];
else y[0]=0;
6986 z=0;
for(
int i=dN;
i<2*dN; ++
i) z[
i] = y[
i-dN];
6993 for(
int n=1;
n<=
N; ++
n ) {
6995 for(
int i=0;
i<
N; ++
i) y[
i] = z[dN+2*
n+
i]*z[dN+2*
n-
i];
6996 for(
int i=-N;
i<0; ++
i) y[dN+
i] = z[dN+2*
n+
i]*z[dN+2*
n-
i];
7000 for(
int m=0;
m<
N;
m++) wv[
m*N+(
n-1)] = y[2*
m];
7022 float cutoff = M_PI;
7027 for (j = 0; j < N-1; j++) dp[j] = p[j+1] - p[j];
7031 for (j = 0; j < N-1; j++)
7032 dps[j] = (dp[j]+M_PI) - floor((dp[j]+M_PI) / (2*M_PI))*(2*M_PI) - M_PI;
7036 for (j = 0; j < N-1; j++)
7037 if ((dps[j] == -M_PI) && (dp[
j] > 0))
7042 for (j = 0; j < N-1; j++)
7043 dp_corr[j] = dps[j] - dp[j];
7047 for (j = 0; j < N-1; j++)
7048 if (
fabs(dp[j]) < cutoff)
7053 cumsum[0] = dp_corr[0];
7054 for (j = 1; j < N-1; j++)
7055 cumsum[j] = cumsum[j-1] + dp_corr[j];
7059 for (j = 1; j <
N; j++) p[j] += cumsum[j-1];
7077 if(((a<b && b>=c && b>0)||(a>b && b<=c && b<0)) && (
fabs((a+c)/(2*b))<1)){
7079 omega = acos((a+c)/(2*b))/
dt;
7080 phase = atan2(a*a+c*a-2*b*b, 2*a*b*sin(omega*dt));
7081 amplitude = (a!=0)?(a/cos(phase)):(b/sin(omega*dt));
7108 cout <<
"CWB::Toolbox::WriteFrameFile - Error : wavearray size=0 or rate=0 " << endl;
7116 if(fmod(gps,1)!=0) {
7117 cout <<
"CWB::Toolbox::WriteFrameFile - Error : start gps time is not integer" << endl;
7121 if(fmod(length,1)!=0) {
7122 cout <<
"CWB::Toolbox::WriteFrameFile - Error : length is not integer" << endl;
7126 if(frDir==
"") frDir=
".";
7129 sprintf(frFile,
"%s/%s-%lu-%lu.gwf",frDir.Data(),frLabel.Data(),(
int)gps,(
int)
length);
7130 cout << frFile << endl;
7131 FrFile *ofp = FrFileONew(frFile,1);
7135 FrameH* simFrame = FrameNew(const_cast<char*>(frLabel.Data()));
7136 simFrame->frame = 0;
7139 simFrame->GTimeS =
gps;
7140 simFrame->GTimeN = 0;
7142 cout <<
"Size (sec) " << x.
size()/x.
rate() << endl;
7143 FrProcData* proc = FrProcDataNew(simFrame,const_cast<char*>(chName.Data()),x.
rate(),x.
size(),-64);
7144 if(proc == NULL) {cout <<
"CWB::Toolbox::WriteFrameFile - Cannot create FrProcData" << endl;
exit(-1);}
7145 proc->timeOffset = 0;
7146 proc->tRange = simFrame->dt;
7149 for (
int i=0;
i<(
int)proc->data->nData;
i++) proc->data->dataD[
i] = x[
i];
7151 int err=FrameWrite(simFrame,ofp);
7152 if (err) {cout <<
"CWB::Toolbox::WriteFrameFile - Error writing frame" << endl;
exit(1);}
7153 FrameFree(simFrame);
7155 if (ofp!=NULL) FrFileOEnd(ofp);
7197 for(
int i=0;
i<hi.
size();
i++) {
7202 int Lt =
int(R*(
int(T)+1)+0.01);
7209 for(
int m=HI.
size()-1;
m>0;
m--) {
7225 j=0; dt=1; am=0; l=1; nl=0;
7226 for(
int m=HI.
size()-1;
m>0;
m--) {
7228 if(t<=0 || t>T)
continue;
7229 if(j==0 && am==0.) am=wi.
data[
m];
7230 if(wi.
data[
m]*am <= 0) {
7244 if(a<0) p.
data[
j]+=a*nr;
7258 for(
int j=0;j<s.
size();j++) s[j]=sqrt(s[j]);
7260 for(
int m=p.
size()-2;
m>0;
m--) {
7261 if(p.
data[
m]-p.
data[
m+1] > 1) mm=2;
else mm=0;
7265 for(
int j=0;j<p.
size();j++) {
7266 if(p[j]<-1) p[
j]+=2;
7284 case 0:
return TColor::GetColor( 0, 107, 164);
7285 case 1:
return TColor::GetColor(255, 128, 14);
7286 case 2:
return TColor::GetColor(171, 171, 171);
7287 case 3:
return TColor::GetColor( 89, 89, 89);
7288 case 4:
return TColor::GetColor( 95, 158, 209);
7289 case 5:
return TColor::GetColor(200, 82, 0);
7290 case 6:
return TColor::GetColor(137, 137, 137);
7291 case 7:
return TColor::GetColor(163, 200, 236);
7292 case 8:
return TColor::GetColor(255, 188, 121);
7293 case 9:
return TColor::GetColor(207, 207, 207);
7294 default:
return kBlack;
7311 if ( name ==
"DeepSkyBlue4" )
return getTableau10BlindColor(0);
7312 else if ( name ==
"DarkOrange1" )
return getTableau10BlindColor(1);
7313 else if ( name ==
"DarkGray" )
return getTableau10BlindColor(2);
7314 else if ( name ==
"DimGray" )
return getTableau10BlindColor(3);
7315 else if ( name ==
"SkyBlue3" )
return getTableau10BlindColor(4);
7316 else if ( name ==
"Chocolate3" )
return getTableau10BlindColor(5);
7317 else if ( name ==
"Gray" )
return getTableau10BlindColor(6);
7318 else if ( name ==
"SlateGray1" )
return getTableau10BlindColor(7);
7319 else if ( name ==
"SandyBrown" )
return getTableau10BlindColor(8);
7320 else if ( name ==
"LightGray" )
return getTableau10BlindColor(9);
7336 const int nColors = 4;
7337 TColor * col[nColors] = {
7338 gROOT->GetColor(getTableau10BlindColor(
"SlateGray1")),
7339 gROOT->GetColor(getTableau10BlindColor(
"DeepSkyBlue4")),
7340 gROOT->GetColor(getTableau10BlindColor(
"SandyBrown")),
7341 gROOT->GetColor(getTableau10BlindColor(
"Crimson"))
7343 Double_t
stop[nColors] = {0., .25, .65, 1.0};
7346 Double_t
r[nColors],
g[nColors], b[nColors];
7347 for (
int c = 0;
c < nColors;
c++ ) {
7348 r[
c] = col[
c]->GetRed();
7349 g[
c] = col[
c]->GetGreen();
7350 b[
c] = col[
c]->GetBlue();
7354 Int_t FI = TColor::CreateGradientColorTable(nColors, stop, r, g, b, nSteps);
7355 for (
int i = 0;
i < nSteps;
i++ ) palette[
i] = FI+
i;
wavearray< double > t(hp.size())
size_t append(const wavearray< DataType_t > &)
Float_t * rho
biased null statistics
static double g(double e)
TB createDagFile(rslagList, full_condor_dir, data_label, jobFiles, cwb_stage_name)
virtual void rate(double r)
wavearray< double > a(hp.size())
wavearray< double > Trun(500000)
cout<< "slagList size : "<< slagList.size()<< endl;cout<< endl<< "Start segments selection from dq cat1 list ..."<< endl<< endl;rslagList=TB.getSlagList(slagList, ifos, segLen, segMLS, segEdge, nDQF, DQF, CWB_CAT1);cout<< "Number of selected jobs after cat1 : "<< rslagList.size()<< endl;cout<< endl<< "Start segments selection from dq cat2 list ..."<< endl<< endl;rslagList=TB.getSlagList(rslagList, ifos, segLen, segTHR, segEdge, nDQF, DQF, CWB_CAT2);cout<< "Number of selected jobs after cat2 : "<< rslagList.size()<< endl;vector< TString > jobFiles
cout<< cwb_condor_flen<< endl;double xtime;char jobListFile[256];sprintf(jobListFile,"%s/%s.sjob", dump_dir, data_label);char dq1ListFile[256];sprintf(dq1ListFile,"%s/%s.cat1", dump_dir, data_label);cout<< endl<<"-------------------------------------------------------------------------------------"<< endl<< endl;CWB_CAT dqcat=CWB_CAT1;vector< waveSegment > dq1List
bool TypeAllowed(char *Name)
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 > psd(33)
CWB::frame fr(FRLIST_NAME)
std::vector< waveSegment > olist
TTree * fChain
root input file cointainig the analyzed TTree
cout<< "baudline_FFL : "<< baudline_FFL<< endl;ofstream out;out.open(baudline_FFL, ios::out);if(!out.good()) {cout<< "Error Opening File : "<< baudline_FFL<< endl;exit(1);} ifstream in;in.open(frFiles[ifoID], ios::in);if(!in.good()) {cout<< "Error Opening File : "<< frFiles[ifoID]<< endl;exit(1);} TString pfile_path="";char istring[1024];while(1) { in > istring
virtual void start(double s)
vector< waveSegment > detSegs_dq2
wavearray< int > Wsel(ntrg)
TB mergeTrees(fileList, TREE_NAME, OUTPUT_MERGE_DIR, OFILE_NAME, false)
std::vector< waveSegment > mlist
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
cout<< "ctime : "<< int(ctime)<< " sec "<< ctime/3600.<< " h "<< ctime/86400.<< " day"<< endl;std::vector< waveSegment > jlist
fprintf(stdout,"start=%f duration=%f rate=%f\, x.start(), x.size()/x.rate(), x.rate())
virtual size_t size() const
char * GetLog(char *Stage, int index)
int GetLogSize(char *Stage)
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
printf("total live time: non-zero lags = %10.1f \, liveTot)
TIter next(twave->GetListOfBranches())
Float_t * lag
time between consecutive events
wavearray< double > Tdlag
TString sel("slag[1]:slag[2]")
bool StageAlreadyPresent(char *Name)
vector< TString > getFileListFromDir(TString dir_name, TString endString="", TString beginString="")
vector< waveSegment > cat1List
bool SetHistoryModify(bool Replace=true)
vector< waveSegment > detSegs
vector< waveSegment > iseg
Double_t * time
beam pattern coefficients for hx
condor_log_dir ReplaceAll("X_HOME", uhome.Data())
cout<< "Starting reading output directory ..."<< endl;vector< TString > fileList
virtual void stop(double s)
double fabs(const Complex &x)
void resample(const wavearray< DataType_t > &, double, int=6)
void AddLog(char *Stage, char *Log, TDatime *Time=NULL)
cout<< "Fixed live file : "<< fix_liv_file_name.Data()<< endl;TFile *flive_fix=new TFile(fix_liv_file_name,"RECREATE");if(flive_fix->IsZombie()) { cout<< "CWB::Toolbox::setVeto - Error opening file "<< fix_liv_file_name.Data()<< endl;exit(1);} TTree *tlive_fix=(TTree *) tlive-> CloneTree(0)
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
strcpy(RunLabel, RUN_LABEL)
cout<< "total cat1 livetime : "<< int(cat1_time)<< " sec "<< cat1_time/3600.<< " h "<< cat1_time/86400.<< " day"<< endl;cout<< endl;vector< waveSegment > cat2List
virtual DataType_t max() const
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
char * AddType(char *TypeName)
wavearray< double > Wlag[NIFO_MAX+1]
char * GetHistory(char *StageName, char *Type)
const CWB::HistoryStage * GetStage(char *Name)
for(int i=0;i< 101;++i) Cos2[2][i]=0
pointers to detectors
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
TString frLabel[NIFO_MAX]
void AddHistory(char *Stage, char *Type, char *History, TDatime *Time=NULL)
TB checkFile(gSystem->Getenv("CWB_ROOTLOGON_FILE"))
char * AddStage(char *StageName)
TB dumpSegList(dq1List, dq1ListFile, false)
wavearray< double > Wslag[NIFO_MAX+1]