28 #include <lal/Interpolate.h> 31 #define FAILMSG( stat, func, file, line, id ) \ 33 if ( lalDebugLevel & LALERROR ) \ 35 LALPrintError( "Error[0]: file %s, line %d, %s\n" \ 36 "\tLAL_CALL: Function call `%s' failed.\n", file, line, id, func ); \ 40 #define LAL_CALL( function, statusptr ) \ 41 ((function),lal_errhandler(statusptr,#function,__FILE__,__LINE__,"$Id$")) 43 LALStatus blank_status;
45 typedef int ( *lal_errhandler_t )(
50 volatile const char *
id 58 volatile const char *
id 61 if ( stat->statusCode )
63 FAILMSG( stat, func, file, line,
id );
74 volatile const char *
id 77 if ( stat->statusCode )
79 FAILMSG( stat, func, file, line,
id );
82 return stat->statusCode;
86 lal_errhandler_t lal_errhandler = LAL_ERR_ABRT;
90 ProcessParamsTable **add_process_param(ProcessParamsTable **proc_param,
const ProcessTable *process,
91 const char *type,
const char *param,
const char *
value) {
92 *proc_param = XLALCreateProcessParamsTableRow(process);
93 snprintf((*proc_param)->program,
sizeof((*proc_param)->program),
"%s",
"cWB");
94 snprintf((*proc_param)->type,
sizeof((*proc_param)->type),
"%s", type);
95 snprintf((*proc_param)->param,
sizeof((*proc_param)->param),
"--%s", param);
96 snprintf((*proc_param)->value,
sizeof((*proc_param)->value),
"%s", value);
97 return &(*proc_param)->next;
100 #define ADD_PROCESS_PARAM(proc_param, process, type, param, value) \ 101 do { proc_param = add_process_param(proc_param, process, type, param, value); } while(0) 241 if(pD[
n]->isBuiltin()) _pD =
new detector(pD[
n]->Name);
242 else _pD =
new detector(pD[
n]->getDetectorParams());
310 if(value.
net!=NULL) {
325 if(value.
inj!=NULL) {
409 gRandom->SetSeed(seed);
443 for(
int i=0;
i<(
int)par.size();
i++) {
444 if(par[
i].name.CompareTo(name)==0)
return par[
i].
value;
463 for(
int i=0;
i<(
int)par.size();
i++) {
464 if(par[
i].name.CompareTo(name)==0)
return par[
i].svalue;
581 if(mdc_type<0 || mdc_type>=
MDC_USER) {
582 cout <<
"CWB::mdc::AddWaveform - Error : mdc type " << mdc_type <<
" not allowed !!!" << endl;
588 mdcid waveid = {
"",0,0};
592 int decimals = (
int)
GetPar(
"decimals",par,error);
593 if(error) decimals = -1;
603 if(mdc_type==
MDC_RDE) iota = 0.;
604 if(mdc_type==
MDC_RDC) iota = 0.;
605 if(mdc_type==
MDC_RD) iota = 90.;
608 cout <<
"CWB::mdc::AddWaveform - " 609 <<
"Error : num par must be at least 2 [frequency,tau]" << endl;
613 int d1 = decimals==-1 ? 0 : decimals;
614 int d2 = decimals==-1 ? 1 : decimals;
616 if(mdc_type==
MDC_RD)
sprintf(wf_name,
"RD_%.*f_%.*f",d1,frequency,d2,tau);
617 if(mdc_type==
MDC_RDC)
sprintf(wf_name,
"RDC_%.*f_%.*f",d1,frequency,d2,tau);
618 if(mdc_type==
MDC_RDE)
sprintf(wf_name,
"RDE_%.*f_%.*f",d1,frequency,d2,tau);
622 wf.
name.ReplaceAll(
".",
"d");
623 if(uname!=
"") wf.
name = uname;
624 wf.
hp =
GetRD(frequency,tau,iota,0);
625 wf.
hx =
GetRD(frequency,tau,iota,1);
638 double Q =
GetPar(
"Q" ,par,error);
641 cout <<
"CWB::mdc::AddWaveform - Error : num par must be 2 [frequency,Q]" << endl;
645 int d1 = decimals==-1 ? 0 : decimals;
646 int d2 = decimals==-1 ? 1 : decimals;
648 if(mdc_type==
MDC_SG)
sprintf(wf_name,
"SG%.*fQ%.*f",d1,frequency,d2,Q);
649 if(mdc_type==
MDC_SGC)
sprintf(wf_name,
"SGC%.*fQ%.*f",d1,frequency,d2,Q);
650 if(mdc_type==
MDC_SGE)
sprintf(wf_name,
"SGE%.*fQ%.*f",d1,frequency,d2,Q);
654 wf.
name.ReplaceAll(
".",
"d");
655 if(decimals==-1) wf.
name.ReplaceAll(
"d0",
"");
656 if(uname!=
"") wf.
name = uname;
678 double Q =
GetPar(
"Q" ,par,error);
681 cout <<
"CWB::mdc::AddWaveform - Error : num par must be 2 [frequency,Q]" << endl;
685 int d1 = decimals==-1 ? 0 : decimals;
686 int d2 = decimals==-1 ? 1 : decimals;
688 if(mdc_type==
MDC_CG)
sprintf(wf_name,
"CG%.*fQ%.*f",d1,frequency,d2,Q);
689 if(mdc_type==
MDC_CGC)
sprintf(wf_name,
"CGC%.*fQ%.*f",d1,frequency,d2,Q);
690 if(mdc_type==
MDC_CGE)
sprintf(wf_name,
"CGE%.*fQ%.*f",d1,frequency,d2,Q);
694 wf.
name.ReplaceAll(
".",
"d");
695 if(decimals==-1) wf.
name.ReplaceAll(
"d0",
"");
696 if(uname!=
"") wf.
name = uname;
716 double bandwidth =
GetPar(
"bandwidth",par,error);
720 cout <<
"CWB::mdc::AddWaveform - " 721 <<
"Error : WNB par must at least " 722 <<
"[frequency,bandwidth,duration]" << endl;
726 error=
false;
int pseed = (
int)
GetPar(
"pseed",par,error);
if(error) pseed=0;
727 error=
false;
int xseed = (
int)
GetPar(
"xseed",par,error);
if(error) xseed=0;
728 error=
false;
bool mode = (bool)
GetPar(
"mode",par,error);
if(error) mode=0;
730 int d1 = decimals==-1 ? 0 : decimals;
731 int d2 = decimals==-1 ? 0 : decimals;
732 int d3 = decimals==-1 ? 3 : decimals;
734 sprintf(wf_name,
"WNB%.*f_%.*f_%.*f",d1,frequency,d2,bandwidth,d3,duration);
738 wf.
name.ReplaceAll(
".",
"d");
739 if(uname!=
"") wf.
name = uname;
740 wf.
hp =
GetWNB(frequency,bandwidth,duration,pseed,mode);
742 wf.
hx =
GetWNB(frequency,bandwidth,duration,xseed,mode);
754 cout <<
"CWB::mdc::AddWaveform - " 755 <<
"Error : num par must at least 1 " 756 <<
"[duration]" << endl;
760 int d1 = decimals==-1 ? 1 : decimals;
762 sprintf(wf_name,
"GA%.*f",d1,1000.*duration);
766 wf.
name.ReplaceAll(
".",
"d");
767 if(uname!=
"") wf.
name = uname;
779 #if LAL_VERSION_MAJOR > 6 || (LAL_VERSION_MAJOR == 6 && \ 780 (LAL_VERSION_MINOR > 15 || (LAL_VERSION_MINOR == 15 && \ 781 LAL_VERSION_MICRO >= 0 ))) // LAL_VERSION >= 6.15.0 787 cout <<
"CWB::mdc::AddWaveform - " 788 <<
"Error : wrong input LAL GA parameters : " 789 <<
"[duration,decimals(opt)]" << endl;
794 SimBurst *sim_burst = XLALCreateSimBurst();
795 strcpy(sim_burst->waveform,
"Gaussian");
797 mdcpar mpar={
"normalization",1};par.push_back(mpar);
799 XLALDestroySimBurst(sim_burst);
802 cout <<
"CWB::mdc::AddWaveform - " 803 <<
"Error : MDC_GA_LAL is enabled only with LAL" << endl;
807 cout <<
"CWB::mdc::AddWaveform - MDC_GA_LAL can not be used with LAL ver < 6.15.0" << endl;
817 double Q =
GetPar(
"Q",par,error);
820 cout <<
"CWB::mdc::AddWaveform - " 821 <<
"Error : wrong input LAL SG parameters : " 822 <<
"[frequency,Q,decimals(opt)]" << endl;
827 SimBurst *sim_burst = XLALCreateSimBurst();
828 strcpy(sim_burst->waveform,
"SineGaussian");
834 sim_burst->pol_ellipse_e=0.;
837 sim_burst->pol_ellipse_angle=0.;
838 mdcpar mpar={
"normalization",1};par.push_back(mpar);
840 XLALDestroySimBurst(sim_burst);
843 cout <<
"CWB::mdc::AddWaveform - " 844 <<
"Error : MDC_SGE_LAL is enabled only with LAL" << endl;
854 double bandwidth =
GetPar(
"bandwidth",par,error);
857 int lseed = (
int)
GetPar(
"lseed" ,par,error);
858 int hseed = (
int)
GetPar(
"hseed" ,par,error);
861 cout <<
"CWB::mdc::AddWaveform - " 862 <<
"Error : wrong input LAL WNB parameters : " 863 <<
"[frequency,bandwidth,duration,lseed,hseed,decimals(opt)]" << endl;
868 SimBurst *sim_burst = XLALCreateSimBurst();
869 strcpy(sim_burst->waveform,
"BTLWNB");
872 sim_burst->bandwidth = bandwidth;
873 sim_burst->egw_over_rsquared = 1;
874 sim_burst->waveform_number = hseed*10000000000+lseed;
875 mdcpar mpar={
"normalization",1};par.push_back(mpar);
877 XLALDestroySimBurst(sim_burst);
880 cout <<
"CWB::mdc::AddWaveform - " 881 <<
"Error : MDC_WNB_LAL is enabled only with LAL" << endl;
888 #if LAL_VERSION_MAJOR > 6 || (LAL_VERSION_MAJOR == 6 && \ 889 (LAL_VERSION_MINOR > 15 || (LAL_VERSION_MINOR == 15 && \ 890 LAL_VERSION_MICRO >= 0 ))) // LAL_VERSION >= 6.15.0 897 cout <<
"CWB::mdc::AddWaveform - " 898 <<
"Error : wrong input LAL SC parameters : " 899 <<
"[frequency,amplitude,decimals(opt)]" << endl;
904 SimBurst *sim_burst = XLALCreateSimBurst();
905 strcpy(sim_burst->waveform,
"StringCusp");
910 sim_burst->amplitude = 1;
911 mdcpar mpar={
"normalization",1};par.push_back(mpar);
914 XLALDestroySimBurst(sim_burst);
917 cout <<
"CWB::mdc::AddWaveform - " 918 <<
"Error : MDC_SC_LAL is enabled only with LAL" << endl;
922 cout <<
"CWB::mdc::AddWaveform - MDC_SC_LAL can not be used with LAL ver < 6.15.0" << endl;
932 cout <<
"CWB::mdc::AddWaveform - " 933 <<
"Error : num par must at least 1 " 934 <<
"[file name (.lst/.root) : list of eBBH mdc]" << endl;
940 if(fName.EndsWith(
".lst")) {
945 cout <<
"CWB::mdc::AddWaveform - Error Opening File : " << fName << endl;
953 in.getline(str,1024);
954 if (!in.good())
break;
955 if(str[0] !=
'#') entries++;
957 cout <<
"entries " << entries << endl;
958 in.clear(ios::goodbit);
959 in.seekg(0, ios::beg);
966 in.getline(str,1024);
967 if(str[0] ==
'#')
continue;
968 if (!in.good())
break;
972 std::stringstream linestream(str);
973 if(!(linestream >>
id >> m1 >> m2 >> rp0 >> e0 >> dist >> redshift)) {
976 if(!(linestream >>
id >> m1 >> m2 >> rp0 >> e0 >> dist)) {
979 if(!(linestream >>
id >> m1 >> m2 >> rp0 >> e0)) {
980 cout <<
"CWB::mdc::AddWaveform - Wrong Format for File : " << fName << endl;
981 cout <<
"input line : " << endl;
983 cout <<
"must be : " << endl;
984 cout <<
"event# " <<
"id " <<
" m1 " <<
" m2 " <<
" rp0 " <<
" e0 " << endl;
985 cout <<
"or : " << endl;
986 cout <<
"event# " <<
"id " <<
" m1 " <<
" m2 " <<
" rp0 " <<
" e0 " <<
" dist " << endl;
987 cout <<
"or : " << endl;
988 cout <<
"event# " <<
"id " <<
" m1 " <<
" m2 " <<
" rp0 " <<
" e0 " <<
" dist " <<
" redshift " << endl;
1008 if(fName.EndsWith(
".root")) {
1010 TFile* efile =
new TFile(fName);
1012 cout <<
"CWB::mdc::AddWaveform - Error opening root file : " << fName.Data() << endl;
1021 TTree* etree = (TTree *) efile->Get(
"ebbh");
1023 cout <<
"CWB::mdc::AddWaveform - file : " << fName.Data()
1024 <<
" not contains tree ebbh" << endl;
1028 etree->SetBranchAddress(
"id",&
id);
1029 etree->SetBranchAddress(
"m1",&m1);
1030 etree->SetBranchAddress(
"m2",&m2);
1031 etree->SetBranchAddress(
"rp0",&rp0);
1032 etree->SetBranchAddress(
"e0",&e0);
1033 etree->SetBranchAddress(
"hp",&hp);
1034 etree->SetBranchAddress(
"hx",&hx);
1035 int dstatus = (etree->GetBranch(
"dist")==NULL) ? 0 : etree->SetBranchAddress(
"dist",&dist);
1036 int rstatus = (etree->GetBranch(
"redshift")==NULL) ? 0 : etree->SetBranchAddress(
"redshift",&redshift);
1039 if(par.size()==2) ecut = par[1].
name;
1040 etree->Draw(
"Entry$",ecut,
"goff");
1041 double*
entry = etree->GetV1();
1042 int esize = etree->GetSelectedRows();
1043 for(
int i=0;
i<esize;
i++) {
1044 etree->GetEntry(entry[
i]);
1046 if(rstatus)
sprintf(str,
"%d %f %f %f %f %f %f",
id,m1,m2,rp0,e0,dist,redshift);
1047 else if(dstatus)
sprintf(str,
"%d %f %f %f %f %f",
id,m1,m2,rp0,e0,dist);
1048 else sprintf(str,
"%d %f %f %f %f",
id,m1,m2,rp0,e0);
1060 wf.
par[1].value=entry[
i];
1076 std::stringstream linestream(str);
1077 if(!(linestream >>
id >> m1 >> m2 >> rp0 >> e0 >> dist >> redshift)) {
1078 linestream.str(str);
1080 if(!(linestream >>
id >> m1 >> m2 >> rp0 >> e0 >> dist)) {
1081 linestream.str(str);
1083 if(!(linestream >>
id >> m1 >> m2 >> rp0 >> e0)) {
1084 cout <<
"CWB::mdc::AddWaveform - Wrong Input Parameter Format : " << str << endl;
1085 cout <<
"input line : " << endl;
1086 cout << str << endl;
1087 cout <<
"must be : " << endl;
1088 cout <<
"event# " <<
"id " <<
" m1 " <<
" m2 " <<
" rp0 " <<
" e0 " << endl;
1089 cout <<
"or : " << endl;
1090 cout <<
"event# " <<
"id " <<
" m1 " <<
" m2 " <<
" rp0 " <<
" e0 " <<
" dist " << endl;
1091 cout <<
"or : " << endl;
1092 cout <<
"event# " <<
"id " <<
" m1 " <<
" m2 " <<
" rp0 " <<
" e0 " <<
" dist " <<
" redshift " << endl;
1109 cout <<
"CWB::mdc::AddWaveform - Error : MDC_EBBH not enabled !!!" << endl;
1116 cout <<
"CWB::mdc::AddWaveform - Warning : waveform not added !!!" << endl;
1156 if(mdc_type<0 || mdc_type>=
MDC_USER) {
1157 cout <<
"CWB::mdc::AddWaveform - Error : mdc type " << mdc_type <<
" not allowed !!!" << endl;
1160 if(sim_burst==NULL) {
1161 cout << endl <<
"CWB::mdc::AddWaveform - " 1162 <<
"Error in LAL AddWaveform : SimBurst is NULL" << endl;
1171 int decimals = (
int)
GetPar(
"decimals",par,error);
1172 if(error) decimals = -1;
1174 int normalization = (
int)
GetPar(
"normalization",par,error);
1175 if(error) normalization = 0;
1177 REAL8TimeSeries *
hp=NULL;
1178 REAL8TimeSeries *
hx=NULL;
1183 int ret = XLALGenerateSimBurst(&hp, &hx, sim_burst, deltaT);
1184 if( ret==XLAL_FAILURE ) {
1185 cout << endl <<
"CWB::mdc::AddWaveform - " 1186 <<
"Error in LAL AddWaveform : check sim burst parameters" << endl;
1189 if(hp->data->length!=hx->data->length) {
1190 cout <<
"CWB::mdc::AddWaveform - Error : LAL hp,hx size not equal !!!" << endl;
1209 XLALDestroyREAL8TimeSeries(hp);
1210 XLALDestroyREAL8TimeSeries(hx);
1214 int d1 = decimals==-1 ? 0 : decimals;
1215 int d2 = decimals==-1 ? 1 : decimals;
1218 wf.
par[0].name=
"frequency"; wf.
par[0].value=sim_burst->frequency;
1219 wf.
par[1].name=
"Q"; wf.
par[1].value=sim_burst->q;
1220 wf.
par[2].name=
"eccentricity"; wf.
par[2].value=sim_burst->pol_ellipse_e;
1221 wf.
par[3].name=
"phase"; wf.
par[3].value=sim_burst->pol_ellipse_angle;
1223 sprintf(wf_name,
"LAL_SGE%.*fQ%.*f",d1,wf.
par[0].value,d2,wf.
par[1].value);
1227 wf.
name.ReplaceAll(
".",
"d");
1228 if(decimals==-1) wf.
name.ReplaceAll(
"d0",
"");
1229 if(uname!=
"") wf.
name = uname;
1238 int d1 = decimals==-1 ? 0 : decimals;
1239 int d2 = decimals==-1 ? 0 : decimals;
1240 int d3 = decimals==-1 ? 3 : decimals;
1243 wf.
par[0].name=
"frequency"; wf.
par[0].value=sim_burst->frequency;
1244 wf.
par[1].name=
"bandwidth"; wf.
par[1].value=sim_burst->bandwidth;
1245 wf.
par[2].name=
"duration"; wf.
par[2].value=sim_burst->duration;
1247 sprintf(wf_name,
"LAL_WNB%.*f_%.*f_%.*f",d1,wf.
par[0].value,d2,wf.
par[1].value,d3,wf.
par[2].value);
1251 wf.
name.ReplaceAll(
".",
"d");
1252 if(uname!=
"") wf.
name = uname;
1261 int d1 = decimals==-1 ? 1 : decimals;
1264 wf.
par[0].name=
"duration"; wf.
par[0].value=sim_burst->duration;
1266 sprintf(wf_name,
"LAL_GA%.*f",d1,1000.*wf.
par[0].value);
1270 wf.
name.ReplaceAll(
".",
"d");
1271 if(uname!=
"") wf.
name = uname;
1279 int d1 = decimals==-1 ? 1 : decimals;
1282 wf.
par[0].name=
"frequency"; wf.
par[0].value=sim_burst->frequency;
1283 wf.
par[1].name=
"amplitude"; wf.
par[1].value=sim_burst->amplitude;
1285 sprintf(wf_name,
"LAL_SC%.*f",d1,wf.
par[0].value);
1290 wf.
name.ReplaceAll(
".",
"d");
1291 if(uname!=
"") wf.
name = uname;
1298 cout <<
"CWB::mdc::AddWaveform - Warning : LAL waveform not added !!!" << endl;
1299 mdcid waveid = {
"",0,0};
1340 if(hx_fName.Sizeof()>1) {
1349 cout <<
"CWB::mdc::AddWaveform - Error : hp,hx size not equal !!!" << endl;
1353 cout <<
"CWB::mdc::AddWaveform - Error : hp,hx rate not equal !!!" << endl;
1393 double srate, vector<mdcpar> par) {
1417 if(hx_fName.Sizeof()>1) {
1426 cout <<
"CWB::mdc::AddWaveform - Error : hp,hx size not equal !!!" << endl;
1430 cout <<
"CWB::mdc::AddWaveform - Error : hp,hx rate not equal !!!" << endl;
1436 double hrss_par =
GetPar(
"hrss",par,error);
1437 if(error) hrss_par=0.;
1442 hrssp=sqrt(hrssp/wf.
hp.
rate());
1443 hrssc=sqrt(hrssc/wf.
hx.
rate());
1444 double hrss=sqrt(hrssp*hrssp+hrssc*hrssc);
1455 for(
int i=0;
i<(
int)wf.
hp.
size();
i++) {wf.
hp[
i]/=hrss/hrss_par;wf.
hx[
i]/=hrss/hrss_par;}
1456 hrssp /= hrss/hrss_par;
1457 hrssc /= hrss/hrss_par;
1463 vector<mdcpar> wfpar;
1464 for(
int i=0;
i<par.size();
i++) {
1465 if(par[
i].
name==
"hrss") {par[
i].value=
hrss;bhrss=
true;}
1466 wfpar.push_back(par[
i]);
1468 mdcpar upar = {
"hrss",1.,
""};
if(!bhrss) wfpar.push_back(upar);
1469 mdcpar ppar = {
"hrssp",hrssp,
""}; wfpar.push_back(ppar);
1470 mdcpar cpar = {
"hrssc",hrssc,
""}; wfpar.push_back(cpar);
1496 cout <<
"CWB::mdc::AddWaveform - Error : hp,hx size not equal !!!" << endl;
1500 cout <<
"CWB::mdc::AddWaveform - Error : hp,hx rate not equal !!!" << endl;
1506 cout <<
"CWB::mdc::AddWaveform - Error : mdc type not allowed !!!" << endl;
1523 waveid.
id = ID==-1 ? 0 :
wfList[
ID].list.size();
1571 cout <<
"CWB::mdc::GetBurst - Error : Dummy method : network is not initialized " << endl;
1575 cout <<
"CWB::mdc::GetBurst - Error : x.rate() != " <<
MDC_SAMPLE_RATE << endl;
1633 tShift =
GetDelay(ifo,
"",phi,theta);
1640 bool ellipticity=
false;
1648 double ePlus = ellipticity ? (1+cos(iota*deg2rad)*cos(iota*deg2rad))/2 : 1.;
1649 double eCross = ellipticity ? cos(iota*deg2rad) : 1.;
1651 if(!ellipticity)
srcList[
k].iota=90;
1666 iShift = tShift<0 ? iShift : 0;
1668 w[
i+iShift] = ePlus*fPlus*wf.
hp[
i]+eCross*fCross*wf.
hx[
i];
1669 SimHpHp+=wf.
hp[
i]*wf.
hp[
i];
1670 SimHcHc+=wf.
hx[
i]*wf.
hx[
i];
1671 SimHpHc+=wf.
hp[
i]*wf.
hx[
i];
1676 double SrcHrss=sqrt(SimHpHp+SimHcHc);
1677 std::stringstream linestream;
1682 linestream.str(wf.
par[0].name.Data());
1683 if(!(linestream >>
id >> m1 >> m2 >> rp0 >> e0 >> dist)) {
1703 double eccentricity = iota>1
e-10 ? iota : 1
e-10;
1705 srcList[
k].iota = (1.-sqrt(1-eccentricity*eccentricity))/eccentricity;
1717 double E = eccentricity;
1718 double A = 1./sqrt(2-E*E);
1719 double B = A*sqrt(1-E*E);
1720 SimHpHp *= 1./(A*
A)/2.;
1721 SimHcHc *=
fabs(B)>1
e-5 ? 1./(B*
B)/2. : 0.;
1722 SimHpHc *=
fabs(B)>1
e-5 ? 1./(A*
B)/2. : 0.;
1723 if(SimHcHc==0) SimHcHc=SimHpHp;
1733 SimHpHp *= hrss>0 ? pow(hrss/SrcHrss,2) : pow(
inj_hrss/SrcHrss,2);
1734 SimHcHc *= hrss>0 ? pow(hrss/SrcHrss,2) : pow(
inj_hrss/SrcHrss,2);
1735 SimHpHc *= hrss>0 ? pow(hrss/SrcHrss,2) : pow(
inj_hrss/SrcHrss,2);
1736 w *= hrss>0 ? hrss/SrcHrss :
inj_hrss/SrcHrss;
1744 SimHpHp*=100./pow(rho,2);
1745 SimHcHc*=100./pow(rho,2);
1746 SimHpHc*=100./pow(rho,2);
1754 int sT = TMath::Nint(T*
w.
rate());
1763 if (j>=(
int)x.
size())
break;
1764 if (j>=0) x[
j]=
w[
i];
1769 listLog = listLog+
log;
1772 mdcList.push_back(log.Data());
1804 int estat = gSystem->GetPathInfo(fName.Data(),&
id,&fsize,&
flags,&
mt);
1806 cout <<
"CWB::mdc::ReadWaveform - File : " << fName.Data() <<
" Not Exist" << endl;
1812 in.open(fName.Data(),
ios::in);
1813 if (!in.good()) {cout <<
"CWB::mdc::ReadWaveform - Error Opening File : " << fName.Data() << endl;
exit(1);}
1816 char*
str =
new char[fsize+1];
1819 in.getline(str,fsize+1);
1820 if (!in.good())
break;
1821 if(str[0] !=
'#') size++;
1824 if(tok->GetEntries()!=2) {
1825 cout <<
"CWB::mdc::ReadWaveform - Input file with bad format : must be 2 columns " << endl;
1830 in.clear(ios::goodbit);
1831 in.seekg(0, ios::beg);
1841 int fpos = in.tellg();
1842 in.getline(str,1024);
1843 if (!in.good())
break;
1845 in.seekg(fpos, ios::beg);
1847 if (!in.good())
break;
1848 if(t[cnt]<tMin) tMin=t[
cnt];
1849 if(t[cnt]>tMax) tMax=t[
cnt];
1852 in.seekg(fpos+1, ios::beg);
1892 int estat = gSystem->GetPathInfo(fName.Data(),&
id,&fsize,&
flags,&
mt);
1894 cout <<
"CWB::mdc::ReadWaveform - File : " << fName.Data() <<
" Not Exist" << endl;
1900 in.open(fName.Data(),
ios::in);
1901 if (!in.good()) {cout <<
"CWB::mdc::ReadWaveform - Error Opening File : " << fName.Data() << endl;
exit(1);}
1904 char*
str =
new char[fsize+1];
1907 in.getline(str,fsize+1);
1908 if (!in.good())
break;
1909 if(str[0] !=
'#') size++;
1911 if((tok->GetEntries()!=1)&&(size>1)) {
1912 cout <<
"CWB::mdc::ReadWaveform - Input file with bad format : must be 1 column " << endl;
1915 if((tok->GetEntries()>1)&&(size==1)) {
1916 x.
resize(tok->GetEntries());
1926 in.clear(ios::goodbit);
1927 in.seekg(0, ios::beg);
1935 in.getline(str,1024);
1936 if (!in.good())
break;
1961 double& iota,
double&
hrss,
int&
ID,
int&
id) {
1990 double& iota,
double&
hrss,
int&
ID,
int&
id) {
2010 cout <<
"CWB::mdc::GetSourceCoordinates - Error : injections not defined !!!" << endl;
2014 int id = gRandom->Uniform(0,
thList.size());
2030 cout <<
"CWB::mdc::GetSourceCoordinates - Error : distribution not defined !!!" << endl;
2045 ID = (
int)gRandom->Uniform(0,
wfList.size());
2046 id = (
int)gRandom->Uniform(0,1+
wfList[ID].list.size());
2067 cout <<
"CWB::mdc::GetSourceList - Warning : buffer too small (stop-start)<=2*inj_length !!!" << endl;
2078 int iStart=
int(0.5+start/timeStep);
2079 int iStop=
int(stop/timeStep);
2081 if(iStart==0) iStart+=1;
2082 vector<source> src_list;
2087 cout <<
"CWB::mdc::GetSourceList - Error : injection object is NULL" << endl;
2094 inj_tree->Draw(
"Entry$",cut,
"goff");
2095 int entries =
inj_tree->GetSelectedRows();
2102 inj_tree->SetBranchAddress(
"phi",phi);
2103 inj_tree->SetBranchAddress(
"theta",theta);
2104 inj_tree->SetBranchAddress(
"psi",psi);
2105 inj_tree->SetBranchAddress(
"time",time);
2106 inj_tree->SetBranchAddress(
"type",type);
2109 for(
int i=0;
i<entries;
i++) {
2111 src.
theta = theta[0];
2118 int TYPE = (type[0]-1)<
mdcName.size() ? type[0]-1 : 0;
2120 int id = (
int)gRandom->Uniform(0,1+
wfList[ID].list.size());
2122 src_list.push_back(src);
2145 src.
id = (
int)gRandom->Uniform(0,1+
wfList[src.
ID].list.size());
2150 src_list.push_back(src);
2153 for(
int i=iStart;
i<=iStop;
i++) {
2165 std::stringstream linestream(src.
wf.
par[0].name.Data());
2167 if((linestream >>
id >> m1 >> m2 >> rp0 >> e0 >> dist)) src.
rho=dist;
2170 src_list.push_back(src);
2192 cout <<
"CWB::mdc::GetBurstLog - Error : Dummy method : network is not initialized " << endl;
2199 char logString[10000]=
"";
2204 sprintf(GravEn_SimID,
"%s",hpPath.Data());
2208 double SimHrss = sqrt(SimHpHp+SimHcHc);
2209 double SimEgwR2 = 0.0;
2210 double GravEn_Ampl = src.
rho>0 ? 10.*hrss/src.
rho :
hrss;
2211 double Internal_x = cos(src.
iota*deg2rad);
2212 double Internal_phi = 0.0;
2213 double External_x = cos(src.
theta*deg2rad);
2217 double EarthCtrGPS = src.
gps;
2220 sprintf(logString,
"%s",GravEn_SimID);
2221 sprintf(logString,
"%s %e",logString,SimHrss);
2222 sprintf(logString,
"%s %e",logString,SimEgwR2);
2223 sprintf(logString,
"%s %e",logString,GravEn_Ampl);
2224 sprintf(logString,
"%s %e",logString,Internal_x);
2225 sprintf(logString,
"%s %e",logString,Internal_phi);
2226 sprintf(logString,
"%s %e",logString,External_x);
2227 sprintf(logString,
"%s %e",logString,External_phi);
2228 sprintf(logString,
"%s %e",logString,External_psi);
2230 sprintf(logString,
"%s %10.6f",logString,FrameGPS);
2231 sprintf(logString,
"%s %10.6f",logString,EarthCtrGPS);
2232 sprintf(logString,
"%s %s",logString,SimName);
2233 sprintf(logString,
"%s %e",logString,SimHpHp);
2234 sprintf(logString,
"%s %e",logString,SimHcHc);
2235 sprintf(logString,
"%s %e",logString,SimHpHc);
2241 double IFOctrGPS = EarthCtrGPS;
2252 if(src.
wf.
name==
"eBBH") {
2254 std::stringstream linestream(src.
wf.
par[0].name.Data());
2255 linestream >>
id >> m1 >> m2 >> rp0 >> e0 >>
distance;
2257 double cosiota = cos(src.
iota*deg2rad);
2258 double eff_dist = distance / sqrt(pow(IFOfPlus*(1.+pow(cosiota,2)),2)/4.+pow(IFOfCross*cosiota,2));
2259 sprintf(logString,
"%s %s %10.6f %e %e %g",logString,ifo.Data(),IFOctrGPS,IFOfPlus,IFOfCross,eff_dist/1000.);
2261 sprintf(logString,
"%s %s %10.6f %e %e",logString,ifo.Data(),IFOctrGPS,IFOfPlus,IFOfCross);
2265 if(src.
wf.
name==
"eBBH") {
2267 std::stringstream linestream(src.
wf.
par[0].name.Data());
2268 linestream >>
id >> m1 >> m2 >> rp0 >> e0 >> distance >> redshift;
2269 sprintf(logString,
"%s ebbh ",logString);
2270 sprintf(logString,
"%s mass1 %g ",logString, m1);
2271 sprintf(logString,
"%s mass2 %g ",logString, m2);
2272 sprintf(logString,
"%s rp0 %g ",logString, rp0);
2273 sprintf(logString,
"%s e0 %g ",logString, e0);
2274 sprintf(logString,
"%s distance %g ",logString, distance/1000.);
2275 sprintf(logString,
"%s redshift %g ",logString, redshift);
2278 double mu = (m1*
m2)/(m1+m2);
2280 double mchirp = pow(mu,3./5)*pow(M,2./5);
2281 sprintf(logString,
"%s mchirp %g ",logString, mchirp);
2284 if(src.
rho>0)
sprintf(logString,
"%s distance %g ",logString, src.
rho/1000.);
2306 polarization.ToUpper();
2309 if(polarization.Contains(
"HP")) plot=
Draw(wf.
hp,type,options,color);
2310 if(polarization.Contains(
"HX")) plot=
Draw(wf.
hx,type,options,color);
2330 polarization.ToUpper();
2333 if(polarization.Contains(
"HP")) plot=
Draw(wf.
hp,type,options,color);
2334 if(polarization.Contains(
"HX")) plot=
Draw(wf.
hx,type,options,color);
2363 cout <<
"CWB::mdc::Draw : Error - No events in the selected period" << endl;
2369 int iStart = (tOffset-tWindow/2)*y.
rate();
if(iStart<0) iStart=0;
2370 int iEnd = (tOffset+tWindow/2)*y.
rate();
if(iEnd>y.
size()) iEnd=y.
size();
2375 for(
int i=iStart;
i<iEnd;
i++)
x[
i-iStart]=y[
i];
2421 cout <<
"CWB::mdc::DrawTime : Error - rate must be >0" << endl;
2426 options.ReplaceAll(
" ",
"");
2428 if(options.Contains(
"SAME")&&(
pts!=NULL)) {
2430 if(
pts!=NULL)
delete pts;
2431 char name[32];
sprintf(name,
"TIME-gID:%d",
int(gRandom->Rndm(13)*1.e9));
2432 pts =
new watplot(const_cast<char*>(name),200,20,800,500);
2434 double tStart,tStop;
2435 if(options.Contains(
"ZOOM")) {
2436 options.ReplaceAll(
"ZOOM",
"");
2444 if(!options.Contains(
"SAME")) {
2445 if(!options.Contains(
"A")) options=options+
" A";
2446 if(!options.Contains(
"L")) options=options+
" L";
2447 if(!options.Contains(
"P")) options=options+
" P";
2449 pts->
plot(x, const_cast<char*>(options.Data()), color, tStart, tStop);
2468 options.ReplaceAll(
" ",
"");
2470 if(options.Contains(
"SAME")&&(
pts!=NULL)) {
2472 if(
pts!=NULL)
delete pts;
2473 char name[32];
sprintf(name,
"FREQ-gID:%d",
int(gRandom->Rndm(13)*1.e9));
2474 pts =
new watplot(const_cast<char*>(name),200,20,800,500);
2479 double tStart,tStop;
2480 if(options.Contains(
"ZOOM")) {
2481 options.ReplaceAll(
"ZOOM",
"");
2490 if(options.Contains(
"NOLOGX")) {logx=
false;options.ReplaceAll(
"NOLOGX",
"");}
2492 if(options.Contains(
"NOLOGY")) {logy=
false;options.ReplaceAll(
"NOLOGY",
"");}
2493 pts->
plot(x, const_cast<char*>(options.Data()), color, tStart, tStop,
true, fLow, fHigh);
2514 double fparm=nfact*6;
2516 double tStart,tStop;
2517 if(options.Contains(
"ZOOM")) {
2518 options.ReplaceAll(
"ZOOM",
"");
2528 if(
pts!=NULL)
delete pts;
2529 stft =
new CWB::STFT(x,nfft,noverlap,
"amplitude",
"gauss",fparm);
2533 stft->
Draw(tStart,tStop,fLow,fHigh,0,0,1);
2555 if(ID<0||ID>=(
int)
wfList.size()) {
2556 cout <<
"CWB::mdc::GetWaveform - Error : ID " << ID <<
" not in the list" << endl;
2562 if((
int)
wfList[
ID].list.size()==0) {
2587 std::stringstream linestream(wf.
par[0].name.Data());
2588 if(!(linestream >>
id >> m1 >> m2 >> rp0 >> e0 >> dist >> redshift)) {
2589 linestream.str(wf.
par[0].name.Data());
2591 if(!(linestream >>
id >> m1 >> m2 >> rp0 >> e0 >> dist)) {
2592 linestream.str(wf.
par[0].name.Data());
2594 if(!(linestream >>
id >> m1 >> m2 >> rp0 >> e0)) {
2595 cout <<
"CWB::mdc::GetWaveform - Wrong eBBH parameter format : " << endl;
2596 cout <<
"input line : " << endl;
2597 cout << wf.
par[0].name.Data() << endl;
2598 cout <<
"must be : " << endl;
2599 cout <<
"event# " <<
"id " <<
" m1 " <<
" m2 " <<
" rp0 " <<
" e0 " << endl;
2600 cout <<
"or : " << endl;
2601 cout <<
"event# " <<
"id " <<
" m1 " <<
" m2 " <<
" rp0 " <<
" e0 " <<
" dist " << endl;
2602 cout <<
"or : " << endl;
2603 cout <<
"event# " <<
"id " <<
" m1 " <<
" m2 " <<
" rp0 " <<
" e0 " <<
" dist " <<
" redshift " << endl;
2609 if(wf.
par.size()==1) {
2621 wf.
hp*=distance_source_Kpc/10.;
2622 wf.
hx*=distance_source_Kpc/10.;
2625 if(wf.
par.size()==2) {
2630 TFile* efile =
new TFile(fName);
2632 cout <<
"CWB::mdc::GetWaveform - Error opening root file : " << fName.Data() << endl;
2639 TTree* etree = (TTree *) efile->Get(
"ebbh");
2641 cout <<
"CWB::mdc::GetWaveform - file : " << fName.Data()
2642 <<
" not contains tree ebbh" << endl;
2646 etree->SetBranchAddress(
"hp",&hp);
2647 etree->SetBranchAddress(
"hx",&hx);
2649 etree->GetEntry(entry);
2652 cout <<
"CWB::mdc::GetWaveform - Error : hp,hx size not equal !!!" << endl;
2656 cout <<
"CWB::mdc::GetWaveform - Error : hp,hx rate not equal !!!" << endl;
2668 cout <<
"CWB::mdc::GetWaveform - Error : wrong number of parameters not !!!" << endl;
2674 cout <<
"CWB::mdc::GetWaveform - Error : hp,hp are empty !!!" << endl;
2703 if((
int)
wfList[
id].list.size()==0) {
2809 T = E>0 ?
T/E : 0.5*size/
rate;
2847 double dF=(rate/(double)size)/2.;
2848 for(
int j=0;
j<size/2;
j+=2) {
2849 a = x[
j]*x[
j]+x[
j+1]*x[
j+1];
2853 F = E>0 ?
F/E : 0.5*
rate;
2871 if(efraction<0) efraction=0;
2872 if(efraction>1) efraction=1;
2878 for(
int i=0;
i<
N;
i++) {avr+=
i*x[
i]*x[
i]; E+=x[
i]*x[
i];}
2885 double sum = ((M>=0)&&(M<N)) ? x[M]*x[M] : 0.;
2886 for(
int j=1;
j<
N;
j++) {
2887 a = ((M-
j>=0)&&(M-j<N)) ? x[M-j] : 0.;
2888 b = ((M+j>=0)&&(M+j<N)) ? x[M+j] : 0.;
2892 if(sum/E > efraction)
break;
2912 if(tShift==0)
return;
2915 int ibeg=0;
int iend=0;
2917 if(x[
i]!=0 && ibeg==0) ibeg=
i;
2920 int ilen=iend-ibeg+1;
2922 int ishift =
fabs(tShift)*x.
rate();
2923 int isize = 2*ilen+2*ishift;
2924 isize = isize + (isize%4 ? 4 - isize%4 : 0);
2926 w.rate(x.
rate()); w=0;
2928 for(
int i=0;
i<ilen;
i++) {w[
i+ishift+ilen/2]=x[ibeg+
i];x[ibeg+
i]=0;}
2934 double df = w.rate()/w.size();
2936 for (
int ii=0;ii<(
int)w.size()/2;ii++) {
2937 TComplex X(w[2*ii],w[2*ii+1]);
2938 X=X*C.Exp(TComplex(0.,-2*pi*ii*df*tShift));
2945 for(
int i=0;
i<(
int)w.size();
i++) {
2946 int j=ibeg-(ishift+ilen/2)+
i;
2947 if((j>=0)&&(j<(
int)x.
size())) x[
j]=w[
i];
2964 if(pShift==0)
return;
2967 int ibeg=0;
int iend=0;
2969 if(x[
i]!=0 && ibeg==0) ibeg=
i;
2972 int ilen=iend-ibeg+1;
2975 isize = isize + (isize%4 ? 4 - isize%4 : 0);
2977 w.rate(x.
rate()); w=0;
2979 for(
int i=0;
i<ilen;
i++) {w[
i+isize/4]=x[ibeg+
i];x[ibeg+
i]=0;}
2986 for (
int ii=0;ii<(
int)w.size()/2;ii++) {
2987 TComplex X(w[2*ii],w[2*ii+1]);
2988 X=X*C.Exp(TComplex(0.,-pShift));
2995 for(
int i=0;
i<(
int)w.size();
i++) {
2996 int j=ibeg-isize/4+
i;
2997 if((j>=0)&&(j<(
int)x.
size())) x[
j]=w[
i];
3017 x.
rate(MDC_SAMPLE_RATE);
3021 AddSGBurst(x, amplitude, frequency, duration,0);
3040 x.
rate(MDC_SAMPLE_RATE);
3044 AddCGBurst(x, amplitude, frequency, duration,0);
3049 #define GA_FORMULA_WNB "[0]*TMath::Exp(-TMath::Power((x-[2])/[1],2)/2)" 3068 gRandom->SetSeed(seed);
3072 x.
rate(MDC_SAMPLE_RATE);
3073 for (
int i=0;
i<(
int)x.
size();
i++) x[
i]=gRandom->Gaus(0,1);
3083 int bFrequency = frequency/
df;
3084 int bBandwidth = (bandwidth/2.)/
df;
3086 int bin = 2*(bFrequency+bBandwidth);
3089 for (
int i=1;
i<bBandwidth;
i++) {
3091 y[bin+2*
i+1]=x[2*
i+1];
3094 y[bin-2*
i+1]=-x[2*
i+1];
3100 int bLow = frequency/
df;
3101 int bHigh = (frequency+bandwidth)/df;
3102 for (
int i=0;
i<bLow;
i++) {x[2*
i]=0;x[2*
i+1]=0;}
3103 for (
int i=bHigh;
i<(
int)x.
size()/2;
i++) {x[2*
i]=0;x[2*
i+1]=0;}
3108 double function_range = 1.;
3109 TF1* ga_function =
new TF1(
"Gaussian",
GA_FORMULA_WNB,-function_range/2,function_range/2);
3110 ga_function->SetParameter(0,1);
3111 ga_function->SetParameter(1,duration);
3112 ga_function->SetParameter(2,function_range/2);
3113 for (
int i=0;
i<(
int)x.
size();
i++) x[
i]*=ga_function->Eval(dt*(
i+1));
3119 for (
int i=0;
i<(
int)x.
size();
i++) x[
i]*=(1./sqrt(2.))*1./hrss;
3126 #define GA_FORMULA "[0]*TMath::Exp(-TMath::Power((x-[2])/[1],2))" 3140 double function_range = 1.;
3141 TF1* ga_function =
new TF1(
"Gaussian",
GA_FORMULA,-function_range/2,function_range/2);
3142 ga_function->SetParameter(0,1);
3143 ga_function->SetParameter(1,duration);
3144 ga_function->SetParameter(2,function_range/2);
3150 x.
rate(MDC_SAMPLE_RATE);
3151 for (
int i=0;
i<(
int)x.
size();
i++) x[
i]=ga_function->Eval(dt*(
i+1));
3164 #define RD_FORMULA "(((x-1./4./[2]+TMath::Abs(x-1./4./[2]))/2./(x-1./4./[2]))*[0]*TMath::Cos(TMath::TwoPi()*[2]*x)+[1]*TMath::Sin(TMath::TwoPi()*[2]*x))*TMath::Exp(-x/[3])" // Heaviside in cos like Andrea Vicere' 3184 double c1 = (1+c2*
c2)/2.;
3186 if (c1/c2<1
e-10) c1=0;
3189 double function_range = 1.;
3190 TF1* rd_function =
new TF1(
"RingDown",rd_formula,0.,function_range);
3191 rd_function->SetParameter(2,frequency);
3192 rd_function->SetParameter(3,tau);
3194 double time_offset=0.05;
3199 x.
rate(MDC_SAMPLE_RATE);
3200 rd_function->SetParameter(0,c1);
3201 rd_function->SetParameter(1,0);
3203 double time = dt*(
i+1)-time_offset;
3204 if (time>=0) x[
i]=rd_function->Eval(time);
else x[
i]=0;
3209 y.
rate(MDC_SAMPLE_RATE);
3210 rd_function->SetParameter(0,0);
3211 rd_function->SetParameter(1,c2);
3213 double time = dt*(
i+1)-time_offset;
3214 if (time>=0) y[
i]=rd_function->Eval(time);
else y[
i]=0;
3225 return polarization==0 ?
x :
y;
3243 for (
int i=0;
i <
n;
i++){
3244 td.
data[
i] += v*gRandom->Gaus(0.,1.)+u;
3266 for (i=0; i<
n; i++){
3268 for (j=0; j<
m; j++){
3269 x = gRandom->Exp(v);
3294 double r = td.
rate();
3301 if(m > n/2-1) m = n/2-2;
3303 for (
int i=0;
i <
m;
i++){
3308 a *= TMath::Sqrt(r/sum);
3310 td.
data[n/2+
int(delay*r)] += 0;
3311 for (
int i=1;
i <
m;
i++){
3336 double r = td.
rate();
3343 if(m > n/2-1) m = n/2-2;
3345 for (
int i=0;
i <
m;
i++){
3350 a *= TMath::Sqrt(r/sum);
3353 for (
int i=1;
i <
m;
i++){
3376 double r = td.
rate();
3382 if(m > n/2-1) m = n/2-2;
3386 for (
int i=0;
i <
m;
i++){
3389 gn.
data[m+
i] = g*gRandom->Gaus(0.,1.);
3390 gn.
data[m-
i-1] = g*gRandom->Gaus(0.,1.);
3393 a *= TMath::Sqrt(r/sum);
3395 for (
int i=0;
i <
m;
i++){
3410 #define MNGD_NUMERATOR "([0]*((x-[3])*(x-[3])+y*y)+([0]+3*sqrt(z*z+[1]*[1]))*pow([0]+sqrt(z*z+[1]*[1]),2))" 3411 #define MNGD_DENOMINATOR "(pow(((x-[3])*(x-[3])+y*y)+pow([0]+sqrt(z*z+[1]*[1]),2),5./2.)*pow(z*z+[1]*[1],3./2.))" 3413 #define MNGD_B 0.3 // Kpc 3414 #define MNGD_Md1 6.6e10 // Msol 3415 #define MNGD_A1 5.81 // Kpc 3416 #define MNGD_Md2 -2.9e10 // Msol 3417 #define MNGD_A2 17.43 // Kpc 3418 #define MNGD_Md3 3.3e9 // Msol 3419 #define MNGD_A3 34.86 // Kpc 3421 #define MNGD_XMAX 40 3424 #define MNGD_SOLAR_SISTEM_DISTANCE_FROM_GC 7.62 // 7.62 [+/-0.32] Kpc it.wikipedia.org/wiki/Via_Lattea 3563 gRandom->SetSeed(seed);
3587 int entries =
GetPar(
"entries",par,error);
3589 if(error) rho_min=0.;
3591 if(error) rho_max=0.;
3594 cout <<
"CWB::mdc::SetSkyDistribution - " 3595 <<
"Error : entries must be positive" << endl;
3598 if(rho_max>0 && rho_min<=0) {
3599 cout <<
"CWB::mdc::SetSkyDistribution - " <<
"Error : rho_min must be > 0" << endl;
3602 if(rho_min>rho_max) {
3603 cout <<
"CWB::mdc::SetSkyDistribution - " <<
"Error : rho_min must be <= rho_max" << endl;
3607 cout <<
"CWB::mdc::SetSkyDistribution - All Sky random distribution" << endl;
3611 if(rho_min!=rho_max) {
3612 double rho_dist =
GetPar(
"rho_dist",par,error);
3613 char rho_dist_func[256];
3614 if((!error)&&(fName==
""))
sprintf(rho_dist_func,
"pow(x,%f)",rho_dist);
3615 else if(( error)&&(fName==
""))
sprintf(rho_dist_func,
"pow(x,2)");
3616 else if(( error)&&(fName!=
""))
strcpy(rho_dist_func,fName.Data());
3617 else if((!error)&&(fName!=
"")) {
3618 cout <<
"CWB::mdc::SetSkyDistribution - " 3619 <<
"Error : ambiguous rho distribution - is defined twice " 3620 <<
" 1) as rho_dist params, 2) as formula" << endl;
3623 rd =
new TF1(
"rd",rho_dist_func,rho_min,rho_max);
3624 TF1* rdcheck = (TF1*)gROOT->GetListOfFunctions()->FindObject(
"rd");
3626 cout <<
"CWB::mdc::SetSkyDistribution - " 3627 <<
"Error : wrong formula format : " << rho_dist_func << endl;
3632 for(
int n=0;
n<entries;
n++) {
3633 thList.push_back(rad2deg*acos(gRandom->Uniform(-1,1))-90.);
3634 phList.push_back(gRandom->Uniform(-180,180));
3635 if(rd!=NULL)
rhoList.push_back(rd->GetRandom());
3636 else rhoList.push_back(rho_min);
3639 double psi = error ? gRandom->Uniform(0,180) :
value;
3642 error=
false;
double iota =
GetPar(
"iota",par,error);
3644 else if(iota>=0&&iota<=180)
iotaList.push_back(iota);
3645 else iotaList.push_back(rad2deg*acos(gRandom->Uniform(-1,1)));
3652 if(rd!=NULL)
delete rd;
3658 double entries =
GetPar(
"entries",par,error);
3661 cout <<
"CWB::mdc::SetSkyDistribution - " 3662 <<
"Error : num par must be 1 [entries]" << endl;
3666 cout <<
"CWB::mdc::SetSkyDistribution - the Miyamoto-Nagai Galactic Disk Model" << endl;
3673 gd1->SetParameter(1,
MNGD_B);
3679 gd2->SetParameter(1,
MNGD_B);
3685 gd3->SetParameter(1,
MNGD_B);
3699 xyz.SetXYZ(xgc,ygc,
zgc);
3707 double gc_rho = sqrt(xyz.mag2());
3709 cout <<
"gc_phi : " << gc_phi <<
" gc_theta : " << gc_theta <<
" " << gc_rho << endl;
3712 for(
int n=0;
n<entries;
n++) {
3714 gd->GetRandom3(x,y,z);
3717 double ilongitude = xyz.Phi()*
rad2deg;
3721 phList.push_back(olongitude);
3722 thList.push_back(olatitude);
3723 psiList.push_back(gRandom->Uniform(0,180));
3724 rhoList.push_back(sqrt(xyz.mag2()));
3726 error=
false;
double iota =
GetPar(
"iota",par,error);
3728 else if(iota>=0&&iota<=180)
iotaList.push_back(iota);
3729 else iotaList.push_back(rad2deg*acos(gRandom->Uniform(-1,1)));
3745 double distance_thr =
GetPar(
"distance_thr",par,error);
3748 cout <<
"CWB::mdc::SetSkyDistribution - " 3749 <<
"Error : num par must be 1 [distance_thr]" << endl;
3753 cout <<
"CWB::mdc::SetSkyDistribution - Gravitational Wave Galaxy Catalog" << endl;
3754 cout <<
"CWB::mdc::SetSkyDistribution - Distance Threshold " << distance_thr <<
" Kpc" << endl;
3758 if (!in.good()) {cout <<
"CWB::mdc::SetSkyDistribution - Error Opening File : " << fName << endl;
exit(1);}
3764 in.getline(str,1024);
3765 if (!in.good())
break;
3766 if(str[0] !=
'#') entries++;
3768 cout <<
"entries " << entries << endl;
3769 in.clear(ios::goodbit);
3770 in.seekg(0, ios::beg);
3773 in.getline(iline,1024);
3776 in.getline(iline,1024);
3777 if (!in.good())
break;
3780 TObjString*
tra = (TObjString*)tok->At(2);
3781 TObjString*
tdec = (TObjString*)tok->At(3);
3782 TObjString*
tdist = (TObjString*)tok->At(14);
3786 double ra = tra->GetString().Atof();
3787 double dec = tdec->GetString().Atof();
3788 double dist = tdist->GetString().Atof();
3791 if (dist<distance_thr) {
3793 phList.push_back(ra*360./24.);
3794 psiList.push_back(gRandom->Uniform(0,180));
3797 error=
false;
double iota =
GetPar(
"iota",par,error);
3799 else if(iota>=0&&iota<=180)
iotaList.push_back(iota);
3800 else iotaList.push_back(rad2deg*acos(gRandom->Uniform(-1,1)));
3811 cout <<
"CWB::mdc::SetSkyDistribution - User Custom Distribution" << endl;
3815 if (!in.good()) {cout <<
"CWB::mdc::SetSkyDistribution - Error Opening File : " << fName << endl;
exit(1);}
3821 in.getline(str,1024);
3822 if (!in.good())
break;
3823 if(str[0] !=
'#') entries++;
3825 cout <<
"entries " << entries << endl;
3826 in.clear(ios::goodbit);
3827 in.seekg(0, ios::beg);
3836 in.getline(str,1024);
3837 if(str[0] ==
'#')
continue;
3838 if (!in.good())
break;
3840 std::stringstream linestream(str);
3841 if(!(linestream >> gps >> name >> theta >> phi >> psi >> rho >> iota >> hrss >> ID >>
id)) {
3842 cout <<
"CWB::mdc::SetSkyDistribution - Wrong Format for File : " << fName << endl;
3843 cout <<
"input line : " << endl;
3844 cout << str << endl;
3845 cout <<
"must be : " << endl;
3846 cout <<
"gps " <<
"name " <<
"theta " <<
"phi " <<
"psi " 3847 <<
"rho " <<
"iota " <<
"hrss " <<
"ID " <<
"id " <<
"..." << endl;
3871 cout <<
"CWB::mdc::SetSkyDistribution - User Distribution from XML" << endl;
3873 #if LAL_VERSION_MAJOR > 6 || (LAL_VERSION_MAJOR == 6 && \ 3874 (LAL_VERSION_MINOR > 14 || (LAL_VERSION_MINOR == 14 && \ 3875 LAL_VERSION_MICRO >= 0 ))) // LAL_VERSION >= 6.14.0 3877 cout <<
"CWB::mdc::SetSkyDistribution - LAL XML File can not be used with LAL ver < 6.14.0" << endl;
3882 int estat = gSystem->GetPathInfo(fName.Data(),&
id,&fsize,&
flags,&
mt);
3884 cout <<
"CWB::mdc::SetSkyDistribution - XML File : " << fName.Data() <<
" Not Exist" << endl;
3892 int decimals =
GetPar(
"decimals",par,error);
3893 if(error) decimals=1;
3897 double hrss_factor =
GetPar(
"hrss_factor",par,error);
3898 if(error) hrss_factor=1;
3902 int engine =
GetPar(
"engine",par,error);
3907 int AUTO =
GetPar(
"auto",par,error);
3916 if(error) stop=std::numeric_limits<double>::max();;
3920 for(
int n=0;
n<par.size();
n++) {
3922 vector<mdcpar> xpar(1); xpar[0]=par[
n];
3927 if(type==
xmlType[
m]) {present=
true;
break;}
3929 if(!present)
xmlType.push_back(type.Data());
3935 LIGOTimeGPS gpsStartTime = {(
int)start, 0};
3936 LIGOTimeGPS gpsStopTime = {(
int)stop, 0};
3937 SimBurst* sim_burst = XLALSimBurstTableFromLIGOLw(fName.Data(),&gpsStartTime,&gpsStopTime);
3938 if(sim_burst==NULL) {
3939 cout <<
"CWB::mdc::SetSkyDistribution - XML File : " << fName.Data() << endl;
3940 cout <<
" No Events present in the range : " << (
int)start <<
":" << (
int)stop <<endl;
3943 int sim_burst_index=0;
3944 for(; sim_burst; sim_burst = sim_burst->next) {
3950 double gps = sim_burst->time_geocent_gps.gpsSeconds;
3951 gps += 1.e-9*sim_burst->time_geocent_gps.gpsNanoSeconds;
3953 if(gps<start || gps>stop)
continue;
3958 double phi = sm.
RA2phi(rad2deg*sim_burst->ra,gps);
3959 double theta = rad2deg*(TMath::PiOver2()-sim_burst->dec);
3960 double psi = rad2deg*sim_burst->psi;
3961 double frequency = sim_burst->frequency;
3962 double Q = sim_burst->q;
3963 double duration = sim_burst->duration;
3964 double bandwidth = sim_burst->bandwidth;
3966 double phase = rad2deg*sim_burst->pol_ellipse_angle;
3967 double eccentricity = sim_burst->pol_ellipse_e;
3968 double hrss = sim_burst->hrss * hrss_factor;
3969 double amplitude = sim_burst->amplitude;
3970 double egw_over_rsquared = sim_burst->egw_over_rsquared;
3971 double waveform_number = sim_burst->waveform_number;
3982 double iota = eccentricity;
3987 vector<mdcpar> wpar;
3990 if(
TString(sim_burst->waveform)==
"Gaussian") {
3994 wpar[0].name=
"duration"; wpar[0].value=
duration;
3995 wpar[1].name=
"decimals"; wpar[1].value=decimals;
4001 if(
TString(sim_burst->waveform)==
"SineGaussian") {
4005 wpar[0].name=
"frequency"; wpar[0].value=
frequency;
4006 wpar[1].name=
"Q"; wpar[1].value=
Q;
4007 wpar[2].name=
"decimals"; wpar[2].value=decimals;
4013 if(
TString(sim_burst->waveform)==
"BTLWNB") {
4017 wpar[0].name=
"frequency"; wpar[0].value=
frequency;
4018 wpar[1].name=
"bandwidth"; wpar[1].value=bandwidth;
4019 wpar[2].name=
"duration"; wpar[2].value=
duration;
4020 wpar[3].name=
"pseed"; wpar[3].value=seed+2*sim_burst_index;
4021 wpar[4].name=
"xseed"; wpar[4].value=seed+2*sim_burst_index+1;
4022 wpar[5].name=
"mode"; wpar[5].value=0;
4023 wpar[6].name=
"decimals"; wpar[6].value=decimals;
4029 if(
TString(sim_burst->waveform)==
"StringCusp") {
4033 wpar[0].name=
"frequency"; wpar[0].value=
frequency;
4034 wpar[1].name=
"amplitude"; wpar[1].value=
amplitude;
4035 wpar[2].name=
"decimals"; wpar[2].value=decimals;
4059 cout <<
"CWB::mdc::SetSkyDistribution - " 4060 <<
"Error : mdc type " << waveid.
name <<
" is not in the xmlType list" << endl;
4061 cout<<endl<<
"xmlType list : "<<endl<<endl;
4063 cout <<
"xmlType[" <<
n <<
"] : " <<
xmlType[
n] << endl;
4088 cout <<
"CWB::mdc::SetSkyDistribution - " 4089 <<
"Error : User Distribution from XML is enabled only with LAL" << endl;
4100 if(gerror) gerror=
error;
4102 if(gerror) gerror=
error;
4105 cout <<
"CWB::mdc::SetSkyDistribution - " 4106 <<
"Error : num par must be at least 2 [theta,phi,(psi,rho,gps)]" << endl;
4111 value=psi=rho=iota=0.;
4112 error=
false;value=
GetPar(
"entries",par,error);
4113 int entries = error ? 10000 :
int(value);
4114 for(
int n=0;
n<entries;
n++) {
4115 error=
false;value=
GetPar(
"psi",par,error);
4116 psi = error ? gRandom->Uniform(0,180) :
value;
4117 error=
false;value=
GetPar(
"rho",par,error);
4118 rho = error ? 0. :
value;
4119 error=
false;value=
GetPar(
"gps",par,error);
4120 if(!error)
gpsList.push_back(value);
4122 error=
false;iota =
GetPar(
"iota",par,error);
4124 else if(iota>=0&&iota<=180)
iotaList.push_back(iota);
4125 else iotaList.push_back(rad2deg*acos(gRandom->Uniform(-1,1)));
4138 cout <<
"CWB::mdc::SetSkyDistribution - Error : Dummy method : network is not initialized " << endl;
4141 int nInj=
net->
readMDClog(const_cast<char*>(fName.Data()),0.);
4142 printf(
"CWB::mdc::SetSkyDistribution - injections loaded : %d\n",nInj);
4146 cout << N <<
" " << M << endl;
4176 if(
psp!=NULL)
delete psp;
4191 int size=360*2*resolution*180*2*resolution;
4198 cout <<
"CWB::mdc::DrawSkyDistribution - Error : injection object is NULL" << endl;
4202 inj_tree->Draw(
"theta[0]:phi[0]:distance[0]:time[0]",
"",
"goff");
4203 entries =
inj_tree->GetSelectedRows();
4210 coordinate.ToUpper();
4211 if(coordinate.CompareTo(
"CELESTIAL")==0) {
4214 for(
int n=0;
n<entries;
n++) x[
n] = phi[
n];
4216 for(
int n=0;
n<entries;
n++) {
4219 if(zmax<z[
n]) zmax=z[
n];
4226 for(
int n=0;
n<entries;
n++) {
4230 if(zmax<z[
n]) zmax=z[
n];
4235 cout <<
"CWB::mdc::DrawSkyDistribution - Warning : no entries in sky distribution " << endl;
4242 for (
int i=0;
i<360*2*resolution;
i++) {
4243 for (
int j=0;
j<180*2*resolution;
j++) {
4244 double ph =
i/(double)(2*resolution);
4245 double th =
j/(double)(2*resolution);
4257 TMath::Sort(size,z.
data,index,
true);
4259 for (
int i=0;
i<
size;
i++) T[
i]=x[index[
i]];
4261 for (
int i=0;
i<
size;
i++) T[
i]=y[index[
i]];
4263 for (
int i=0;
i<
size;
i++) T[
i]=z[index[
i]];
4302 cout <<
"CWB::mdc::WriteFrameFile - Error : Dummy method : network is not initialized " << endl;
4311 if(chName.size()!=0 && chName.size()!=
nIFO) {
4312 cout <<
"CWB::mdc::WriteFrameFile - Error : chName list size must be equals to nIFO " << nIFO << endl;
4318 sprintf(sdir,
"%s-%s-%04d",ifoLabel,frLabel.Data(),
int(gps/100000));
4319 char cmd[128];
sprintf(cmd,
"mkdir -p %s/%s",frDir.Data(),sdir);
4320 cout << cmd << endl;
4325 cout << frFile << endl;
4326 FrFile *ofp = FrFileONew(frFile,1);
4330 FrameH* simFrame = FrameNew(const_cast<char*>(ifoLabel));
4331 simFrame->frame = 0;
4334 simFrame->GTimeS =
gps;
4335 simFrame->GTimeN = 0;
4345 if(chName[
i]!=
"")
sprintf(chNAME,
"%s",chName[
i].Data());
4353 cout <<
"Size (sec) " << x.
size()/x.
rate() << endl;
4354 FrProcData* proc = FrProcDataNew(simFrame,chNAME,x.
rate(),x.
size(),-64);
4355 if(proc == NULL) {cout <<
"CWB::mdc::WriteFrameFile - Cannot create FrProcData" << endl;
exit(-1);}
4356 proc->timeOffset = 0;
4357 proc->tRange = simFrame->dt;
4360 for (
int i=0;
i<(
int)proc->data->nData;
i++) proc->data->dataD[
i] = x[
i];
4363 int err=FrameWrite(simFrame,ofp);
4364 if (err) {cout <<
"CWB::mdc::WriteFrameFile - Error writing frame" << endl;
exit(1);}
4365 FrameFree(simFrame);
4367 if (ofp!=NULL) FrFileOEnd(ofp);
4377 logFile.ReplaceAll(
".gwf",
"-Log.txt");
4405 polarization.ToUpper();
4408 if(polarization.Contains(
"HP"))
Dump(fname, wf.
hp);
4409 if(polarization.Contains(
"HX"))
Dump(fname, wf.
hx);
4427 polarization.ToUpper();
4430 if(polarization.Contains(
"HP"))
Dump(fname, wf.
hp);
4431 if(polarization.Contains(
"HX"))
Dump(fname, wf.
hx);
4449 if (!out.good()) {cout <<
"CWB::mdc::WriteFrameFile - Error Opening File : " << fname.Data() << endl;
exit(1);}
4451 for(
int i=0;
i<(
int)x.
size();
i++) out << x[
i] << endl;
4470 sprintf(cmd,
"mkdir -p %s",dname.Data());
4513 else pD =
new detector(const_cast<char*>(ifo.Data()));
4516 if(ifoId<0)
delete pD;
4517 polarization.ToUpper();
4518 if(polarization.Contains(
"HP"))
return F.
real();
else return F.
imag();
4535 if(ifo1.CompareTo(ifo2)==0)
return 0.0;
4542 if(ifoId1>=0) pD1 =
net->
getifo(ifoId1);
4543 else pD1 =
new detector(const_cast<char*>(ifo1.Data()));
4546 if(ifoId1<0)
delete pD1;
4547 double delay = pD1->
getTau(theta,phi);
4554 if(ifoId2>=0) pD2 =
net->
getifo(ifoId2);
4555 else pD2 =
new detector(const_cast<char*>(ifo2.Data()));
4557 double delay = pD1->
getTau(theta,phi)-pD2->
getTau(theta,phi);
4559 if(ifoId1<0)
delete pD1;
4560 if(ifoId2<0)
delete pD2;
4577 if(size<=0 &&
inj==NULL) {
4578 cout <<
"CWB::mdc::DumpLogHeader - Error : dump needs injection object " << endl;
4586 if (!out.good()) {cout <<
"CWB::mdc::DumpLogHeader - Error Opening File : " << fName.Data() << endl;
exit(1);}
4590 out <<
"# Log File for Burst " << label.Data() << endl;
4592 out <<
"# The following detectors were simulated" << endl;
4595 out <<
"# There were " << size <<
" injections" << endl;
4597 out <<
"# There were " <<
net->
mdcList.size() <<
" injections" << endl;
4605 out <<
"# Burst MDC Sim type " <<
net->
mdcType[
i].data() <<
"\t occurred " << ninj <<
" times" << endl;
4609 out <<
"# GravEn_SimID SimHrss SimEgwR2 GravEn_Ampl ";
4610 out <<
"Internal_x Internal_phi ";
4611 out <<
"External_x External_phi External_psi ";
4612 out <<
"FrameGPS EarthCtrGPS SimName SimHpHp SimHcHc SimHpHc ";
4629 CWB::mdc::CreateBurstXML(
TString fName, vector<mdcpar> xml_parms) {
4639 cout <<
"CWB::mdc::CreateBurstXML : Error : This function is allowed only for Burst !!!" << endl;
4644 cout <<
"CWB::mdc::CreateBurstXML : Error : file XML already opened !!!" << endl;
4650 double start_time =
GetPar(
"start_time",xml_parms,error);
if(error) start_time = 0;
4652 double end_time =
GetPar(
"end_time",xml_parms,error);
if(error) end_time = 0;
4655 xml_process_table_head = NULL;
4656 xml_process_params_table_head = NULL;
4657 xml_search_summary_table_head = NULL;
4658 xml_time_slide_table_head = NULL;
4659 xml_sim_burst_table_head = NULL;
4660 xml_sim_burst = &xml_sim_burst_table_head;
4664 if(ifos!=
"") ifos.Resize(ifos.Sizeof()-2);
4667 xml_process_table_head = xml_process = XLALCreateProcessTableRow();
4668 sprintf(xml_process->program,
"cWB");
4669 XLALGPSTimeNow(&xml_process->start_time);
4671 sprintf(xml_process->cvs_repository,
"https://git.ligo.org/cWB/library/");
4673 for(
int i=0;
i<gApplication->Argc();
i++)
sprintf(cmd_line,
"%s %s",cmd_line,gApplication->Argv(
i));
4674 sprintf(xml_process->comment, cmd_line);
4675 sprintf(xml_process->domain,
"cWB");
4677 UserGroup_t*
uinfo = gSystem->GetUserInfo();
4678 xml_process->unix_procid = gSystem->GetPid();
4679 sprintf(xml_process->node, gSystem->HostName());
4680 sprintf(xml_process->username, uinfo->fUser);
4681 sprintf(xml_process->ifos,ifos.Data());
4682 xml_process->process_id = 0;
4685 ProcessParamsTable **paramaddpoint = &xml_process_params_table_head;
4687 ADD_PROCESS_PARAM(paramaddpoint, xml_process,
"string",
"instrument",
net->
getifo(
n)->
Name);
4689 ADD_PROCESS_PARAM(paramaddpoint, xml_process,
"double",
"inj-hrss",
4691 ADD_PROCESS_PARAM(paramaddpoint, xml_process,
"double",
"time-step",
4692 TString::Format(
"%f",1./
GetInjRate()).Data());
4693 ADD_PROCESS_PARAM(paramaddpoint, xml_process,
"double",
"jitter",
4695 ADD_PROCESS_PARAM(paramaddpoint, xml_process,
"enum",
"sky-distribution",
4697 ADD_PROCESS_PARAM(paramaddpoint, xml_process,
"string",
"sky-file",
sky_file.Data());
4699 ADD_PROCESS_PARAM(paramaddpoint, xml_process,
"double",
4702 for(
int i=0;
i<xml_parms.size();
i++) {
4703 if(xml_parms[
i].svalue!=
"") {
4704 ADD_PROCESS_PARAM(paramaddpoint, xml_process,
"string",
4705 xml_parms[
i].
name.Data(), xml_parms[
i].svalue.Data());
4707 ADD_PROCESS_PARAM(paramaddpoint, xml_process,
"double",
4708 xml_parms[
i].
name.Data(), TString::Format(
"%f",xml_parms[
i].
value).Data());
4713 xml_search_summary_table_head = xml_search_summary = XLALCreateSearchSummaryTableRow(xml_process);
4714 sprintf(xml_search_summary->comment,
"");
4715 sprintf(xml_search_summary->ifos,ifos.Data());
4716 xml_search_summary->nnodes = 1;
4718 xml_search_summary->in_start_time.gpsSeconds = in_start_time.
GetSec();
4719 xml_search_summary->in_start_time.gpsNanoSeconds = in_start_time.
GetNSec();
4720 xml_search_summary->out_start_time.gpsSeconds = in_start_time.
GetSec();
4721 xml_search_summary->out_start_time.gpsNanoSeconds = in_start_time.
GetNSec();
4723 xml_search_summary->in_end_time.gpsSeconds = in_end_time.
GetSec();
4724 xml_search_summary->in_end_time.gpsNanoSeconds = in_end_time.
GetNSec();
4725 xml_search_summary->out_end_time.gpsSeconds = in_end_time.
GetSec();
4726 xml_search_summary->out_end_time.gpsNanoSeconds = in_end_time.
GetNSec();
4734 CWB::mdc::FillBurstXML(
bool verbose) {
4742 cout <<
"CWB::mdc::FillBurstXML : Error : file XML not declared !!! use CreateBurstXML" << endl;
4746 double NaN = std::numeric_limits<double>::quiet_NaN();
4750 cout <<
"gps" <<
"\t\t" <<
"wf-name" <<
"\t\t\t" <<
"theta" <<
"\t" <<
"phi" <<
"\t" <<
"psi" 4751 <<
"\t" <<
"rho" <<
"\t" <<
"iota" <<
"\t" <<
"hrss" <<
"\t" <<
"ID" <<
"\t" <<
"id" << endl;
4759 cout << std::setprecision(13) <<
srcList[
i].gps <<
"\t" << std::setprecision(6)
4765 SimBurst *sim_burst = XLALCreateSimBurst();
4766 if(!sim_burst) {*xml_sim_burst=NULL;
return;}
4772 sim_burst->time_geocent_gps.gpsSeconds = time.GetSec();
4773 sim_burst->time_geocent_gps.gpsNanoSeconds = time.GetNSec();
4776 sim_burst->dec = TMath::PiOver2()-deg2rad*
srcList[
i].theta;
4777 sim_burst->psi = deg2rad*
srcList[
i].psi;
4779 sim_burst->frequency =
GetPar(
"frequency",wf.
par,error);
if(error) sim_burst->frequency=NaN;
4781 sim_burst->q =
GetPar(
"Q",wf.
par,error);
if(error) sim_burst->q=NaN;
4783 sim_burst->duration =
GetPar(
"duration",wf.
par,error);
if(error) sim_burst->duration=NaN;
4785 sim_burst->bandwidth =
GetPar(
"bandwidth",wf.
par,error);
if(error) sim_burst->bandwidth=NaN;
4786 sim_burst->pol_ellipse_angle = 0;
4787 sim_burst->pol_ellipse_e = cos(deg2rad*
srcList[
i].iota);
4788 sim_burst->hrss =
hrss;
4789 sim_burst->amplitude = 0;
4790 sim_burst->egw_over_rsquared = 0;
4791 sim_burst->waveform_number =
srcList[
i].ID;
4793 *xml_sim_burst = sim_burst;
4794 xml_sim_burst = &(*xml_sim_burst)->next;
4803 CWB::mdc::CloseBurstXML() {
4811 cout <<
"CWB::mdc::CloseBurstXML : Error : file XML not declared !!! use CreateBurstXML" << endl;
4815 XLALGPSTimeNow(&xml_process->end_time);
4818 LIGOLwXMLStream *xml;
4823 if(XLALWriteLIGOLwXMLProcessTable(xml, xml_process_table_head)) {
4824 cout <<
"CWB::mdc::CloseBurstXML : Error in XLALWriteLIGOLwXMLProcessTable" << endl;
4828 if(XLALWriteLIGOLwXMLProcessParamsTable(xml, xml_process_params_table_head)) {
4829 cout <<
"CWB::mdc::CloseBurstXML : Error in XLALWriteLIGOLwXMLProcessParamsTable" << endl;
4833 if(XLALWriteLIGOLwXMLSearchSummaryTable(xml, xml_search_summary_table_head)) {
4834 cout <<
"CWB::mdc::CloseBurstXML : Error in XLALWriteLIGOLwXMLSearchSummaryTable" << endl;
4838 if(XLALWriteLIGOLwXMLTimeSlideTable(xml, xml_time_slide_table_head)) {
4839 cout <<
"CWB::mdc::CloseBurstXML : Error in XLALWriteLIGOLwXMLTimeSlideTable" << endl;
4843 if(XLALWriteLIGOLwXMLSimBurstTable(xml, xml_sim_burst_table_head)) {
4844 cout <<
"CWB::mdc::CloseBurstXML : Error in XLALWriteLIGOLwXMLSimBurstTable" << endl;
4848 XLALCloseLIGOLwXMLFile(xml);
4851 XLALDestroyProcessTable(xml_process_table_head);
4852 XLALDestroyProcessParamsTable(xml_process_params_table_head);
4853 XLALDestroyTimeSlideTable(xml_time_slide_table_head);
4854 XLALDestroySearchSummaryTable(xml_search_summary_table_head);
4855 XLALDestroySimBurstTable(xml_sim_burst_table_head);
4880 if(name==
"") error=
true;
4881 if((name==
"Gaussian")&&(name!=
"SineGaussian")&&(name==
"BTLWNB")) error=
true;
4883 cout <<
"CWB::mdc::GetBurstNameLAL : Error - Not valid --name option" << endl;
4884 cout <<
"available --name options are : Gaussian/SineGaussian/BTLWNB" << endl;
4891 if(sdecimals.IsDigit()) decimals = sdecimals.Atoi();
4893 if(name==
"Gaussian") {
4895 int d1 = decimals==-1 ? 1 : decimals;
4899 if(sduration.IsFloat()) {
4900 duration = sduration.Atof();
4902 cout <<
"CWB::mdc::GetBurstNameLAL : Error : --duration not defined or not a number!!! " << endl;
4906 wfname = TString::Format(
"LAL_GA%.*f",d1,1000.*duration);
4907 wfname.ReplaceAll(
".",
"d");
4910 }
else if(name==
"SineGaussian") {
4912 int d1 = decimals==-1 ? 0 : decimals;
4913 int d2 = decimals==-1 ? 1 : decimals;
4917 if(sfrequency.IsFloat()) {
4918 frequency = sfrequency.Atof();
4920 cout <<
"CWB::mdc::GetBurstNameLAL : Error : --frequency not defined or not a number!!! " << endl;
4929 cout <<
"CWB::mdc::GetBurstNameLAL : Error : --q not defined or not a number!!! " << endl;
4933 wfname = TString::Format(
"LAL_SGE%.*fQ%.*f",d1,frequency,d2,q);
4934 wfname.ReplaceAll(
".",
"d");
4935 if(decimals==-1) wfname.ReplaceAll(
"d0",
"");
4938 }
else if(name==
"BTLWNB") {
4940 int d1 = decimals==-1 ? 0 : decimals;
4941 int d2 = decimals==-1 ? 0 : decimals;
4942 int d3 = decimals==-1 ? 3 : decimals;
4946 if(sfrequency.IsFloat()) {
4947 frequency = sfrequency.Atof();
4949 cout <<
"CWB::mdc::GetBurstNameLAL : Error : --frequency not defined or not a number!!! " << endl;
4955 if(sbandwidth.IsFloat()) {
4956 bandwidth = sbandwidth.Atof();
4958 cout <<
"CWB::mdc::GetBurstNameLAL : Error : --bandwidth not defined or not a number!!! " << endl;
4964 if(sduration.IsFloat()) {
4965 duration = sduration.Atof();
4967 cout <<
"CWB::mdc::GetBurstNameLAL : Error : --duration not defined or not a number!!! " << endl;
4971 wfname = TString::Format(
"LAL_WNB%.*f_%.*f_%.*f",d1,frequency,d2,bandwidth,d3,duration);
4972 wfname.ReplaceAll(
".",
"d");
4975 }
else if(name==
"StringCusp") {
4977 int d1 = decimals==-1 ? 1 : decimals;
4981 if(sfrequency.IsFloat()) {
4982 frequency = sfrequency.Atof();
4984 cout <<
"CWB::mdc::GetBurstNameLAL : Error : --frequency not defined or not a number!!! " << endl;
4988 wfname = TString::Format(
"LAL_SC%.*f",d1,1000.*frequency);
4989 wfname.ReplaceAll(
".",
"d");
5022 cout <<
"CWB::mdc::e2cosi - Error : eccentricity must be defined in the range [0:1] " << endl;
5026 double E = e<0.9999999999 ?
e : 0.9999999999;
5027 double A = 1./sqrt(2-E*E);
5028 double B = A*sqrt(1-E*E);
5029 double x1 = (A+sqrt(A*A-B*B))/(B*B);
5030 double x2 = (A-sqrt(A*A-B*B))/(B*B);
5031 double x =
fabs(x1*B)>1 ? x2 : x1;
5051 cout <<
"CWB::mdc::cosi2e - Error : cosi must be defined in the range [-1:1] " << endl;
5055 double A = (1+cosi*cosi)/2;
5058 double C = sqrt(A*A+B*B);
5061 double E = sqrt(1-(B*B)/(A*A));
5078 if(fName.Contains(
".root")) {
5081 cout <<
"CWB::mdc::DumpLog - Error : dump to root needs injection object " << endl;
5085 TFile
froot(fName,
"RECREATE");
5086 if(label.Sizeof()>1) this->
Write(label);
else this->
Write(
"WaveMDC");
5088 cout <<
"CWB::mdc::DumpLog - inj entries written : " <<
inj_tree->GetEntries() << endl;
5093 if(fName.Contains(
".txt")) {
5101 if(append) out.open(fName.Data(),ios::app);
5102 else out.open(fName.Data(),
ios::out);
5103 if (!out.good()) {cout <<
"CWB::mdc::DumpLog - Error Opening File : " << fName.Data() << endl;
exit(1);}
5109 if(fName.Contains(
".lst")) {
5113 if(append) out.open(fName.Data(),ios::app);
5114 else out.open(fName.Data(),
ios::out);
5115 if (!out.good()) {cout <<
"CWB::mdc::DumpLog - Error Opening File : " << fName.Data() << endl;
exit(1);}
5117 out <<
"#gps\t\t" <<
"name\t\t" <<
"theta\t" <<
"phi\t" <<
"psi\t" 5118 <<
"rho\t" <<
"iota\t" <<
"hrss\t" <<
"ID\t" <<
"id\t" << endl;
5129 cout <<
"CWB::mdc::DumpLog - Error : file extention must be .root" << endl;
5154 cout <<
"CWB::mdc::GetInspiral - Error : Dummy method : network is not initialized " << endl;
5158 cout <<
"CWB::mdc::GetInspiral - Error : x.rate() != " <<
MDC_SAMPLE_RATE << endl;
5165 int gpsStartSec = x.
start();
5166 int gpsEndSec = gpsStartSec + x.
size()/x.
rate();
5169 INT4 numInjections = 0;
5170 SimInspiralTable *injections = NULL;
5171 TString injectionFile = GenInspiralXML(0., 0.,
false);
5174 lal_errhandler = LAL_ERR_EXIT;
5176 XLALSetSilentErrorHandler();
5179 numInjections = SimInspiralTableFromLIGOLw(&injections, injectionFile.Data(), gpsStartSec, gpsEndSec);
5181 if ( numInjections == 0 ) {
5182 fprintf( stderr,
"CWB::mdc::GetInspiral - Warning : No injections in specified time\n");
5185 fprintf( stdout,
"CWB::mdc::GetInspiral : Read %d injection(s) from the file '%s'\n",
5186 numInjections, injectionFile.Data() );
5194 fprintf(stdout,
"CWB::mdc::GetInspiral : Generating injection for: %s",ifo.Data());
5196 FindChirpInjectSignals(x, injections, ifo);
5201 SimInspiralTable *thisInj = NULL;
5202 for (thisInj = injections; thisInj; thisInj = thisInj->next)
5209 mdcList.push_back(log.Data());
5210 double gps = thisInj->geocent_end_time.gpsSeconds+1.e-9*thisInj->geocent_end_time.gpsNanoSeconds;
5217 while ( injections ) {
5218 SimInspiralTable *thisInj = NULL;
5219 thisInj = injections;
5220 injections = injections->next;
5224 LALCheckMemoryLeaks();
5232 CWB::mdc::GetInspiralLog(
TString inspName,
double FrameGPS, SimInspiralTable *thisInj) {
5242 cout <<
"CWB::mdc::GetInspiralLog - Error : Dummy method : network is not initialized " << endl;
5256 if (strncmp(thisInj->waveform,
"NumRel", LIGOMETA_WAVEFORM_MAX) == 0)
5257 sprintf(output,
"%s ", thisInj->numrel_data);
5261 sprintf(output,
"%s 1 ",output);
5263 sprintf(output,
"%s 0 ",output);
5265 sprintf(output,
"%s 0 ",output);
5267 sprintf(output,
"%s %g ",output, cos(thisInj->inclination));
5269 sprintf(output,
"%s %g ",output, thisInj->coa_phase);
5271 sprintf(output,
"%s %g ",output, cos(thisInj->latitude - LAL_PI_2));
5273 gmst = fmod(XLALGreenwichMeanSiderealTime(&thisInj->geocent_end_time), LAL_TWOPI);
5274 longitude = fmod(thisInj->longitude - gmst, LAL_TWOPI);
5276 longitude += LAL_TWOPI;
5277 sprintf(output,
"%s %g ",output, longitude);
5279 sprintf(output,
"%s %g ",output, thisInj->polarization);
5281 sprintf(output,
"%s %d ",output,
int(FrameGPS));
5283 sprintf(output,
"%s %d.%09d ",output, thisInj->geocent_end_time.gpsSeconds, thisInj->geocent_end_time.gpsNanoSeconds);
5285 sprintf(output,
"%s %s ",output, inspName.Data());
5287 sprintf(output,
"%s 0 ",output);
5289 sprintf(output,
"%s 0 ",output);
5291 sprintf(output,
"%s 0 ",output);
5301 latitude=thisInj->latitude;
5302 longitude=thisInj->longitude;
5303 polarization=thisInj->polarization;
5304 theta = LAL_PI_2-latitude;
5306 longitude = fmod(longitude - gmst, LAL_TWOPI);
5307 if (longitude < 0) longitude += LAL_TWOPI;
5308 phi = longitude > 0 ? longitude : 2*
TMath::Pi()+longitude;
5322 double D_R = sqrt(D->
Rv[0]*D->
Rv[0]+D->
Rv[1]*D->
Rv[1]+D->
Rv[2]*D->
Rv[2]);
5325 dV.SetTheta(D_theta*d2r);
5326 dV.SetPhi(D_phi*d2r);
5329 sV.SetTheta(theta*d2r);
5332 double cos_dVdS = dV.Dot(sV);
5333 double D_geocent_delay=-cos_dVdS*D_R/
speedlight;
5334 int tShift_sec =
int(
fabs(D_geocent_delay));
5335 int tShift_nsec =
int(1e9*(
fabs(D_geocent_delay)-tShift_sec));
5338 int end_time_gpsSeconds = thisInj->geocent_end_time.gpsSeconds;
5339 int end_time_gpsNanoSeconds = thisInj->geocent_end_time.gpsNanoSeconds;
5340 if(D_geocent_delay>=0) {
5341 end_time_gpsSeconds += tShift_sec;
5342 end_time_gpsNanoSeconds += tShift_nsec;
5343 if ( end_time_gpsNanoSeconds >= 1000000000 ) {
5344 end_time_gpsSeconds +=
INT_4S( end_time_gpsNanoSeconds / 1000000000 );
5345 end_time_gpsNanoSeconds = end_time_gpsNanoSeconds % 1000000000;
5348 end_time_gpsSeconds -= tShift_sec;
5349 if ( end_time_gpsNanoSeconds >= tShift_nsec ) {
5350 end_time_gpsNanoSeconds -= tShift_nsec;
5352 --end_time_gpsSeconds;
5353 end_time_gpsNanoSeconds += 1000000000 - tShift_nsec;
5358 double cosiota = cos(thisInj->inclination);
5359 double eff_dist = thisInj->distance / sqrt(pow(fp*(1.+pow(cosiota,2)),2)/4.+pow(fc*cosiota,2));
5361 sprintf(output,
"%s %s %d.%09d %e %e %g ",output, ifo.Data(), end_time_gpsSeconds,
5362 end_time_gpsNanoSeconds, fp, fc, eff_dist);
5367 f_isco = 205 * (20 / (thisInj->mass1 + thisInj->mass2));
5370 sprintf(output,
"%s insp ",output);
5371 sprintf(output,
"%s distance %g ",output, thisInj->distance);
5372 sprintf(output,
"%s mass1 %g ",output, thisInj->mass1);
5373 sprintf(output,
"%s mass2 %g ",output, thisInj->mass2);
5374 sprintf(output,
"%s mchirp %g ",output, thisInj->mchirp);
5375 sprintf(output,
"%s spin1 %g %g %g ",output, thisInj->spin1x, thisInj->spin1y, thisInj->spin1z);
5376 sprintf(output,
"%s spin2 %g %g %g ",output, thisInj->spin2x, thisInj->spin2y, thisInj->spin2z);
5377 sprintf(output,
"%s freq %g %g",output, thisInj->f_lower, f_isco);
5378 sprintf(output,
"%s redshift %g",output, thisInj->alpha3);
5379 sprintf(output,
"%s alpha1 %g",output, thisInj->alpha1);
5380 sprintf(output,
"%s alpha2 %g\n",output, thisInj->alpha2);
5413 TString option[nOpt] = {
"--help",
"--verbose",
"--user-tag",
"--write-compress",
5414 "--ipn-gps-time",
"--t-distr",
5415 "--write-sim-ring"};
5419 if(inspName.Sizeof()<=1) {
5420 cout <<
"CWB::mdc::SetInspiral - Error : inspName not declared" << endl;
5424 for(
int i=0;
i<nOpt;
i++) {
5425 if(inspOptions.Contains(option[
i])) {
5426 cout <<
"CWB::mdc::SetInspiral - Error : option " << option[
i].Data()
5427 <<
" not allowed in CWB" << endl;
5434 inspOptions.ReplaceAll(
"--approximant",
"--waveform");
5441 for(
int i=0;
i<token->GetEntries();
i++) {
5442 TObjString* otoken = (TObjString*)token->At(
i);
5443 TString stoken = otoken->GetString();
5444 if(stoken==
"--xml") bxml++;
5445 if(stoken==
"--waveform") bwaveform++;
5449 cout <<
"CWB::mdc::SetInspiral - Error : ";
5450 if(bwaveform>1) cout <<
"imultiple declaration of waveform option" << endl;
5451 if(!bwaveform) cout <<
"missing waveform option" << endl;
5460 cout <<
"CWB::mdc::SetInspiral - Read options ..." << endl;
5461 for(
int i=0;
i<token->GetEntries();
i++) {
5462 TObjString* otoken = (TObjString*)token->At(
i);
5463 TString stoken = otoken->GetString();
5464 if(stoken==
"--xml") {
5465 otoken = (TObjString*)token->At(
i+1);
5466 stoken = otoken->GetString();
5467 cout <<
"--xml = " << stoken.Data() << endl;
5472 if(stoken==
"--dir") {
5474 otoken = (TObjString*)token->At(
i+1);
5475 stoken = otoken->GetString();
5476 cout <<
"--dir = " << stoken.Data() << endl;
5481 inspOptions.ReplaceAll(
"--dir",
"");
5482 inspOptions.ReplaceAll(stoken,
"");
5484 if(stoken==
"--output") {
5485 otoken = (TObjString*)token->At(
i+1);
5486 stoken = otoken->GetString();
5487 cout <<
"--output = " << stoken.Data() << endl;
5490 inspOptions.ReplaceAll(
"--output",
"");
5491 inspOptions.ReplaceAll(stoken,
"");
5493 if(stoken==
"--dump") {
5494 otoken = (TObjString*)token->At(
i+1);
5495 stoken = otoken->GetString();
5496 cout <<
"--dump = " << stoken.Data() << endl;
5499 inspOptions.ReplaceAll(
"--dump",
"");
5500 inspOptions.ReplaceAll(stoken,
"");
5502 if(stoken==
"--time-step") {
5503 otoken = (TObjString*)token->At(
i+1);
5504 stoken = otoken->GetString();
5505 cout <<
"--time-step = " << stoken.Data() << endl;
5508 if(stoken==
"--time-interval") {
5509 otoken = (TObjString*)token->At(
i+1);
5510 stoken = otoken->GetString();
5511 cout <<
"--time-interval = " << stoken.Data() << endl;
5514 if(stoken==
"--waveform") {
5515 otoken = (TObjString*)token->At(
i+1);
5516 stoken = otoken->GetString();
5517 cout <<
"--waveform = " << stoken.Data() << endl;
5519 if(XLALGetApproximantFromString(stoken.Data()) == XLAL_FAILURE) {
5520 cout <<
"CWB::mdc::SetInspiral - waveform : \'" << stoken.Data()
5521 <<
"\' not allowed !!!" << endl;
5531 cout <<
"CWB::mdc::SetInspiral - Error : temporary dir not defined !!! Use --dir option" << endl;
5533 cout <<
" if this option is called from a config plugin add the following line to the inspOptions : " << endl;
5534 cout <<
" inspOptions+= \"--dir \"+TString(cfg->tmp_dir)+\" " << endl;
5535 cout <<
" before to call : 'MDC.SetInspiral(inspOptions);'" << endl << endl;
5536 cout <<
" and these lines to the plugin : " << endl;
5537 cout <<
" // export config object to the config plugin" << endl;
5538 cout <<
" sprintf(cmd,\"CWB::config* cfg = (CWB::config*)%p;\",cfg);" << endl;
5539 cout <<
" gROOT->ProcessLine(cmd);" << endl;
5540 cout <<
" before to execute the config plugin : 'cfg->configPlugin.Exec(NULL,&error);'" << endl;
5550 TString injectionFile = GenInspiralXML(0, 0,
false);
5551 cout <<
"Save inspiral injection file to " << inspDump.Data() << endl;
5552 if(
inspXML==
"") gSystem->Exec(
"mv "+injectionFile+
" "+inspDump);
5553 else gSystem->Exec(
"cp "+injectionFile+
" "+inspDump);
5558 if(inspOutput!=
"") {
5563 TString injectionFile = GenInspiralXML(0, 0,
false);
5564 cout <<
"Save inspiral injection file to " << inspOutput.Data() << endl;
5565 gSystem->Exec(
"cp "+injectionFile+
" "+inspOutput);
5570 if(
inspXML==
"") GenInspiralXML(900000000, 900000000,
true);
5579 CWB::mdc::GetInspiralOption(
TString inspOption) {
5586 if(inspOptions==
"")
return "";
5593 if (!in.good()) {cout <<
"CWB::mdc::GetInspiralOption - Error Opening xml File : " 5594 <<
inspXML.Data() << endl;gSystem->Exit(1);}
5598 in.getline(str,2048);
5599 if (!in.good())
break;
5600 if(
TString(str).Contains(
"inspinj")&&
TString(str).Contains(inspOption)) {
5602 if(token->GetEntries()!=5) {cout <<
"CWB::mdc::GetInspiralOption - bad line format : " 5603 << str << endl;gSystem->Exit(1);}
5604 TObjString* otoken = (TObjString*)token->At(4);
5605 TString stoken = otoken->GetString();
5606 stoken.ReplaceAll(
"\"",
"");
5616 for(
int i=0;
i<token->GetEntries();
i++) {
5617 TObjString* otoken = (TObjString*)token->At(
i);
5618 TString stoken = otoken->GetString();
5619 if(stoken.Contains(inspOption)) {
5620 otoken = (TObjString*)token->At(
i+1);
5621 stoken = otoken->GetString();
5622 stoken.ReplaceAll(
"\"",
"");
5636 CWB::mdc::GenInspiralXML(
int gps_start_time,
int gps_end_time,
bool rmFile) {
5650 if(gSystem->Getenv(
"LALINSPINJ_EXEC")==NULL) {
5651 cout <<
"CWB::mdc::GenInspiralXML - Error : environment LALINSPINJ_EXEC is not defined!!!" << endl;
exit(1);
5653 lalinspinj_exec=
TString(gSystem->Getenv(
"LALINSPINJ_EXEC"));
5660 cmdOptions +=
" --output "+
ofName;
5661 if(gps_start_time!=0) {
5662 cmdOptions +=
" --gps-start-time ";
5663 cmdOptions += gps_start_time;
5665 if(gps_end_time!=0) {
5666 cmdOptions +=
" --gps-end-time ";
5667 cmdOptions += gps_end_time;
5671 sprintf(cmd,
"%s %s",lalinspinj_exec.Data(),cmdOptions.Data());
5672 cout << cmd << endl;
5673 int error = gSystem->Exec(cmd);
5674 if(rmFile) gSystem->Exec(
TString(
"rm ")+ofName);
5679 sprintf(help,
"%s --help",lalinspinj_exec.Data());
5680 gSystem->Exec(help);
5682 cout <<
"CWB::mdc::GenInspiralXML - Error in exececuting :" << endl << endl;
5683 cout << cmd << endl << endl;
5703 gRandom->SetSeed(0);
5704 int rnID =
int(gRandom->Rndm(13)*1.e9);
5708 UserGroup_t*
uinfo = gSystem->GetUserInfo();
5712 if(mkdir) error=gSystem->Exec(
TString(
"mkdir -p ")+dir);
5715 cout <<
"CWB::mdc::GetTemporaryFileName - Error in exececuting :" << endl << endl;
5722 sprintf(fName,
"%s/%s_%d.%s",dir.Data(),tag.Data(),rnID,ext.Data());
5750 REAL8TimeSeries *
hp=NULL;
5751 REAL8TimeSeries *
hx=NULL;
5761 else pD =
new detector(const_cast<char*>(ifo.Data()));
5763 REAL8 location[3] = {pD->
Rv[0],pD->
Rv[1],pD->
Rv[2]};
5764 if(ifoId<0)
delete pD;
5767 SimInspiralTable *thisInj = NULL;
5768 for (thisInj = injections; thisInj; thisInj = thisInj->next)
5770 REAL8 deltaT = 1./w.
rate();
5774 Approximant approximant = (Approximant)XLALSimInspiralGetApproximantFromString(thisInj->waveform);
5786 if(approximant==NR_hdf5) {
5788 #if LAL_VERSION_MAJOR > 6 || (LAL_VERSION_MAJOR == 6 && \ 5789 (LAL_VERSION_MINOR > 18 || (LAL_VERSION_MINOR == 18 && \ 5790 (LAL_VERSION_MICRO > 0 || (LAL_VERSION_MICRO == 0 && LAL_VERSION_DEVEL >= 1))))) // LAL_VERSION >= 6.18.0.1 5793 LALDict *LALpars=XLALCreateDict();
5794 XLALSimInspiralWaveformParamsInsertNumRelData(LALpars, thisInj->numrel_data);
5797 LALValue* ModeArray = XLALSimInspiralCreateModeArray();
5799 int ell_min = thisInj->numrel_mode_min;
5800 int ell_max = thisInj->numrel_mode_max;
5802 if ((ell_min == 0) && (ell_max == 0)){ell_min=2;ell_max=LAL_SIM_L_MAX_MODE_ARRAY;}
5804 for (
int ell=ell_min; ell<=ell_max; ell++){
5805 XLALSimInspiralModeArrayActivateAllModesAtL(ModeArray, ell);
5808 XLALSimInspiralWaveformParamsInsertModeArray(LALpars, ModeArray);
5811 REAL8 f_min = XLALSimInspiralfLow2fStart(thisInj->f_lower, thisInj->amp_order, approximant);
5835 if( thisInj->amp_order!=0 ) {
5836 cout << endl <<
"CWB::mdc::FindChirpInjectSignals - " 5837 <<
"Error : amp_order must be 0 for NR waveform : " 5838 <<
"check inspiral parameters" << endl;
5843 int ret = XLALSimInspiralTD(&hp, &hx, thisInj->mass1*LAL_MSUN_SI, thisInj->mass2*LAL_MSUN_SI,
5844 thisInj->spin1x, thisInj->spin1y, thisInj->spin1z,
5845 thisInj->spin2x, thisInj->spin2y, thisInj->spin2z,
5846 thisInj->distance*LAL_PC_SI*1.0e6, thisInj->inclination,
5847 thisInj->coa_phase, 0., 0., 0.,
5849 LALpars, approximant);
5851 XLALDestroyDict(LALpars);
5853 if( ret==XLAL_FAILURE ) {
5854 cout << endl <<
"CWB::mdc::GetInspiral - " 5855 <<
"Error in XLALInspiralTDWaveformFromSimInspiral(1) : " 5856 <<
"check inspiral parameters" << endl;
5860 cout << endl <<
"CWB::mdc::GetInspiral - " 5861 <<
"Error - NR_hdf5 not supported for LAL_VERSION < 6.18.0.1" << endl;
5867 #if LAL_VERSION_MAJOR > 6 || (LAL_VERSION_MAJOR == 6 && \ 5868 (LAL_VERSION_MINOR > 16 || (LAL_VERSION_MINOR == 16 && \ 5869 LAL_VERSION_MICRO >= 1 ))) // LAL_VERSION >= 6.16.1 5873 if(source(0,5)==
"#CWB:") {
5875 TString version=source(source.Index(
":",4)+1,source.Index(
":",5)-source.Index(
":",4)-1);
5877 if(version.Atoi()<1) {
5878 cout << endl <<
"CWB::mdc::FindChirpInjectSignals - " 5879 <<
"Error: outdated version of XML produced by CWB::mdc::Posterior2XML, must be >= 1" << endl << endl;
5884 #if LAL_VERSION_MAJOR == 6 && LAL_VERSION_MINOR == 16 && LAL_VERSION_MICRO == 1 // used only for tagged version lalinference_o2 used for O2 posteriors 5888 double lambda1 = source(0,5)==
"#CWB:" ? thisInj->alpha : 0;
5889 double lambda2 = source(0,5)==
"#CWB:" ? thisInj->beta : 0;
5892 REAL8 f_min = XLALSimInspiralfLow2fStart(thisInj->f_lower, thisInj->amp_order, approximant);
5895 int pn_order = XLALSimInspiralGetPNOrderFromString(thisInj->waveform);
5898 int ret = XLALSimInspiralTD(&hp, &hx, thisInj->coa_phase, deltaT,
5899 thisInj->mass1*LAL_MSUN_SI, thisInj->mass2*LAL_MSUN_SI,
5900 thisInj->spin1x, thisInj->spin1y, thisInj->spin1z,
5901 thisInj->spin2x, thisInj->spin2y, thisInj->spin2z,
5902 f_min, thisInj->f_final,
5903 thisInj->distance*LAL_PC_SI*1.0e6, 0., thisInj->inclination,
5904 lambda1, lambda2, NULL, NULL,
5905 thisInj->amp_order, pn_order, approximant);
5909 int ret = this->XLALInspiralTDWaveformFromSimInspiral(&hp, &hx, thisInj, deltaT);
5912 if( ret==XLAL_FAILURE ) {
5913 cout << endl <<
"CWB::mdc::GetInspiral - " 5914 <<
"Error in XLALInspiralTDWaveformFromSimInspiral(2) : " 5915 <<
"check inspiral parameters" << endl;
5919 cout << endl <<
"CWB::mdc::GetInspiral - " 5920 <<
"Error - LAL_VERSION < 6.16.1 not supported " << endl;
5931 XLALGPSAddGPS(&hp->epoch, &thisInj->geocent_end_time);
5932 XLALGPSAddGPS(&hx->epoch, &thisInj->geocent_end_time);
5935 wat::Time gpsinj(hp->epoch.gpsSeconds,hp->epoch.gpsNanoSeconds);
5939 double binoff = fmod(gpsoff.
GetDouble(),deltaT);
5941 bininj+=hp->data->length;
5942 if(bininj>w.
size()) {
5943 cout <<
"CWB::mdc::FindChirpInjectSignals - Warning " << endl;
5944 cout <<
" signal length " << hp->data->length*deltaT
5945 <<
" (sec) is too high, signal is cutted to fit the buffer length" << endl;
5946 cout <<
" Decrease the injection rate and/or reduce the waveform length" << endl;
5950 gmst = fmod(XLALGreenwichMeanSiderealTime(&thisInj->geocent_end_time), LAL_TWOPI);
5952 latitude=thisInj->latitude;
5953 longitude=thisInj->longitude;
5954 polarization=thisInj->polarization;
5955 theta = LAL_PI_2-latitude;
5957 longitude = fmod(longitude - gmst, LAL_TWOPI);
5958 if (longitude < 0) longitude += LAL_TWOPI;
5959 phi = longitude > 0 ? longitude : 2*
TMath::Pi()+longitude;
5967 int length = hp->data->length<w.
size() ? hp->data->length : w.
size();
5969 for(
int i=bininj;
i>=0;
i--) {
5970 int j=length-(bininj-
i)-1;
5972 z[
i] = fp*hp->
data->data[
j]+fx*hx->data->data[
j];
5975 int simulation_id = thisInj->bandpass;
5976 wat::Time wtime(thisInj->geocent_end_time.gpsSeconds, thisInj->geocent_end_time.gpsNanoSeconds);
5977 CalibrateInspiral(z, ifo, wtime.
GetDouble(), simulation_id);
5981 double tShift = XLALTimeDelayFromEarthCenter(location, thisInj->longitude, thisInj->latitude, &hp->epoch);
5987 XLALDestroyREAL8TimeSeries( hp );
5988 XLALDestroyREAL8TimeSeries( hx );
5989 hp = NULL; hx = NULL;
5998 CWB::mdc::GetInspiral(
TString pol,
int gps_start_time,
int gps_end_time) {
6013 if(pol!=
"hp" && pol!=
"hx") {
6014 fprintf( stderr,
"CWB::mdc::GetInspiral - Warning : wrong polarization (enter hp or hx)\n");
6019 int gpsStartSec = gps_start_time;
6020 int gpsEndSec = gps_end_time;
6023 INT4 numInjections = 0;
6024 SimInspiralTable *injections = NULL;
6025 TString injectionFile = GenInspiralXML(0., 0.,
false);
6028 lal_errhandler = LAL_ERR_EXIT;
6030 XLALSetSilentErrorHandler();
6033 numInjections = SimInspiralTableFromLIGOLw(&injections, injectionFile.Data(), gpsStartSec, gpsEndSec);
6035 if ( numInjections == 0 ) {
6036 fprintf( stderr,
"CWB::mdc::GetInspiral - Warning : No injections in specified time\n");
6039 fprintf( stdout,
"CWB::mdc::GetInspiral : Read %d injection(s) from the file '%s'\n",
6040 numInjections, injectionFile.Data() );
6043 REAL8TimeSeries *
hp=NULL;
6044 REAL8TimeSeries *
hx=NULL;
6047 SimInspiralTable *thisInj = NULL;
6048 for (thisInj = injections; thisInj; thisInj = thisInj->next)
6050 REAL8 deltaT = 1./w.
rate();
6054 Approximant approximant = (Approximant)XLALSimInspiralGetApproximantFromString(thisInj->waveform);
6066 if(approximant==NR_hdf5) {
6068 #if LAL_VERSION_MAJOR > 6 || (LAL_VERSION_MAJOR == 6 && \ 6069 (LAL_VERSION_MINOR > 18 || (LAL_VERSION_MINOR == 18 && \ 6070 (LAL_VERSION_MICRO > 0 || (LAL_VERSION_MICRO == 0 && LAL_VERSION_DEVEL >= 1))))) // LAL_VERSION >= 6.18.0.1 6073 LALDict *LALpars=XLALCreateDict();
6074 XLALSimInspiralWaveformParamsInsertNumRelData(LALpars, thisInj->numrel_data);
6077 LALValue* ModeArray = XLALSimInspiralCreateModeArray();
6079 int ell_min = thisInj->numrel_mode_min;
6080 int ell_max = thisInj->numrel_mode_max;
6082 if ((ell_min == 0) && (ell_max == 0)){ell_min=2;ell_max=LAL_SIM_L_MAX_MODE_ARRAY;}
6084 for (
int ell=ell_min; ell<=ell_max; ell++){
6085 XLALSimInspiralModeArrayActivateAllModesAtL(ModeArray, ell);
6089 XLALSimInspiralWaveformParamsInsertModeArray(LALpars, ModeArray);
6092 REAL8 f_min = XLALSimInspiralfLow2fStart(thisInj->f_lower, thisInj->amp_order, approximant);
6116 if( thisInj->amp_order!=0 ) {
6117 cout << endl <<
"CWB::mdc::GetInspiral - " 6118 <<
"Error : amp_order must be 0 for NR waveform : " 6119 <<
"check inspiral parameters" << endl;
6124 int ret = XLALSimInspiralTD(&hp, &hx, thisInj->mass1*LAL_MSUN_SI, thisInj->mass2*LAL_MSUN_SI,
6125 thisInj->spin1x, thisInj->spin1y, thisInj->spin1z,
6126 thisInj->spin2x, thisInj->spin2y, thisInj->spin2z,
6127 thisInj->distance*LAL_PC_SI*1.0e6, thisInj->inclination,
6128 thisInj->coa_phase, 0., 0., 0.,
6130 LALpars, approximant);
6132 XLALDestroyDict(LALpars);
6134 if( ret==XLAL_FAILURE ) {
6135 cout << endl <<
"CWB::mdc::GetInspiral - " 6136 <<
"Error in XLALInspiralTDWaveformFromSimInspiral(3) : " 6137 <<
"check inspiral parameters" << endl;
6141 cout << endl <<
"CWB::mdc::GetInspiral - " 6142 <<
"Error - NR_hdf5 not supported for LAL_VERSION < 6.18.0.1" << endl;
6148 #if LAL_VERSION_MAJOR > 6 || (LAL_VERSION_MAJOR == 6 && \ 6149 (LAL_VERSION_MINOR > 16 || (LAL_VERSION_MINOR == 16 && \ 6150 LAL_VERSION_MICRO >= 1 ))) // LAL_VERSION >= 6.16.1 6154 if(source(0,5)==
"#CWB:") {
6156 TString version=source(source.Index(
":",4)+1,source.Index(
":",5)-source.Index(
":",4)-1);
6158 if(version.Atoi()<1) {
6159 cout << endl <<
"CWB::mdc::FindChirpInjectSignals - " 6160 <<
"Error: outdated version of XML produced by CWB::mdc::Posterior2XML, must be >= 1" << endl << endl;
6165 #if LAL_VERSION_MAJOR == 6 && LAL_VERSION_MINOR == 16 && LAL_VERSION_MICRO == 1 // used only for tagged version lalinference_o2 used for O2 posteriors 6169 double lambda1 = source(0,5)==
"#CWB:" ? thisInj->alpha : 0;
6170 double lambda2 = source(0,5)==
"#CWB:" ? thisInj->beta : 0;
6173 REAL8 f_min = XLALSimInspiralfLow2fStart(thisInj->f_lower, thisInj->amp_order, approximant);
6176 int pn_order = XLALSimInspiralGetPNOrderFromString(thisInj->waveform);
6179 int ret = XLALSimInspiralTD(&hp, &hx, thisInj->coa_phase, deltaT,
6180 thisInj->mass1*LAL_MSUN_SI, thisInj->mass2*LAL_MSUN_SI,
6181 thisInj->spin1x, thisInj->spin1y, thisInj->spin1z,
6182 thisInj->spin2x, thisInj->spin2y, thisInj->spin2z,
6183 f_min, thisInj->f_final,
6184 thisInj->distance*LAL_PC_SI*1.0e6, 0., thisInj->inclination,
6185 lambda1, lambda2, NULL, NULL,
6186 thisInj->amp_order, pn_order, approximant);
6190 int ret = this->XLALInspiralTDWaveformFromSimInspiral(&hp, &hx, thisInj, deltaT);
6193 if( ret==XLAL_FAILURE ) {
6194 cout << endl <<
"CWB::mdc::GetInspiral - " 6195 <<
"Error in XLALInspiralTDWaveformFromSimInspiral : " 6196 <<
"check inspiral parameters" << endl;
6200 cout << endl <<
"CWB::mdc::GetInspiral - " 6201 <<
"Error - LAL_VERSION < 6.16.1 not supported " << endl;
6212 XLALGPSAddGPS(&hp->epoch, &thisInj->geocent_end_time);
6213 XLALGPSAddGPS(&hx->epoch, &thisInj->geocent_end_time);
6217 w.
resize(hp->data->length);
6218 wat::Time gps(hp->epoch.gpsSeconds,hp->epoch.gpsNanoSeconds);
6223 w.
resize(hx->data->length);
6224 wat::Time gps(hx->epoch.gpsSeconds,hx->epoch.gpsNanoSeconds);
6229 XLALDestroyREAL8TimeSeries(hp);
6230 XLALDestroyREAL8TimeSeries(hx);
6231 hp = NULL; hx = NULL;
6238 while ( injections ) {
6239 SimInspiralTable *thisInj = NULL;
6240 thisInj = injections;
6241 injections = injections->next;
6245 LALCheckMemoryLeaks();
6264 CWB::mdc::XLALInspiralTDWaveformFromSimInspiral(
6265 REAL8TimeSeries **hplus,
6266 REAL8TimeSeries **hcross,
6267 SimInspiralTable *thisRow,
6271 #if !(LAL_VERSION_MAJOR == 6 && LAL_VERSION_MINOR == 16 && LAL_VERSION_MICRO == 1) // not used for tagged version lalinference_o2 used for O2 posteriors 6274 Approximant approximant;
6275 REAL8 phi0 = thisRow->coa_phase;
6276 REAL8
m1 = thisRow->mass1 * LAL_MSUN_SI;
6277 REAL8
m2 = thisRow->mass2 * LAL_MSUN_SI;
6278 REAL8 S1x = thisRow->spin1x;
6279 REAL8 S1y = thisRow->spin1y;
6280 REAL8 S1z = thisRow->spin1z;
6281 REAL8 S2x = thisRow->spin2x;
6282 REAL8 S2y = thisRow->spin2y;
6283 REAL8 S2z = thisRow->spin2z;
6284 REAL8 f_min = thisRow->f_lower;
6285 REAL8 f_ref = thisRow->f_final;
6287 REAL8
distance = thisRow->distance * LAL_PC_SI * 1e6;
6288 REAL8 inclination = thisRow->inclination;
6289 LALDict *params = XLALCreateDict();
6290 XLALSimInspiralWaveformParamsInsertPNAmplitudeOrder(params, thisRow->amp_order);
6293 approximant = (Approximant)XLALSimInspiralGetApproximantFromString(thisRow->waveform);
6294 if ( (
int) approximant == XLAL_FAILURE)
6295 XLAL_ERROR(XLAL_EFUNC);
6297 if (approximant == NR_hdf5) {
6298 char *filepath = thisRow->numrel_data;
6299 XLALSimInspiralWaveformParamsInsertNumRelData(params, filepath);
6302 order = (LALPNOrder)XLALSimInspiralGetPNOrderFromString(thisRow->waveform);
6303 if ( (
int) order == XLAL_FAILURE)
6304 XLAL_ERROR(XLAL_EFUNC);
6305 XLALSimInspiralWaveformParamsInsertPNPhaseOrder(params, order);
6311 #ifdef NEW_CODE_FOR_SEOBNRv4P_AND_SEOBNRv4PHM 6312 if (approximant == SEOBNRv4P || approximant == SEOBNRv4PHM) {
6313 ret = XLALSimInspiralChooseTDWaveform(hplus, hcross, m1, m2, S1x, S1y, S1z, S2x, S2y, S2z, distance, inclination, phi0, 0.0, 0.0, 0.0, deltaT, f_min, f_ref, params, approximant);
6314 XLALDestroyDict(params);
6315 if( ret == XLAL_FAILURE )
6316 XLAL_ERROR(XLAL_EFUNC);
6318 const double extra_time_fraction = 0.1;
6319 const double extra_cycles = 3.0;
6322 double tchirp = XLALSimInspiralChirpTimeBound(f_min, m1, m2, S1z, S2z);
6327 double textra = extra_cycles / f_min;
6331 XLALSimInspiralTDConditionStage1(*hplus, *hcross, extra_time_fraction * tchirp + textra, f_min);
6338 double fisco = 1.0 / (pow(6.0, 1.5) * LAL_PI * (m1 +
m2) * LAL_MTSUN_SI / LAL_MSUN_SI);
6339 XLALSimInspiralTDConditionStage2(*hplus, *hcross, f_min, fisco);
6341 ret = XLALSimInspiralTD(hplus, hcross, m1, m2, S1x, S1y, S1z, S2x, S2y, S2z, distance, inclination, phi0, 0.0, 0.0, 0.0, deltaT, f_min, f_ref, params, approximant);
6342 XLALDestroyDict(params);
6343 if( ret == XLAL_FAILURE )
6344 XLAL_ERROR(XLAL_EFUNC);
6366 ret = XLALSimInspiralTD(hplus, hcross, m1, m2, S1x, S1y, S1z, S2x, S2y, S2z, distance, inclination, phi0, 0.0, 0.0, 0.0, deltaT, f_min, f_ref, params, approximant);
6367 XLALDestroyDict(params);
6368 if( ret == XLAL_FAILURE )
6369 XLAL_ERROR(XLAL_EFUNC);
6371 #endif // endif lalinference_o2 6373 return XLAL_SUCCESS;
6408 double gps_start_time=0, gps_stop_time=0, time_step = 0;
6409 int time_nsec = -1.;
6411 int ninjections = 0;
6414 if(sgps_start_time.IsDigit()) {
6415 gps_start_time = sgps_start_time.Atoi();
6416 }
else if(sgps_start_time.IsFloat()) {
6417 gps_start_time = sgps_start_time.Atoi();
6418 sgps_start_time.Remove(0,sgps_start_time.Last(
'.'));
6419 time_nsec =
int(1.e9*sgps_start_time.Atof());
6420 }
else if(sgps_start_time!=
"") {
6421 cout << endl <<
"CWB::mdc::Posterior2XML - " 6422 <<
"Error : gps_start_time must be an interger number" << endl;
6427 if(sgps_stop_time.IsDigit()) {
6428 gps_stop_time = sgps_stop_time.Atoi();
6429 }
else if(sgps_stop_time!=
"") {
6430 cout << endl <<
"CWB::mdc::Posterior2XML - " 6431 <<
"Error : gps_stop_time must be an interger number" << endl;
6436 if(stime_step.IsDigit()) {
6437 time_step = stime_step.Atoi();
6438 }
else if(stime_step!=
"") {
6439 cout << endl <<
"CWB::mdc::Posterior2XML - " 6440 <<
"Error : time_step must be an interger number" << endl;
6445 if(sseed.IsDigit()) seed = sseed.Atoi();
6448 if(sninjections.IsDigit()) ninjections = sninjections.Atof();
6453 if(ssource==
"") ssource=
"NULL";
6454 if(ssource.Sizeof()>20) {
6455 cout << endl <<
"CWB::mdc::Posterior2XML - " 6456 <<
"Error : max source length is 20" << endl;
6461 TString extension = clb_file(clb_file.Last(
'.'),clb_file.Sizeof()-clb_file.Last(
'.')-1);
6462 if(extension!=
".clb") {
6463 cout <<
"CWB::mdc::Posterior2XML - Error : Calibration file extension must be .clb" << endl;
6468 if(sample_opt!=
"" && sample_opt!=
"map" && sample_opt!=
"maxl" && !sample_opt.IsDigit()) {
6469 cout << endl <<
"CWB::mdc::Posterior2XML - " 6470 <<
"Error : wrong --sample option type, must be map/maxl or an integer number >=0" << endl;
6475 if(taper!=
"" && taper!=
"TAPER_NONE" && taper!=
"TAPER_START" && taper!=
"TAPER_END" && taper!=
"TAPER_STARTEND") {
6476 cout << endl <<
"CWB::mdc::Posterior2XML - " 6477 <<
"Error : wrong --taper option type, must be TAPER_NONE,TAPER_START,TAPER_END,TAPER_STARTEND" << endl;
6483 if(sf_lower.IsFloat()) {
6484 f_lower = sf_lower.Atof();
6485 }
else if(sf_lower!=
"") {
6486 cout << endl <<
"CWB::mdc::Posterior2XML - " 6487 <<
"Error : f_lower must be a float number" << endl;
6493 if(sf_ref.IsFloat()) {
6494 f_ref = sf_ref.Atof();
6495 }
else if(sf_ref!=
"") {
6496 cout << endl <<
"CWB::mdc::Posterior2XML - " 6497 <<
"Error : f_ref must be a float number" << endl;
6502 int pn_order_opt=-1;
6503 if(spn_order.IsDigit()) {
6504 pn_order_opt = spn_order.Atoi();
6505 }
else if(spn_order!=
"") {
6506 cout << endl <<
"CWB::mdc::Posterior2XML - " 6507 <<
"Error : pn_order must be an interger number" << endl;
6510 if(pn_order_opt!=-1 && swaveform!=
"") {
6511 cout << endl <<
"CWB::mdc::Posterior2XML - " 6512 <<
"Error : only one of the following options can be defined: --pn_order, --waveform" << endl;
6517 int amp_order_opt=-1;
6518 if(samp_order.IsDigit()) {
6519 amp_order_opt = samp_order.Atoi();
6520 }
else if(samp_order!=
"") {
6521 cout << endl <<
"CWB::mdc::Posterior2XML - " 6522 <<
"Error : amp_order must be an interger number" << endl;
6526 if(gps_start_time>0) {
6527 if(gps_stop_time<gps_start_time) {
6528 cout << endl <<
"CWB::mdc::Posterior2XML - " 6529 <<
"Error : if gps_start_time>0 then gps_stop_time must be >= gps_start_time" << endl;
6533 cout << endl <<
"CWB::mdc::Posterior2XML - " 6534 <<
"Error : if gps_start_time>0 then time_step must be > 0" << endl;
6538 cout << endl <<
"CWB::mdc::Posterior2XML - " 6539 <<
"Error : if gps_start_time>0 then ninjections can not be declared" << endl;
6543 ninjections =
int((gps_stop_time-gps_start_time)/time_step);
6545 if(ninjections<=0) {
6546 cout << endl <<
"CWB::mdc::Posterior2XML - " 6547 <<
"Error : if gps_start_time<=0 then --ninjections option must be > 0" << endl;
6553 in.open(sampleFile.Data(),
ios::in);
6555 cout << endl <<
"CWB::mdc::Posterior2XML - " 6556 <<
"Error Opening File : " << sampleFile.Data() << endl;
6562 in.getline(hline,4*1024);
6569 sline.ReplaceAll(
"#",
"\t");
6570 sline.ReplaceAll(
" ",
"\t");
6571 sline.ReplaceAll(
"\treference_frequency\t",
"\tf_ref\t");
6572 sline.ReplaceAll(
"\tlambda_1\t",
"\tlambda1\t");
6573 sline.ReplaceAll(
"\tlambda_2\t",
"\tlambda2\t");
6574 sline.ReplaceAll(
"\tluminosity_distance\t",
"\tdist\t");
6575 sline.ReplaceAll(
"\tmass_ratio\t",
"\tq\t");
6576 sline.ReplaceAll(
"\tchirp_mass\t",
"\tmc\t");
6577 sline.ReplaceAll(
"\trightascension\t",
"\tra\t");
6578 sline.ReplaceAll(
"\tdeclination\t",
"\tdec\t");
6579 sline.ReplaceAll(
"\tcos_theta_jn\t",
"\tcostheta_jn\t");
6580 sline.ReplaceAll(
"\ttilt_1\t",
"\ttilt1\t");
6581 sline.ReplaceAll(
"\ttilt_2\t",
"\ttilt2\t");
6582 sline.ReplaceAll(
"\tcos_tilt_1\t",
"\tcostilt1\t");
6583 sline.ReplaceAll(
"\tcos_tilt_2\t",
"\tcostilt2\t");
6584 sline.ReplaceAll(
"\tphi_12\t",
"\tphi12\t");
6585 sline.ReplaceAll(
"\ta_1\t",
"\ta1\t");
6586 sline.ReplaceAll(
"\ta_2\t",
"\ta2\t");
6587 sline.ReplaceAll(
"\tspin_1z\t",
"\ta1z\t");
6588 sline.ReplaceAll(
"\tspin_2z\t",
"\ta2z\t");
6589 sline.ReplaceAll(
"\ta_1\t",
"\ta1\t");
6590 sline.ReplaceAll(
"\ta_2\t",
"\ta2\t");
6591 sline.ReplaceAll(
"\tmarginalized_phase\t",
"\tphase_maxl\t");
6592 sline.ReplaceAll(
"\tgeocent_time\t",
"\ttime\t");
6593 sline.ReplaceAll(
"\tmarginalized_geocent_time\t",
"\ttime_maxl\t");
6594 sline.ReplaceAll(
"\tlog_likelihood\t",
"\tlogl\t");
6595 sline.ReplaceAll(
"\tlog_prior\t",
"\tlogprior\t");
6596 sline.ReplaceAll(
"\tH1_time\t",
"\th1_end_time\t");
6597 sline.ReplaceAll(
"\tL1_time\t",
"\tl1_end_time\t");
6598 sline.ReplaceAll(
"\tV1_time\t",
"\tv1_end_time\t");
6599 sline.ReplaceAll(
"\tG1_time\t",
"\tg1_end_time\t");
6600 sline=sline(1,sline.Sizeof()-3);
6602 cout << endl <<
"CWB::mdc::Posterior2XML - header converted " << endl;
6603 cout << endl <<
"before convertion:" << endl << hline << endl;
6604 cout << endl <<
"after convertion:" << endl << sline << endl << endl;
6606 strcpy(hline,sline.Data());
6608 std::map<TString, int> hsample;
6610 TString first_sheader=
"";
6611 for(
int j=0;
j<token->GetEntries();
j++) {
6612 TString sheader = ((TObjString*)token->At(
j))->
GetString();
6613 sheader.ReplaceAll(
" " ,
"");
6618 first_sheader=sheader;
6620 hsample[sheader] =
j;
6623 hsample[first_sheader] = token->GetEntries();
6632 std::vector<vector<double> > sample;
6635 int hsize=token->GetEntries();
6636 std::vector<double> vsample(hsample.size());
6638 for(
int i=0;i<hsize;i++) in >> vsample[
i];
6639 vsample[hsize]=vsample[0];
6640 if(!in.good())
break;
6641 sample.push_back(vsample);
6645 cout << endl <<
"CWB::mdc::Posterior2XML - " 6646 <<
"Read\t" << nsample <<
"\tsamples from posterior samples file" << endl;
6650 if(sample_opt==
"map") {
6651 double post_max=-1
e+20;
6652 for(
int i=0;
i<sample.size();
i++) {
6654 if(hsample[
"logpost"]) {
6655 post = sample[
i][hsample[
"logpost"]];
6656 }
else if(hsample[
"logprior"] && hsample[
"logl"]) {
6657 post = sample[
i][hsample[
"logprior"]] + sample[
i][hsample[
"logl"]];
6658 }
else if(hsample[
"post"]) {
6659 post = sample[
i][hsample[
"post"]];
6661 cout << endl <<
"CWB::mdc::Posterior2XML - " 6662 <<
"Error : none of values used by map option are present in the sample header" << endl;
6665 if(post>post_max) {post_max=post;mapSampleID=
i;}
6667 cout << endl <<
"CWB::mdc::Posterior2XML - found map sample id = " << mapSampleID << endl << endl;
6671 int maxlSampleID=-1;
6672 if(sample_opt==
"maxl") {
6673 double maxl_max=-1
e+20;
6674 for(
int i=0;
i<sample.size();
i++) {
6676 if(hsample[
"logl"]) {
6677 maxl = sample[
i][hsample[
"logl"]];
6679 cout << endl <<
"CWB::mdc::Posterior2XML - " 6680 <<
"Error : logl is not present in the sample header" << endl;
6683 if(maxl>maxl_max) {maxl_max=maxl;maxlSampleID=
i;}
6685 cout << endl <<
"CWB::mdc::Posterior2XML - found maxl sample id = " << maxlSampleID << endl << endl;
6690 if(sample_opt.IsDigit()) {
6691 SampleID=sample_opt.Atoi();
6692 if(SampleID<0 || SampleID>=sample.size()) {
6693 cout << endl <<
"CWB::mdc::Posterior2XML - " 6694 <<
"Error : sample index is not valid. Allowed valued for this posterior sample file are [0," << sample.size()-1 <<
"]" << endl;
6702 if(!clb_file.Contains(
".clb")) {
6703 cout << endl <<
"CWB::mdc::Posterior2XML - " 6704 <<
"Error : calibration file extention must be .clb" << endl;
6707 clb.open(clb_file.Data(),
ios::out);
6709 cout << endl <<
"CWB::mdc::Posterior2XML - " 6710 <<
"Error Opening Output Calibration File : " << clb_file.Data() << endl;
6715 std::map<TString, int>::iterator iter;
6716 for (iter=hsample.begin(); iter!=hsample.end(); iter++) {
6717 if(
TString(iter->first).Contains(
"spcal")) clb << iter->first <<
"\t";
6719 clb <<
"simulation_id";
6723 static LALStatus status;
6724 LIGOLwXMLStream xmlfp;
6725 MetadataTable injections;
6726 SimInspiralTable *simTable;
6728 gps_start_time =
int(gps_start_time/time_step)*time_step;
6732 wat::Time currentGpsTime(gps_start_time);
6736 simTable = injections.simInspiralTable = (SimInspiralTable *) calloc( 1,
sizeof(SimInspiralTable) );
6739 TRandom3 random(abs(seed));
6743 int simulation_id = 0;
6744 int initialNanoSeconds=-1;
6747 int id = (seed>=0) ? gRandom->Uniform(0,sample.size()-1) : ninj;
6748 if(sample_opt==
"map")
id=mapSampleID;
6749 if(sample_opt==
"maxl")
id=maxlSampleID;
6750 if(sample_opt.IsDigit())
id=SampleID;
6753 std::map<TString, double> psample = GetPsample(hsample, sample[
id], f_lower, f_ref);
6754 if(psample.size()==0) {
6756 if(seed>=0)
continue;
6757 else {
if(ninj<sample.size())
continue;
else break;}
6761 char source_tag[30];
sprintf(source_tag,
"#CWB:1:%s",ssource.Data());
6762 memcpy( simTable->source, source_tag,
sizeof(CHAR) * LIGOMETA_SOURCE_MAX );
6765 double time = psample[
"time"];
6769 if(gps_start_time>0) {
6770 simTable->geocent_end_time.gpsSeconds = currentGpsTime.
GetSec();
6772 simTable->geocent_end_time.gpsNanoSeconds = time_nsec;
6774 simTable->geocent_end_time.gpsNanoSeconds = iwtime.
GetNSec();
6776 if(initialNanoSeconds==-1) initialNanoSeconds=iwtime.
GetNSec();
6777 int deltaNSec =
int(iwtime.
GetNSec())-initialNanoSeconds;
6778 if(deltaNSec> 500000000) simTable->geocent_end_time.gpsSeconds-=1;
6779 if(deltaNSec<-500000000) simTable->geocent_end_time.gpsSeconds+=1;
6782 simTable->geocent_end_time.gpsSeconds = iwtime.
GetSec();
6783 simTable->geocent_end_time.gpsNanoSeconds = iwtime.
GetNSec();
6785 wat::Time owtime(simTable->geocent_end_time.gpsSeconds, simTable->geocent_end_time.gpsNanoSeconds);
6788 wat::Time h_time(psample[
"h1_end_time"]);
6789 simTable->h_end_time.gpsSeconds = h_time.
GetSec();
6790 simTable->h_end_time.gpsNanoSeconds = h_time.
GetNSec();
6792 wat::Time l_time(psample[
"l1_end_time"]);
6793 simTable->l_end_time.gpsSeconds = l_time.
GetSec();
6794 simTable->l_end_time.gpsNanoSeconds = l_time.
GetNSec();
6796 wat::Time v_time(psample[
"v1_end_time"]);
6797 simTable->v_end_time.gpsSeconds = v_time.
GetSec();
6798 simTable->v_end_time.gpsNanoSeconds = v_time.
GetNSec();
6800 wat::Time g_time(psample[
"g1_end_time"]);
6801 simTable->g_end_time.gpsSeconds = g_time.
GetSec();
6802 simTable->g_end_time.gpsNanoSeconds = g_time.
GetNSec();
6805 simTable->longitude = psample[
"ra"];
6806 simTable->latitude = psample[
"dec"];
6807 simTable->polarization = psample[
"psi"];
6808 simTable->distance = psample[
"dist"];
6810 if(gps_start_time>0) {
6812 double gps_source = time;
6813 double phi_source = sm.
RA2phi(simTable->longitude*180./
TMath::Pi(), gps_source);
6815 if(ra_source>180) ra_source-=360.;
6816 simTable->longitude = ra_source*
TMath::Pi()/180.;
6820 simTable->spin1x = psample[
"s1x"];
6821 simTable->spin1y = psample[
"s1y"];
6822 simTable->spin1z = psample[
"s1z"];
6823 simTable->spin2x = psample[
"s2x"];
6824 simTable->spin2y = psample[
"s2y"];
6825 simTable->spin2z = psample[
"s2z"];
6828 simTable->mass1 = psample[
"m1"];
6829 simTable->mass2 = psample[
"m2"];
6830 simTable->mchirp = psample[
"mchirp"];
6833 int approx = psample[
"approx"];
6834 int pn_order = psample[
"pn_order"];
6835 if(pn_order_opt>=0) pn_order = pn_order_opt;
6838 sprintf(waveform,
"%s",XLALSimInspiralGetStringFromApproximant((Approximant)approx));
6840 sprintf(waveform,
"%s%s",XLALSimInspiralGetStringFromApproximant((Approximant)approx),XLALSimInspiralGetStringFromPNOrder((LALPNOrder)pn_order));
6842 if(swaveform!=
"")
sprintf(waveform,
"%s",swaveform.Data());
6843 memcpy( simTable->waveform, waveform,
sizeof(CHAR) * LIGOMETA_WAVEFORM_MAX );
6846 simTable->f_lower = psample[
"flow"];
6847 if(f_lower>=0) simTable->f_lower = f_lower;
6848 simTable->f_final = psample[
"f_ref"];
6849 if(f_ref>=0) simTable->f_final = f_ref;
6850 simTable->amp_order =
int(psample[
"amp_order"]);
6851 if(amp_order_opt>=0) simTable->amp_order = amp_order_opt;
6852 simTable->inclination = psample[
"iota"];
6853 simTable->coa_phase = psample[
"phase"];
6856 simTable->alpha = psample[
"lambda1"];
6857 simTable->beta = psample[
"lambda2"];
6862 simTable->alpha1 = psample[
"logl"];
6863 simTable->alpha2 = psample[
"post"];
6870 if(
TString(simTable->taper)==
"") memcpy( simTable->taper,
"TAPER_START",
sizeof(CHAR) * LIGOMETA_INSPIRALTAPER_MAX );
6872 if(taper!=
"") memcpy( simTable->taper, taper.Data(),
sizeof(CHAR) * LIGOMETA_INSPIRALTAPER_MAX );
6877 if(gps_start_time>0) {
6878 currentGpsTime+=timeStep;
6879 if((currentGpsTime>gpsEndTime) || (seed<0 && ninj>=sample.size())) isBreak=
true;
6881 if(ninj>=ninjections) isBreak=
true;
6886 clb << std::scientific;
6888 std::map<TString, int>::iterator iter;
6889 for (iter=hsample.begin(); iter!=hsample.end(); iter++) {
6890 if(
TString(iter->first).Contains(
"spcal")) clb << psample[iter->first] <<
"\t";
6892 simTable->bandpass = simulation_id;
6895 clb << simulation_id++;
6903 simTable = simTable->next = (SimInspiralTable *) calloc( 1,
sizeof(SimInspiralTable) );
6907 if(clb_file!=
"") clb.close();
6910 memset( &xmlfp, 0,
sizeof(LIGOLwXMLStream) );
6911 LAL_CALL( LALOpenLIGOLwXMLFile( &status, &xmlfp, xmlFile.Data() ), &status );
6912 if(injections.simInspiralTable) {
6913 LAL_CALL( LALBeginLIGOLwXMLTable( &status, &xmlfp, sim_inspiral_table ), &status );
6914 LAL_CALL( LALWriteLIGOLwXMLTable( &status, &xmlfp, injections, sim_inspiral_table ), &status );
6915 LAL_CALL( LALEndLIGOLwXMLTable ( &status, &xmlfp ), &status );
6917 LAL_CALL( LALCloseLIGOLwXMLFile ( &status, &xmlfp ), &status );
6918 LALCheckMemoryLeaks();
6919 cout << endl <<
"CWB::mdc::Posterior2XML - " 6920 <<
"Write\t" << ninj <<
"\tsamples to output xml file" << endl;
6926 std::map<TString, double>
6927 CWB::mdc::GetPsample(std::map<TString, int> hsample, vector<double> sample,
float if_lower,
float if_ref) {
6931 std::map<TString, double> psample;
6932 std::map<TString, double> psample_null;
6936 int pn_order = hsample[
"lal_pnorder"] ?
int(sample[hsample[
"lal_pnorder"]]) : -1;
6937 int amp_order = hsample[
"lal_amporder"] ?
int(sample[hsample[
"lal_amporder"]]) : -1;
6938 int approx = hsample[
"lal_approximant"] ?
int(sample[hsample[
"lal_approximant"]]) : IMRPhenomPv2;
6940 psample[
"pn_order"] = pn_order;
6941 psample[
"amp_order"] = amp_order;
6942 psample[
"approx"] = approx;
6946 double flow = hsample[
"flow"] ? sample[hsample[
"flow"]] : 30.;
6947 double f_ref = hsample[
"f_ref"] ? sample[hsample[
"f_ref"]] : 30.;
6948 #ifdef NEW_CODE_FOR_SEOBNRv4P_AND_SEOBNRv4PHM 6949 flow = 2*flow /(2+psample[
"amp_order"]);
6950 f_ref = 2*f_ref/(2+psample[
"amp_order"]);
6951 psample[
"amp_order"] = 0;
6953 if(if_lower>=0) flow = if_lower;
6954 if(if_ref>=0) f_ref = if_ref;
6956 psample[
"flow"] =
flow;
6957 psample[
"f_ref"] = f_ref;
6959 double lambda1 = hsample[
"lambda1"] ? sample[hsample[
"lambda1"]] : 0.;
6960 double lambda2 = hsample[
"lambda2"] ? sample[hsample[
"lambda2"]] : 0.;
6962 psample[
"lambda1"] = lambda1;
6963 psample[
"lambda2"] = lambda2;
6965 psample[
"logl"] = hsample[
"logl"] ? sample[hsample[
"logl"]] : 0.;
6966 if(hsample[
"logpost"]) {
6967 psample[
"post"] = sample[hsample[
"logpost"]];
6968 }
else if(hsample[
"logprior"] && hsample[
"logl"]) {
6969 psample[
"post"] = sample[hsample[
"logprior"]] + sample[hsample[
"logl"]];
6970 }
else if(hsample[
"post"]) {
6971 psample[
"post"] = sample[hsample[
"post"]];
6972 }
else psample[
"post"] = 0.;
6976 if(hsample[
"distance"]) {dist = sample[hsample[
"distance"]];dist_check++;}
6977 if(hsample[
"dist"]) {dist = sample[hsample[
"dist"]];dist_check++;}
6978 if(hsample[
"logdistance"]) {dist =
TMath::Exp(sample[hsample[
"logdistance"]]);dist_check++;}
6980 cout << endl <<
"CWB::mdc::Posterior2XML - " 6981 <<
"multiple definitions of distance in sample header " << endl;
6982 return psample_null;
6985 cout << endl <<
"CWB::mdc::Posterior2XML - " 6986 <<
"distance not defined in sample header " << endl;
6987 return psample_null;
6989 psample[
"dist"] =
dist;
6991 double q = sample[hsample[
"q"]];
6992 double factor = sample[hsample[
"mc"]] * pow(1 + q, 1.0/5.0);
6994 double m1 = factor * pow(q, -3.0/5.0);
6995 double m2 = factor * pow(q, 2.0/5.0);
6997 double q1 = hsample[
"q1"] ? sample[hsample[
"q1"]] : 0.;
6998 double q2 = hsample[
"q2"] ? sample[hsample[
"q2"]] : 0.;
7000 double ra = sample[hsample[
"ra"]];
7001 double dec = sample[hsample[
"dec"]];
7003 double psi = hsample[
"psi"] ? sample[hsample[
"psi"]] : sample[hsample[
"polarization"]];
7006 psample[
"dec"] = dec;
7007 psample[
"psi"] =
psi;
7014 psample[
"mchirp"] = pow(m1*m2,3./5.)/pow(m1+m2,1./5.);
7018 if(hsample[
"costheta_jn"]) {
7019 iota = TMath::ACos(sample[hsample[
"costheta_jn"]]);
7020 }
else if(hsample[
"theta_jn"]) {
7021 iota = sample[hsample[
"theta_jn"]];
7023 cout << endl <<
"CWB::mdc::Posterior2XML - " 7024 <<
"Error: costheta_jn or theta_jn not defined in posterior samples header file" << endl;
7025 return psample_null;
7029 static bool phase_warning=
true;
7030 if(hsample[
"phi_orb"]) {
7031 phase = sample[hsample[
"phi_orb"]];
7033 if(hsample[
"phase"]) {
7034 phase = sample[hsample[
"phase"]];
7036 phase = sample[hsample[
"phase_maxl"]];
7038 cout <<
"WARNING: Samples already marginalized over phase." << endl;
7039 cout <<
" Begrudgingly using max loglikelihood estimates." << endl;
7040 phase_warning=
false;
7044 psample[
"phase"] =
phase;
7046 double s1x=0, s1y=0, s1z = 0.;
7047 double s2x=0, s2y=0, s2z = 0.;
7050 double theta_jn, phi_jl, tilt1, tilt2, phi12, a1,
a2;
7051 if((hsample[
"theta_jn"]||hsample[
"costheta_jn"]) && hsample[
"phi_jl"]) {
7053 theta_jn = hsample[
"theta_jn"] ? sample[hsample[
"theta_jn"]] : TMath::ACos(sample[hsample[
"costheta_jn"]]);
7054 phi_jl = sample[hsample[
"phi_jl"]];
7056 if(hsample[
"tilt1"] && hsample[
"tilt2"]) {
7057 tilt1 = sample[hsample[
"tilt1"]];
7058 tilt2 = sample[hsample[
"tilt2"]];
7060 tilt1 = TMath::ACos(sample[hsample[
"costilt1"]]);
7061 tilt2 = TMath::ACos(sample[hsample[
"costilt2"]]);
7064 phi12 = sample[hsample[
"phi12"]];
7065 a1 = sample[hsample[
"a1"]];
7066 a2 = sample[hsample[
"a2"]];
7071 double phiRef=
phase;
7072 #if LAL_VERSION_MAJOR == 6 && LAL_VERSION_MINOR == 16 && LAL_VERSION_MICRO == 1 // used only for tagged version lalinference_o2 used for O2 posteriors 7073 int ret = XLALSimInspiralTransformPrecessingNewInitialConditions(&iota, &s1x, &s1y, &s1z, &s2x, &s2y, &s2z, theta_jn,
7074 phi_jl, tilt1, tilt2, phi12, a1, a2, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, f_ref);
7076 int ret = XLALSimInspiralTransformPrecessingNewInitialConditions(&iota, &s1x, &s1y, &s1z, &s2x, &s2y, &s2z, theta_jn,
7077 phi_jl, tilt1, tilt2, phi12, a1, a2, m1*LAL_MSUN_SI, m2*LAL_MSUN_SI, f_ref, phiRef);
7079 if( ret==XLAL_FAILURE ) {
7080 cout << endl <<
"CWB::mdc::Posterior2XML - " 7081 <<
"Error in XLALSimInspiralTransformPrecessingNewInitialConditions" << endl;
7082 return psample_null;
7088 if(hsample[
"a1z"] && hsample[
"a2z"]) {
7089 s1z = sample[hsample[
"a1z"]];
7090 s2z = sample[hsample[
"a2z"]];
7092 if(hsample[
"a1"] && hsample[
"a2"]) {
7093 s1z = sample[hsample[
"a1"]];
7094 s2z = sample[hsample[
"a2"]];
7096 cout << endl <<
"CWB::mdc::Posterior2XML - " 7097 <<
"Error: It must be non-spinning !!!" << endl;
7102 psample[
"iota"] = iota;
7104 psample[
"s1x"] = s1x;
7105 psample[
"s1y"] = s1y;
7106 psample[
"s1z"] = s1z;
7107 psample[
"s2x"] = s2x;
7108 psample[
"s2y"] = s2y;
7109 psample[
"s2z"] = s2z;
7112 static bool time_warning=
true;
7113 if(hsample[
"time"]) {
7114 time = sample[hsample[
"time"]];
7116 time = sample[hsample[
"time_maxl"]];
7118 cout <<
"WARNING: Samples already marginalized over time." << endl;
7119 cout <<
" Begrudgingly using max loglikelihood estimates." << endl;
7123 psample[
"time"] = time;
7125 if(hsample[
"h1_end_time"]) psample[
"h1_end_time"]=sample[hsample[
"h1_end_time"]];
else psample[
"h1_end_time"]=0.;
7126 if(hsample[
"l1_end_time"]) psample[
"l1_end_time"]=sample[hsample[
"l1_end_time"]];
else psample[
"l1_end_time"]=0.;
7127 if(hsample[
"v1_end_time"]) psample[
"v1_end_time"]=sample[hsample[
"v1_end_time"]];
else psample[
"v1_end_time"]=0.;
7128 if(hsample[
"g1_end_time"]) psample[
"g1_end_time"]=sample[hsample[
"g1_end_time"]];
else psample[
"g1_end_time"]=0.;
7131 std::map<TString, int>::iterator iter;
7132 for (iter=hsample.begin(); iter!=hsample.end(); iter++) {
7133 if(
TString(iter->first).Contains(
"spcal")) psample[iter->first] = sample[hsample[iter->first]];
7157 if(extension!=
".clb") {
7158 cout <<
"CWB::mdc::CalibrateInspiral - Error : Calibration file extension must be .clb" << endl;
7165 cout << endl <<
"CWB::mdc::CalibrateInspiral - " 7166 <<
"Error Opening Calibration File : " <<
inspCLB.Data() << endl;
7170 std::map<TString, int> hsample;
7174 in.getline(hline,4*1024);
7178 for(
int j=0;
j<token->GetEntries();
j++) {
7180 sheader.ReplaceAll(
" " ,
"");
7185 first_sheader=sheader;
7187 hsample[sheader] =
j;
7190 hsample[first_sheader] = token->GetEntries();
7194 int hsize=token->GetEntries();
7195 std::vector<double> sample(hsample.size());
7197 for(
int i=0;i<hsize;i++) in >> sample[
i];
7198 sample[hsize]=sample[0];
7199 if(!in.good())
break;
7200 if((
fabs(time-sample[hsample[
"time"]])<0.001)&&(simulation_id==sample[hsample[
"simulation_id"]])) {spcal=
true;
break;}
7205 printf(
"\nCWB::mdc::CalibrateInspiral - Calibration Found @ GPS time %10.6f\n",sample[hsample[
"time"]]);
7215 ifo_label.ToLower();
7216 if(sample[hsample[
"spcal_active"]]) {
7217 int spcal_npts = sample[hsample[
"spcal_npts"]];
7219 logfreq.
resize(spcal_npts);
7221 phase.
resize(spcal_npts);
7222 for(
int j=0;
j<spcal_npts;
j++) {
7223 logfreq[
j] =
log(sample[hsample[TString::Format(
"%s_spcal_freq_%d",ifo_label.Data(),
j)]]);
7224 amp[
j] = sample[hsample[TString::Format(
"%s_spcal_amp_%d",ifo_label.Data(),
j)]];
7225 phase[
j] = sample[hsample[TString::Format(
"%s_spcal_phase_%d",ifo_label.Data(),
j)]];
7234 TSpline3 spamp(
"spamp",logfreq.
data,amp.
data,logfreq.
size());
7235 TSpline3 spphs(
"spphs",logfreq.
data,phase.
data,logfreq.
size());
7242 double logF =
log(
i*df);
7243 if(logF<logfreq[0] || logF>logfreq[logfreq.
size()-1])
continue;
7244 double dA = spamp.Eval(logF);
7245 double dP = spphs.Eval(logF);
7248 TComplex cal_factor = (1.0+dA)*C.Exp(TComplex(0.,dP));
7250 TComplex W(w[2*
i],w[2*i+1]);
7272 cout <<
"CwB::mdc::GetXCorr error: start w1 != start w2" << endl;
7277 cout <<
"CwB::mdc::GetXCorr error: size w1 != size w2" << endl;
7289 complex<double>
A(W1[
i], W1[i+1]);
7290 complex<double>
B(W2[i],-W2[i+1]);
7291 complex<double>
C = A*
B;
7324 if(w1==NULL)
return w;
7325 if(w1->
size()==0)
return w;
7327 double R=w1->
rate();
7329 double b_w1 = w1->
start();
7331 double b_w2 = w2->
start();
7334 int o_w1 = b_w1>b_w2 ? 0 :
int((b_w2-b_w1)*R+0.5);
7335 int o_w2 = b_w1<b_w2 ? 0 :
int((b_w1-b_w2)*R+0.5);
7337 double start = b_w1>b_w2 ? b_w1 : b_w2;
7338 double stop = e_w1<e_w2 ? e_w1 : e_w2;
7339 int size =
int((stop-start)*R+0.5);
7393 double R=w1->
rate();
7395 double b_w1 = w1->
start();
7397 double b_w2 = w2->
start();
7400 int o_w1 = b_w1>b_w2 ? 0 :
int((b_w2-b_w1)*R+0.5);
7401 int o_w2 = b_w1<b_w2 ? 0 :
int((b_w1-b_w2)*R+0.5);
7403 double start = b_w1>b_w2 ? b_w1 : b_w2;
7404 double stop = e_w1<e_w2 ? e_w1 : e_w2;
7405 int size =
int((stop-start)*R+0.5);
7410 w.
start(b_w1+
double(o_w1)/R);
7465 for(
int i=0;
i<
size;
i++) {ET+=x[
i]*x[
i];}
7469 if((T>=start)&&(T<=stop)) {
7472 if(t<T) EL+=x[
i]*x[
i];
else ER+=x[
i]*x[
i];
7478 double el = EL*(1-P);
7479 double er = ER*(1-
Q);
7493 for(
int i=size-1;
i>0;
i--) {
7529 num+=w1[
i]*W2[
i]-W1[
i]*w2[
i];
7530 den+=w1[
i]*w2[
i]+W1[
i]*W2[
i];
7532 sync_phase = TMath::ATan2(num,den);
7537 double sync_xcor=0.0;
7538 for(
int i=0;
i<
size;
i++) sync_xcor+=w1[
i]*w2[
i];
7539 sync_xcor*=1./(w1.
rms()*sqrt(size));
7540 sync_xcor*=1./(w2.
rms()*sqrt(size));
7579 double sync_xcor=-1.e-20;
7580 double isync_xcor=0;
7581 for(
int i=0;i<size;i++) if(X[i]>sync_xcor) {isync_xcor=
i;sync_xcor=X[
i];}
7582 sync_xcor*=1./(w1.
rms()*sqrt(size));
7583 sync_xcor*=1./(w2.
rms()*sqrt(size));
7587 int isync_time = isync_xcor<size/2 ? -isync_xcor : size-isync_xcor;
7588 sync_time = isync_time/w1.
rate();
7635 for(
int i=0;
i<
size;
i++) phase[
i] = TMath::ATan2(w1W2[
i]-W1w2[
i],w1w2[i]+W1W2[i]);
7638 for(
int i=0;i<
size;i++) z1[i] = w1w2[i]*cos(phase[i])-W1w2[
i]*sin(phase[i]);
7640 double sync_xcor=-1.e-20;
7641 double isync_xcor=0;
7642 for(
int i=0;i<size;i++) if(z1[i]>sync_xcor) {isync_xcor=
i;sync_xcor=z1[
i];}
7644 int isync_time = isync_xcor<size/2 ? -isync_xcor : size-isync_xcor;
7645 sync_time = isync_time/w1.
rate();
7648 int isync_phase = isync_xcor;
7649 sync_phase = phase[isync_phase]*180./
TMath::Pi();
7653 for(
int i=0;i<
size;i++) sync_xcor+=w1[i]*w2[i];
7654 sync_xcor*=1./(w1.
rms()*sqrt(size));
7655 sync_xcor*=1./(w2.
rms()*sqrt(size));
7703 if(w1.
size()==0) {cout <<
"CWB::mdc::Align - Error: w1 size is 0" << endl;w1=
w;w2=
w;
return -1;}
7704 if(w2.
size()==0) {cout <<
"CWB::mdc::Align - Error: w2 size is 0" << endl;w1=
w;w2=
w;
return -1;}
7705 if(w1.
rate()!=w2.
rate()) {cout <<
"CWB::mdc::Align - Error: w1 rate != w2 rate" << endl;w1=
w;w2=
w;
return -1;}
7709 double b_w1 = w1.
start();
7711 double b_w2 = w2.
start();
7714 double b_w = b_w1<b_w2 ? b_w1 : b_w2;;
7715 double e_w = e_w1>e_w2 ? e_w1 : e_w2;;
7717 int o_w1 = b_w1<b_w2 ? 0 :
int((b_w1-b_w2)*R+0.5);
7718 int o_w2 = b_w2<b_w1 ? 0 :
int((b_w2-b_w1)*R+0.5);
7721 w.
resize(
int(R*((e_w-b_w)+0.5)));
7725 for(
int i=0;
i<w1.
size();
i++) w[
i+o_w1] = w1[
i];
7729 for(
int i=0;
i<w2.
size();
i++) w[
i+o_w2] = w2[
i];
7752 int vsize=w1.size();
7754 if(w2.size()!=vsize) {
7755 cout <<
"CWB::mdc::GetMatchFactor - Error: input vectors wavearray are inconsistents, different size" << endl;
7756 return std::numeric_limits<double>::max();
7759 if((tstart.size()>0) && (tstart.size()!=vsize)) {
7760 cout <<
"CWB::mdc::GetMatchFactor - Error: input vectors tstart are inconsistents, different size" << endl;
7761 return std::numeric_limits<double>::max();
7764 if((tstop.size()>0) && (tstop.size()!=vsize)) {
7765 cout <<
"CWB::mdc::GetMatchFactor - Error: input vectors tstop are inconsistents, different size" << endl;
7766 return std::numeric_limits<double>::max();
7770 cout <<
"CWB::mdc::GetMatchFactor - Error: input vectors are empty" << endl;
7771 return std::numeric_limits<double>::max();
7773 double rate=w1[0].rate();
7775 for(
int n=0;
n<vsize;
n++) {
7776 if((w1[
n].
rate()!=rate)||(w2[
n].
rate()!=rate)) {
7777 cout <<
"CWB::mdc::GetMatchFactor - Error: input wavearray are inconsistents, different rate" << endl;
7778 return std::numeric_limits<double>::max();
7781 double size=w1[0].size();
7783 for(
int n=0;
n<vsize;
n++) {
7784 if((w1[
n].
size()!=size)||(w2[
n].
size()!=size)) {
7785 cout <<
"CWB::mdc::GetMatchFactor - Error: input wavearray are inconsistents, different size" << endl;
7786 return std::numeric_limits<double>::max();
7789 double start=w1[0].start();
7791 for(
int n=0;
n<vsize;
n++) {
7793 cout <<
"CWB::mdc::GetMatchFactor - Error: input wavearray are inconsistents, different start time" << endl;
7794 return std::numeric_limits<double>::max();
7800 if(tstart.size()==0) {tstart.resize(vsize);
for(
int n=0;
n<vsize;
n++) tstart[
n]=start;}
7801 if(tstop.size()==0) {tstop.resize(vsize);
for(
int n=0;
n<vsize;
n++) tstop[
n] =start+size*dt;}
7807 for(
int n=0;
n<vsize;
n++) {
7810 if((t<=tstart[
n])||(t>tstop[n]))
continue;
7811 ew1 += w1[
n][
j]*w1[
n][
j];
7812 ew2 += w2[
n][
j]*w2[
n][
j];
7813 ew12 += w1[
n][
j]*w2[
n][
j];
7814 re += pow(w1[n][
j]-w2[n][
j],2);
7818 if((match==
"ff")&&(ew1*ew2==0))
return std::numeric_limits<double>::max();
7819 if((match==
"of")&&(ew2==0))
return std::numeric_limits<double>::max();
7821 if(match==
"ff")
return ew12/sqrt(ew1*ew2);
7822 if(match==
"of")
return ew12/sqrt(ew2*ew2);
7823 if(match==
"re")
return re;
7840 if(x==NULL)
return xq;
7841 if(x->
size()==0)
return xq;
7848 for(
int i=0;
i<x->
size();
i++) xq[
i] = sqrt(pow(x->
data[
i],2)+pow(xq[
i],2));
7868 if(x==NULL)
return xs;
7869 if(x->
size()==0)
return xs;
7875 double df=(double)xs.
rate()/(double)xs.
size();
7879 if(oneside)
for(
int i=0;
i<size/2;
i++) xs.
data[
i]*=sqrt(2.);
7904 double dF=(y.
rate()/(double)y.
size())/2.;
7906 for(
int j=0;
j<y.
size();
j+=2) {
7908 if(F<bF || F>eF) {y[
j]=0;y[
j+1]=0;}
7937 for(
int i=0;
i<
N;
i++) {ET+=x[
i]*x[
i];}
7939 double EF = ET*(1-P);
7944 for(
int i=0;
i<
N;
i++) {
7953 for(
int i=N-1;
i>0;
i--) {
7959 double dF=(x.
rate()/(double)x.
size())/2.;
7969 CWB::mdc::SimIMRSEOBNRv4ROMTimeOfFrequency(
double freq,
double m1,
double m2,
double chi1,
double chi2) {
7988 if(XLALSimIMRSEOBNRv4ROMTimeOfFrequency(&time, freq, m1*Mo, m2*Mo, chi1, chi2) == XLAL_FAILURE) {
7989 cout <<
"CWB::mdc::XLALSimIMRSEOBNRv4ROMTimeOfFrequency error" << endl;
7997 CWB::mdc::SimIMRSEOBNRv4ROMFrequencyOfTime(
double time,
double m1,
double m2,
double chi1,
double chi2) {
8017 if(XLALSimIMRSEOBNRv4ROMFrequencyOfTime(&freq, time, m1*Mo, m2*Mo, chi1, chi2) == XLAL_FAILURE) {
8018 cout <<
"CWB::mdc::XLALSimIMRSEOBNRv4ROMFrequencyOfTime error" << endl;
wavearray< double > t(hp.size())
std::vector< char * > ifoName
detector * getifo(size_t n)
param: detector index
void Dump(TString fname, int ID, int id, TString polarization)
detectorParams getDetectorParams()
waveform GetSourceWaveform(int &ID, int &id)
static double g(double e)
watplot * Draw(TString name, int id=0, TString polarization="hp", MDC_DRAW type=MDC_TIME, TString options="ALP", Color_t color=kBlack)
vector< mdcpar > sky_parms
size_t add(detector *)
param: detector structure return number of detectors in the network
static double GetTimeRange(wavearray< double > x, double &tMin, double &tMax, double efraction=EPZOOM)
std::vector< float > psiList
size_t readMDClog(char *, double=0., int=11, int=12)
param: MDC log file param: approximate gps time
TString Get(wavearray< double > &x, TString ifo)
double min(double x, double y)
std::vector< std::string > mdcList
virtual void rate(double r)
wavearray< double > a(hp.size())
static double GetCentralTime(wavearray< double > x)
char * watversion(char c='s')
std::vector< std::string > mdcType
void output(TTree *, network *, double, bool=true)
wavearray< double > GetSGQ(double frequency, double Q)
static void PhaseShift(wavearray< double > &x, double pShift=0.)
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1) *itime+ifreq;double time=itime *dT;double freq=(ifreq >0) ? ifreq *dF :dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
unsigned int srcList_seed
void GetSourceCoordinates(double &theta, double &phi, double &psi, double &rho, double &iota, double &hrss, int &ID, int &id)
std::vector< double > gpsList
Complex Exp(double phase)
static void AddExp(wavearray< double > &td, double v, int M)
wavearray< double > GetGA(double duration)
std::vector< float > phList
std::vector< std::string > mdcType
void Draw(int dpaletteId=1, Option_t *option="colfz")
std::vector< int > IDList
void setPolarization(POLARIZATION polarization=TENSOR)
virtual void start(double s)
std::vector< size_t > mdc__ID
void SetGridxColor(Color_t colorGridx=kBlack)
void plot(wavearray< double > &, char *=NULL, int=1, double=0., double=0., bool=false, float=0., float=0., bool=false, float=0., bool=false)
std::vector< detector * > ifoList
double GetDelay(TString ifo1, TString ifo2, double phi, double theta)
watplot * DrawTime(wavearray< double > &x, TString options="ALP", Color_t color=kBlack)
mdc & operator=(const mdc &)
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
mdcid AddWaveform(MDC_TYPE mdc_type, vector< mdcpar > par, TString uname="")
r add(tfmap, const_cast< char *>("hchannel"))
std::vector< std::string > nameList
network ** net
NOISE_MDC_SIMULATION.
void Draw(double t1=0.0, double t2=0.0, double f1=0.0, double f2=0.0, double z1=0.0, double z2=0.0, int dpaletteId=DUMMY_PALETTE_ID, Option_t *option="colfz")
void SetSkyDistribution(MDC_DISTRIBUTION sky_distribution, vector< mdcpar > par, int seed=0, bool add=false)
fprintf(stdout,"start=%f duration=%f rate=%f\, x.start(), x.size()/x.rate(), x.rate())
TString GetTemporaryFileName(TString tag="mdc", TString ext="txt", TString dir="/tmp", bool mkdir=false)
double GetPar(TString name, vector< mdcpar > par, bool &error)
int GetWaveformID(TString name)
POLARIZATION getPolarization()
virtual size_t size() const
void FillData(int size, double *phi, double *theta, double *binc)
std::vector< std::string > xmlType
void DrawSkyDistribution(TString name="skymap", TString projection="", TString coordinate="Geographic", double resolution=2, bool background=true)
wavearray< double > GetCGQ(double frequency, double Q)
static void AddGauss(wavearray< double > &td, double v, double u=0.)
wavecomplex antenna(double, double, double=0.)
param: source theta,phi, polarization angle psi in degrees
static double GetMatchFactor(TString match, vector< wavearray< double > > &w1, vector< wavearray< double > > &w2, vector< double > tstart=DEFAULT_VECtOR_DOUBLE, vector< double > tstop=DEFAULT_VECtOR_DOUBLE)
TString WriteFrameFile(TString frDir, TString frLabel, size_t gps, size_t length=1000, bool log=false, vector< TString > chName=vector< TString >())
static wavearray< double > GetEnvelope(wavearray< double > *x)
static wavearray< double > GetAdd(wavearray< double > *w1, wavearray< double > *w2)
void DumpLogHeader(TString fName, TString label="", int size=0)
static wavearray< double > GetAligned(wavearray< double > *w1, wavearray< double > *w2)
std::vector< std::string > mdcList
const char * DistributionToString(MDC_DISTRIBUTION n)
TString GetBurst(wavearray< double > &x, TString ifo)
static wavearray< double > GetXCorr(wavearray< double > &w1, wavearray< double > &w2)
printf("total live time: non-zero lags = %10.1f \, liveTot)
cout<< "Injected signals: "<< mdc.GetEntries()<< endl;cout<< "Injected signals in histogram factor_events_inj: "<< NEVTS<< endl;float myifar, ecor, m1, m2, netcc[3], neted, penalty;float rho[2];float chirp[6];float range[2];float frequency[2];float iSNR[3], sSNR[3];sim.SetBranchAddress("mass", mass);sim.SetBranchAddress("factor", &factor);sim.SetBranchAddress("range", range);sim.SetBranchAddress("chirp", chirp);sim.SetBranchAddress("rho", rho);sim.SetBranchAddress("netcc", netcc);sim.SetBranchAddress("neted", &neted);sim.SetBranchAddress("ifar", &myifar);sim.SetBranchAddress("ecor", &ecor);sim.SetBranchAddress("penalty", &penalty);sim.SetBranchAddress("time", mytime);sim.SetBranchAddress("iSNR", iSNR);sim.SetBranchAddress("sSNR", sSNR);sim.SetBranchAddress("spin", spin);sim.SetBranchAddress("frequency", frequency);float **volume=new float *[NBINS_mass1];float **volume_first_shell=new float *[NBINS_mass1];float **radius=new float *[NBINS_mass1];float **error_volume=new float *[NBINS_mass1];float **error_volume_first_shell=new float *[NBINS_mass1];float **error_radius=new float *[NBINS_mass1];for(int i=0;i< NBINS_mass1;i++) { volume[i]=new float[NBINS_mass2];volume_first_shell[i]=new float[NBINS_mass2];radius[i]=new float[NBINS_mass2];error_volume[i]=new float[NBINS_mass2];error_volume_first_shell[i]=new float[NBINS_mass2];error_radius[i]=new float[NBINS_mass2];for(int j=0;j< NBINS_mass2;j++) { volume[i][j]=0.;volume_first_shell[i][j]=0.;radius[i][j]=0.;error_volume[i][j]=0.;error_volume_first_shell[i][j]=0.;error_radius[i][j]=0.;} } float **spin_mtot_volume=new float *[NBINS_MTOT+1];float **spin_mtot_radius=new float *[NBINS_MTOT+1];float **error_spin_mtot_volume=new float *[NBINS_MTOT+1];float **error_spin_mtot_radius=new float *[NBINS_MTOT+1];for(int i=0;i< NBINS_MTOT+1;i++) { spin_mtot_volume[i]=new float[NBINS_SPIN+1];spin_mtot_radius[i]=new float[NBINS_SPIN+1];error_spin_mtot_volume[i]=new float[NBINS_SPIN+1];error_spin_mtot_radius[i]=new float[NBINS_SPIN+1];for(int j=0;j< NBINS_SPIN+1;j++) { spin_mtot_volume[i][j]=0.;error_spin_mtot_volume[i][j]=0.;spin_mtot_radius[i][j]=0.;error_spin_mtot_radius[i][j]=0.;} } char fname[1024];sprintf(fname, "%s/recovered_signals.txt", netdir);ofstream fev;fev.open(fname, std::ofstream::out);sprintf(line, "#GPS@L1 FAR[Hz] eFAR[Hz] Pval " "ePval factor rho frequency iSNR sSNR \");fev<< line<< endl;ofstream *fev_single=new ofstream[nfactor];for(int l=1;l< nfactor+1;l++) { sprintf(fname, "%s/recovered_signals_%d.txt", netdir, l);fev_single[l - 1].open(fname, std::ofstream::out);fev_single[l - 1]<< line<< endl;} double Vrho[RHO_NBINS], eVrho[RHO_NBINS], Rrho[RHO_NBINS], eRrho[RHO_NBINS], Trho[RHO_NBINS];for(int i=0;i< RHO_NBINS;i++) { Vrho[i]=0.;eVrho[i]=0.;Rrho[i]=0.;eRrho[i]=0.;Trho[i]=RHO_MIN+i *RHO_BIN;} double dV, dV1, dV_spin_mtot, nevts, internal_volume;int nT;int countv=0;int cnt=0;int cnt2=0;int cntfreq=0;bool bcut=false;double liveTot=sim.GetMaximum("ifar");double BKG_LIVETIME_yr=liveTot/CYS;double BKG_LIVETIME_Myr=BKG_LIVETIME_yr/(1.e6);cout.precision(14);cout<< "Total live time ---> background
wavearray< double > GetWNB(double frequency, double bandwidth, double duration, int seed=0, bool mode=0)
double RA2phi(double ph, double gps)
std::vector< int > idList
void SetGridyColor(Color_t colorGridy=kBlack)
static wavearray< double > GetBandpass(wavearray< double > x, double bF, double eF)
void DrawTF(wavearray< double > &x, TString options="")
static double TimeSync(wavearray< double > &w1, wavearray< double > &w2, double &sync_time)
watplot * DrawFFT(wavearray< double > &x, TString options="ALP", Color_t color=kBlack)
static double GetTimeBoundaries(wavearray< double > x, double P, double &bT, double &eT, double T=-1., double Q=-1.)
double phi2RA(double ph, double gps)
vector< waveform > wfList
static double cosi2e(double cosi)
static void AddCGBurst(wavearray< double > &td, double a, double f, double s, double d=0.)
TString GetBurstLog(source src, double FrameGPS, double SimHpHp, double SimHcHc, double SimHpHc)
double GravitationalConstant()
static void TimeShift(wavearray< double > &x, double tShift=0.)
MDC_COORDINATES mdc_coordinates
static void AddWGNoise(wavearray< double > &td, double a, double s)
std::vector< std::string > mdcName
static wavearray< double > GetDiff(wavearray< double > *w1, wavearray< double > *w2)
MDC_DISTRIBUTION sky_distribution
static double PhaseSync(wavearray< double > &w1, wavearray< double > &w2, double &sync_phase)
std::vector< float > rhoList
vector< source > GetSourceList(double start, double stop)
int getEBBH(double m1, double m2, double rmin0, double e0, wavearray< double > &Hp, wavearray< double > &Hx, double t_end)
static int Align(wavearray< double > &w1, wavearray< double > &w2)
double fabs(const Complex &x)
void resample(const wavearray< DataType_t > &, double, int=6)
waveform GetWaveform(int ID, int id=0)
void ReadWaveform(wavearray< double > &x, TString fName, double srate)
strcpy(RunLabel, RUN_LABEL)
static double e2cosi(double e)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
void GalacticToEquatorial(double ilongitude, double ilatitude, double &olongitude, double &olatitude)
double getTau(double, double)
param: source theta,phi angles in degrees
double distance_source_Kpc
TString GetParString(TString name, vector< mdcpar > par, bool &error)
std::vector< double > hrssList
std::vector< float > iotaList
static double TimePhaseSync(wavearray< double > &w1, wavearray< double > &w2, double &sync_time, double &sync_phase)
void DumpLog(TString fName, TString label="", bool append=false)
std::vector< source > srcList
static double GetCentralFrequency(wavearray< double > x)
#define MNGD_SOLAR_SISTEM_DISTANCE_FROM_GC
double GetAntennaPattern(TString ifo, double phi, double theta, double psi=0., TString polarization="hp")
void GeographicToCwb(double ilongitude, double ilatitude, double &olongitude, double &olatitude)
std::vector< double > mdcTime
wavearray< double > GetRD(double frequency, double tau, double iota, bool polarization=0)
static double GetFrequencyBoundaries(wavearray< double > x, double P, double &bF, double &eF)
for(int i=0;i< 101;++i) Cos2[2][i]=0
pointers to detectors
virtual void resize(unsigned int)
TString frLabel[NIFO_MAX]
void SetOptions(TString projection="hammer", TString coordinate="Geographic", double resolution=1, bool goff=false)
static wavearray< double > GetSpectrum(wavearray< double > *x, bool oneside=false)
static void AddSGBurst(wavearray< double > &td, double a, double f, double s, double d=0.)
std::vector< float > thList
void SetZaxisTitle(TString zAxisTitle)
MDC SetInspiral("inspNameTEST", inspOptions)
double SpeedOfLightInVacuo()