Logo coherent WaveBurst  
Library Reference Guide
Logo
GWOSC_Tools.C
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Gabriele Vedovato
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 // this macro is a collection of tools used manage the CWOSC data used for the cwb_gwosc command (see $CWB_SCRIPTS/Makefile.gwosc)
19 // these functions work for the GWOSC GWTC-1 events
20 
21 #include <lal/LALSimInspiral.h>
22 
23 #define MAX_IFOS 3
24 
25 void ConvertSamples(TString ifname, TString ofname, double gps);
26 void Posterior2XML(TString gwname, int seed=-1, TString model="");
27 void ConvertPSD(TString ifname);
28 void ConvertPSD2ASD(TString ifnamePSD, TString ofnameASD);
29 
30 void GWOSC_Tools(TString tool, TString spar1="", TString spar2="", double dpar1=0) {
31 
32  if(tool=="ConvertSamples") ConvertSamples(spar1, spar2, dpar1);
33  if(tool=="Posterior2XML") Posterior2XML(spar1,int(dpar1),spar2);
34  if(tool=="ConvertPSD") ConvertPSD(spar1);
35  if(tool=="ConvertPSD2ASD") {
36  vector<TString> fileList = CWB::Toolbox::getFileListFromDir(".",".dat",spar1,"",true);
37  for(int i=0;i<fileList.size();i++) {
38  TString ifnamePSD = fileList[i];
39  TString ofnameASD = fileList[i]; ofnameASD.ReplaceAll("PSD","ASD");
40  ConvertPSD2ASD(ifnamePSD,ofnameASD);
41  }
42  }
43 }
44 
45 void ConvertSamples(TString ifname, TString ofname, double gps) {
46 //
47 // this macro is used to convert public samples (IMRPhenomPv2) format to the format used by cWB
48 //
49 // input: ifname: name of input posterior sample file, obtained from the public hdf5 file and converted with h5dump command
50 // eg: h5dump -o GW150914_GWTC-1.txt -y -w 400 GW150914_GWTC-1.hdf5
51 // gps: gps time used to generate the xml file (the gps time is not included in the public hdf5 file)
52 // output: ofname: name of the converted output file
53 //
54 
55  std::map<TString, double> psample;
56  std::vector<TString> pname;
57 
58  TString name="";
59 
60  name="dummy"; psample[name] = 0; pname.push_back(name); // defined to avoid a bug in mdc::Posterior2XML
61 
62  name="time"; psample[name] = gps; pname.push_back(name); // coa time
63 
64  name="costheta_jn"; psample[name] = 0; pname.push_back(name); // simTable->inclination
65  name="distance"; psample[name] = 0; pname.push_back(name); // luminosity distance
66  name="ra"; psample[name] = 0; pname.push_back(name);
67  name="dec"; psample[name] = 0; pname.push_back(name);
68  name="m1"; psample[name] = 0; pname.push_back(name); // mass1 in the detector frame
69  name="m2"; psample[name] = 0; pname.push_back(name); // mass2 in the detector frame
70  name="a1"; psample[name] = 0; pname.push_back(name); // spin1
71  name="a2"; psample[name] = 0; pname.push_back(name); // spin2
72  name="costilt1"; psample[name] = 0; pname.push_back(name); // used to compute spin components s1x,s1y,s1z
73  name="costilt2"; psample[name] = 0; pname.push_back(name); // used to compute spin components s2x,s2y,s2z
74 
75  if(ifname.Contains("GW170817")) {
76  name="lambda1"; psample[name] = 0; pname.push_back(name); // lambda1
77  name="lambda2"; psample[name] = 0; pname.push_back(name); // lambda2
78  }
79 
80  name="q"; psample[name] = 0; pname.push_back(name); // mass ratio
81  name="mc"; psample[name] = 0; pname.push_back(name); // chirp mass
82 
83  name="lal_pnorder"; psample[name] = 8; pname.push_back(name); // pseudoFourPN
84  name="lal_amporder"; psample[name] = 0; pname.push_back(name);
85  name="lal_approximant"; psample[name] = IMRPhenomPv2; pname.push_back(name); // IMRPhenomPv2
86  name="flow"; psample[name] = 20; pname.push_back(name); // simTable->f_lower
87  name="f_ref"; psample[name] = 20; pname.push_back(name); // simTable->f_final
88 
89  name="phase"; psample[name] = 0; pname.push_back(name); // simTable->coa_phase
90  name="psi"; psample[name] = 0; pname.push_back(name); // polarization
91  name="phi12"; psample[name] = 0; pname.push_back(name); // used to compute spin components s1x,s1y,s1z
92  name="phi_jl"; psample[name] = 0; pname.push_back(name); // used to compute spin components s2x,s2y,s2z
93 
94 
95  ifstream in;
96  in.open(ifname.Data(),ios::in);
97  if(!in.good()) {cout << "Error Opening File : " << ifname.Data() << endl;exit(1);}
98 
99  ofstream out;
100  out.open(ofname.Data(),ios::out);
101  if(!out.good()) {cout << "Error Opening File : " << ofname.Data() << endl;exit(1);}
102  out.precision(14);
103 
104  // write header
105  for(int i = 0; i < pname.size(); i++) out << pname[i] << "\t";
106  out << endl;
107 
108  int entries=0;
109  char line[1024];
110  bool data=false;
111  while (1) {
112  in >> line; if(!in.good()) break;
113  if(TString(line).Contains("{")) {
114  in >> line; if(!in.good()) break;
115  data=true;
116  }
117  if(data) {
118  psample["costheta_jn"] = TString(line).Atof(); in >> line; if(!in.good()) break;
119  psample["distance"] = TString(line).Atof(); in >> line; if(!in.good()) break;
120  psample["ra"] = TString(line).Atof(); in >> line; if(!in.good()) break;
121  psample["dec"] = TString(line).Atof(); in >> line; if(!in.good()) break;
122  psample["m1"] = TString(line).Atof(); in >> line; if(!in.good()) break;
123  psample["m2"] = TString(line).Atof(); in >> line; if(!in.good()) break;
124  if(ifname.Contains("GW170817")) {
125  psample["lambda1"] = TString(line).Atof(); in >> line; if(!in.good()) break;
126  psample["lambda2"] = TString(line).Atof(); in >> line; if(!in.good()) break;
127  }
128  psample["a1"] = TString(line).Atof(); in >> line; if(!in.good()) break;
129  psample["a2"] = TString(line).Atof(); in >> line; if(!in.good()) break;
130  psample["costilt1"] = TString(line).Atof(); in >> line; if(!in.good()) break;
131  psample["costilt2"] = TString(line).Atof(); in >> line; if(!in.good()) break;
132 
133  double m1 = psample["m1"];
134  double m2 = psample["m2"];
135  psample["mc"] = pow(m1*m2,3./5.)/pow(m1+m2,1./5.);
136  psample["q"] = m2<m1 ? m2/m1 : m1/m2;
137 
138  // write data
139  for(int i = 0; i < pname.size(); i++) out << psample[pname[i]] << "\t";
140  out << endl;
141 
142  entries++;
143 /*
144  cout << "----------------------------------------------------------" << endl;
145  cout << entries << endl;
146  cout << "----------------------------------------------------------" << endl;
147 
148  for (int i = 0; i < pname.size(); i++) {
149  cout << pname[i] << "\t" << psample[pname[i]] << endl;
150  }
151 */
152  }
153  if(TString(line).Contains("}")) data=false;
154 
155 // if(entries>3) break;
156  }
157 
158  cout << endl << "ConvertSamples.C - converted entries: " << entries << endl;
159 
160  cout << endl << "ConvertSamples.C - created file: " << ofname << endl << endl;
161 
162  out.close();
163  in.close();
164 
165  exit(0);
166 }
167 
168 void Posterior2XML(TString gwname, int seed, TString model) {
169 //
170 // extract one sample from the posterior samples file and convert into the xml format
171 // the public files are at: https://dcc.ligo.org/LIGO-P1800370/public
172 // how to download, eg: wget https://dcc.ligo.org/public/0157/P1800370/005/GW150914_GWTC-1.hdf5 -q -O GW150914.hdf5
173 //
174 // input: gwname: GW name
175 // seed: seed used to randomly select a sample from the gwname+"_PE.dat" posterior sample file
176 // if =0 then the maxl sample is selected
177 // if =1 then the map sample is selected
178 // if <0 then the first sample is selected
179 //
180 // model: PE approximant
181 //
182 
183  gRandom->SetSeed(seed);
184 
185  TString options = "";
186 
187  options += "--ninjections 1 ";
188 
189  if(seed==0) options += "--sample maxl ";
190  else if(seed==1) options += "--sample map ";
191  else options += TString::Format("--seed %d ",seed);
192 
193  TString file_tag = model!="" ? gwname+"_PE_"+model : gwname+"_PE";
194  TString approximant = model; approximant.ReplaceAll("C01:","");
195 
196  if(model=="") approximant = "IMRPhenomPv2"; // used for GWTC-1 catalog
197 
198  options += "--source "+gwname+" ";
199  options += "--waveform "+approximant+"pseudoFourPN ";
200  options += "--f_ref 20 ";
201  options += "--f_lower 20 ";
202 
203  TString posteriorFile = file_tag+".dat";
204 
205  TString xmlFile = file_tag+".xml";
206  options += "--clb_file "+file_tag+".clb ";
207 
208  CWB::mdc::Posterior2XML(posteriorFile, xmlFile, options);
209 
210  exit(0);
211 }
212 
213 void ConvertPSD(TString ifname) {
214 //
215 // convert the GWOSC PSD curves to the format used by the the CWB_Plugin_SimNoise.C plugin
216 // the public files are at: https://dcc.ligo.org/LIGO-P1900011/public
217 // how to download, eg: wget https://dcc.ligo.org/public/0158/P1900011/001/GWTC1_GW150914_PSDs.dat -q -O GWTC1_GW150914_PSDs.dat
218 //
219 // input: ifname: input public PSD file
220 //
221 
222  ifstream in;
223  in.open(ifname.Data(), ios::in);
224  if (!in.good()) {cout << "Error Opening File : " << ifname << endl;exit(1);}
225 
226  int size=0;
227  char str[1024];
228  int fpos=0;
229  while(true) {
230  in.getline(str,1024);
231  if (!in.good()) break;
232  if(str[0] != '#') size++;
233  }
234  cout << "size " << size << endl;
235  in.clear(ios::goodbit);
236  in.seekg(0, ios::beg);
237 
238  in.getline(str,1024); // get header
239  TString header = str;
240 
241  int iH1 = header.Index("LIGO_Hanford_PSD");
242  int iL1 = header.Index("LIGO_Livingston_PSD");
243  int iV1 = header.Index("Virgo_PSD");
244 
245  double iIFO[MAX_IFOS];
246  TString sIFO[MAX_IFOS];
247 
248  int nIFO=0;
249  if(iH1>=0) {iIFO[nIFO]=iH1;sIFO[nIFO]="H1";nIFO++;}
250  if(iL1>=0) {iIFO[nIFO]=iL1;sIFO[nIFO]="L1";nIFO++;}
251  if(iV1>=0) {iIFO[nIFO]=iV1;sIFO[nIFO]="V1";nIFO++;}
252 
253  if(nIFO==0 || nIFO>MAX_IFOS) {cout << "Error: no IFOs in PSDs file or not allowed Format" << endl; exit(1);}
254 
255  double* freq = new double[size];
256  double* psd[MAX_IFOS]; for(int n=0;n<MAX_IFOS;n++) psd[n] = new double[size];
257 
258  int lines=0;
259  while(1) {
260  if(nIFO==1) in >> freq[lines] >> psd[0][lines];
261  if(nIFO==2) in >> freq[lines] >> psd[0][lines] >> psd[1][lines];
262  if(nIFO==3) in >> freq[lines] >> psd[0][lines] >> psd[1][lines] >> psd[2][lines];
263  if (!in.good()) break;
264  size=lines;
265  lines++;
266  }
267  in.close();
268 
269  int idx[3];
270  TMath::Sort(nIFO,iIFO,idx,false);
271 
272  for(int n=0;n<nIFO;n++) {
273  int m=idx[n];
274  TString ofname = ifname;
275  ofname.ReplaceAll("_PSDs.dat","_ASD_"+sIFO[m]+".dat");
276  cout << ofname << endl;
277  ofstream out;
278  out.open(ofname.Data(),ios::out);
279  if(!out.good()) {cout << "Error Opening Output File : " << ofname << endl;exit(1);}
280  out << 0 << "\t" << 10*sqrt(psd[m][0]) << endl; // blind the detector at low frequencies
281  for(int i=0;i<size;i++) out << freq[i] << "\t" << sqrt(psd[m][i]) << endl;
282  out << 8192 << "\t" << 10*sqrt(psd[m][size-1]) << endl; // blind the detector at high frequencies
283  out.close();
284  }
285 
286  exit(0);
287 }
288 
289 void ConvertPSD2ASD(TString ifnamePSD, TString ofnameASD) {
290 
291  #define SH_SUP 1e-40
292 
293  ifstream in;
294  in.open(ifnamePSD.Data(),ios::in);
295  if (!in.good())
296  {cout << "ConvertPSD2ASD - Error Opening File : "
297  << ifnamePSD.Data() << endl;gSystem->Exit(1);}
298 
299  cout << "ConvertPSD2ASD - Read File : " << ifnamePSD.Data() << endl;
300 
301  vector<double> vfr;
302  vector<double> vsh;
303  char sfr[1024],ssh[1024];
304  char str[1024];
305  double fr,sh;
306 
307  int fpos;
308  while (1) {
309  fpos=in.tellg();
310  in.getline(str,1024);
311  if(str[0] == '#') continue;
312  in.seekg(fpos, ios::beg);
313 
314  in >> sfr >> ssh;
315  if (!in.good()) break;
316  fpos=in.tellg();
317  in.seekg(fpos+1, ios::beg);
318  fr = TString(sfr).Atof();
319  sh = TString(ssh).Atof();
320  if(sh<=0)
321  {cout << "ConvertPSD2ASD - input sensitivity file : " << ifnamePSD.Data()
322  << " contains zero at frequency : " << fr << " Hz " << endl;gSystem->Exit(1);}
323  vfr.push_back(fr); vsh.push_back(sh);
324  }
325  in.close();
326 
327  // set high sh values at beginning and end frequency range to exclude that frequency from the analysis
328  vsh[0] = vsh[1];
329  fr = 2*vfr[vfr.size()-1] - vfr[vfr.size()-2];
330  sh = vsh[vsh.size()-1];
331  vfr.push_back(fr); vsh.push_back(sh);
332  vfr.push_back(8192.); vsh.push_back(sh);
333 
334  for(int i=0;i<vsh.size();i++) if(vsh[i]>SH_SUP) vsh[i]=SH_SUP;
335 
336  // overwrite input psd with the PSD values ASD=sqrt(PSD)
337  cout << "ConvertPSD2ASD - Write File : " << ofnameASD.Data() << endl;
338  ofstream out;
339  out.open(ofnameASD.Data(),ios::out);
340  for(int i=0;i<vfr.size();i++) out << vfr[i] << "\t" << sqrt(vsh[i]) << endl; // oneside spectrum
341  out.close();
342 
343  return;
344 }
345 
static vector< TString > getFileListFromDir(TString dir_name, TString endString="", TString beginString="", TString containString="", bool fast=false)
Definition: Toolbox.cc:5108
double m1
par [0] name
int n
Definition: cwb_net.C:28
TString("c")
TString posteriorFile
Definition: Posterior2XML.C:33
ofstream out
Definition: cwb_merge.C:214
wavearray< double > psd(33)
CWB::frame fr(FRLIST_NAME)
TString xmlFile
Definition: Posterior2XML.C:34
Long_t size
int m
Definition: cwb_net.C:28
i drho i
model
Definition: h52psd.py:25
#define nIFO
char str[1024]
wavearray< double > freq
Definition: Regression_H1.C:79
void ConvertSamples(TString ifname, TString ofname, double gps)
Definition: GWOSC_Tools.C:37
void ConvertPSD2ASD(TString ifnamePSD, TString ofnameASD)
Definition: GWOSC_Tools.C:289
#define SH_SUP
void GWOSC_Tools(TString tool, TString spar1="", TString spar2="", double dpar1=0)
Definition: GWOSC_Tools.C:29
char options[256]
#define MAX_IFOS
Definition: GWOSC_Tools.C:23
double gps
ifstream in
cout<< "Starting reading output directory ..."<< endl;vector< TString > fileList
double m2
void Posterior2XML(TString gwname, int seed=-1)
Definition: GWOSC_Tools.C:160
char line[1024]
void ConvertPSD(TString ifname)
Definition: GWOSC_Tools.C:193
exit(0)