Skip to content

Commit 4e2d2ee

Browse files
author
Ilya Mandel
committed
First draft of detection code
1 parent e991085 commit 4e2d2ee

File tree

2 files changed

+270
-22
lines changed

2 files changed

+270
-22
lines changed

‎compas_matlab_utils/CosmicHistoryIntegrator.m‎

Lines changed: 111 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -1,28 +1,44 @@
1-
function [Zlist, MergerRateByRedshiftByZ, SFR, Zweight]=...
2-
CosmicHistoryIntegrator(filename, zlistformation, zlistdetection, Msimulated, makeplots)
1+
function [Zlist, MergerRateByRedshiftByZ, SFR, Rdetections, DetectableMergerRate, Mtzlist, etalist, zlistdetection]=...
2+
CosmicHistoryIntegrator(filename, zlistformation, zmaxdetection, Msimulated, makeplots)
33
% Integrator for the binary black hole merger rate over cosmic history
44
% COMPAS (Compact Object Mergers: Population Astrophysics and Statistics)
55
% software package
66
%
77
% USAGE:
8-
% [Zlist, MergerRateByRedshiftByZ]=...
9-
% CosmicHistoryIntegrator(filename, zlistformation, zlistmerger, [,makeplots, filename2, name1, name2])
8+
% [Zlist, MergerRateByRedshiftByZ, SFR, Zweight, Rdetections, DetectableMergerRate, Mtzlist, etalist, zlistdetection]=...
9+
% CosmicHistoryIntegrator(filename, zlistformation, zmaxdetection, Msimulated, [,makeplots])
1010
%
1111
% INPUTS:
1212
% filename: name of population synthesis input file
1313
% should be in COMPAS output h5 format
1414
% zlistformation: vector of redshifts at which the formation rate is
1515
% computed
16-
% zlistdetection: vector of redshifts at which the detection rate is computed
16+
% zmaxdetection: maximum redshift to which the detection rate is computed
1717
% Msimulated: total star forming mass represented by the simulation (for
1818
% normalisation)
1919
% makeplots: if set to 1, generates a set of useful plots (default = 0)
2020
%
21-
% Zlist is a vector of metallicities, taken from input file
22-
% MergerRateByRedshiftByZ is a matrix of size length(Zlist) X length(zlist)
23-
% which contains a merger rate of binary black holes in the given redshift
24-
% and metallicity bin, in units of mergers per Gpc^3 of comoving volume per
21+
% OUTPUTS:
22+
% Zlist is a vector of metallicities, taken from the COMPAS run input file
23+
% MergerRateByRedshiftByZ is a matrix of size length(Zlist) X length(zformationlist)
24+
% which contains a merger rate of merging compact objects in the given redshift
25+
% and metallicity bin, in units of mergers per Mpc^3 of comoving volume per
2526
% year of source time
27+
% SFR is a vector of size length(zlistformation) containing the star formation rate
28+
% (solar masses per Mpc^3 of comoving volume per year of source time)
29+
% Rdetection is a matrix of size length(zlistdetection) X length Mtzlist X
30+
% length(etalist) containing the detection rate per year of observer time
31+
% from a given redshift bin and redshifted total mass and symmetric mass
32+
% ratio pixel
33+
% DetectableMergerRate is a matrix of the same size as Rdetection but
34+
% containing the intrinsic rate of detectable mergers per Mpc^3 of comoving
35+
% volume per year of source time
36+
% Mtzlist is a list of redshifted total mass bins (for computational
37+
% efficiency, very rare sources with Mtz>200 Msun are folded into the
38+
% last bin)
39+
% etalist is a list of symmetric mass ratio bins
40+
% zlistdetection is a vector of redshifts at which detection rates are
41+
% computed (a subset of zlistformation going up to zmaxdetection)
2642
%
2743
% EXAMPLE:
2844
% CosmicHistoryIntegrator('~/Work/COMPASresults/runs/Zdistalpha1-031803.h5', ...
@@ -53,7 +69,7 @@
5369
global yr;
5470
Mpcm=1*10^6 * 3.0856775807e16; %Mpc in meters
5571
c=299792458; %speed of light, m/s
56-
Mpc=Mpcm/c; %Gpc in seconds
72+
Mpc=Mpcm/c; %Mpc in seconds
5773
yr=3.15569e7; %year in seconds
5874

5975
if (nargin<4)
@@ -62,9 +78,9 @@
6278
if (nargin<5), makeplots=0; end;
6379

6480
%cosmology calculator
65-
[tL]=Cosmology(zlistformation);
81+
[tL,Dl,dVc]=Cosmology(zlistformation);
6682
%load COMPAS data
67-
[M1,M2,Z,Tdelay]=DataRead(filename);
83+
[M1,M2,Z,Tdelay,maxNS]=DataRead(filename);
6884
%metallicity-specific SFR
6985
[SFR,Zlist,Zweight]=Metallicity(zlistformation,min(Z),max(Z));
7086

@@ -75,21 +91,37 @@
7591
%metallicity-specific star formation rate
7692
dz=zlistformation(2)-zlistformation(1);
7793
tLmerge=tL(floor(zlistformation/dz)+1);
94+
etavec=0.01:0.01:0.25;
95+
Mtzvec=1:1:200; %ignore things on wrong side of mass gap, go up to z=1
7896
MergerRateByRedshiftByZ=zeros(length(zlistformation),length(Zlist));
97+
MergerRateByRedshiftByMtzByEta=zeros(length(zlistformation),length(etavec),length(Mtzvec));
7998
for(i=1:length(M1)),
99+
i, M1(i), M2(i)
80100
Zcounter=find(Zlist>=Z(i),1);
81101
zformindex=find(tL>=Tdelay(i),1);
102+
etaindex=M1(i)^0.6*M2(i)^0.6/(M1(i)+M2(i))^0.2;
82103
for(k=1:length(zlistformation)),
83-
zformindex=find(tL>=(Tdelay(i)+tLmerge(k)),1);
104+
zformindex=find(tL>=(Tdelay(i)+tLmerge(k)),1)
105+
Mtzindex=ceil((M1(i)+M2(i))*zlistformation(k))
84106
MergerRateByRedshiftByZ(k,Zcounter)=...
85107
MergerRateByRedshiftByZ(k,Zcounter)+...
86-
SFR(zformindex)*Zweight(zformindex,Zcounter)/Msimulated;
108+
SFR(zformindex)*Zweight(zformindex,Zcounter)/Msimulated;
109+
MergerRateByRedshiftByMtzByEta(k,Mtzindex,etaindex) =...
110+
MergerRateByRedshiftByMtzByEta(k,Mtzindex,etaindex) + ...
111+
SFR(zformindex)*Zweight(zformindex,Zcounter)/Msimulated;
87112
end;
88113
end;
89114

115+
zlistdetection=zlistformation(1:find(zlistformation<=zmaxdetection,1,"last"));
116+
fin=load('~/Work/Rai/LIGOfuture_data/freqVector.txt');
117+
%noise=load('~/Work/Rai/LIGOfuture_data/dataNomaLIGO.txt');
118+
noise=load('~/Work/Rai/LIGOfuture_data/dataEarly_low.txt');
119+
[Rdetections,DetectableMergerRate]=...
120+
DetectionRate(zlistformation,Mtzlist,etalist,MergerRateByRedshiftByMtzByEta,zlistdetection,fin,noise,Dl,dVc)
121+
90122
if(makeplots==1), %make a set of default plots
91123
MakePlots(M1,M2,Z,Tdelay,zlistformation,Zlist,SFR,Zweight,...
92-
MergerRateByRedshiftByZ, 1);
124+
MergerRateByRedshiftByZ, Rdetections, DetectableMergerRate, Mtzlist, etalist, 1);
93125
end;
94126

95127
end %end of CosmicHistoryIntegrator
@@ -98,7 +130,7 @@
98130
%Load the data stored in COMPAS .h5 output format from a file
99131
%Select only double compact object mergers of interest, and return the
100132
%component masses, metallicities, and star formation to merger delay times
101-
function [M1,M2,Z,Tdelay]=DataRead(file)
133+
function [M1,M2,Z,Tdelay, maxNS]=DataRead(file)
102134
if(exist(file, 'file')~=2),
103135
error('Input file does not exist');
104136
end;
@@ -119,6 +151,7 @@
119151
%NSBH=(((type1==13) & (type2==14)) | ((type1==14) & (type2==13)));
120152
%mergingDCO=mergingBNS | mergingNSBH | mergingBBH;
121153
%BNScount=sum(mergingBNS); NSBHcount=sum(mergingNSBH); BBHcount=sum(mergingBBH);
154+
maxNS=max(max(mass1(type1==13)), max(mass2(type2==13)));
122155
chirpmass=mass1.^0.6.*mass2.^0.6./(mass1+mass2).^0.2;
123156
q=mass2./mass1;
124157
seedCE=h5read(file,'/BSE_Common_Envelopes/SEED');
@@ -135,9 +168,9 @@
135168
M1=mass1(mergingDCO); M2=mass2(mergingDCO); Z=Zdco(mergingDCO); Tdelay=Ttotal(mergingDCO);
136169
end %end of DataRead
137170

138-
%Compute the star formation rate and lookback time (in years)
139-
%for an array of redshifts
140-
function [tL]=Cosmology(zvec)
171+
%Compute the lookback time (yr), lumionosity distance (Mpc), and comoving
172+
%volume (Mpc^3) for an array of redshifts
173+
function [tL, Dl, dVc]=Cosmology(zvec)
141174
global Mpcm
142175
global Mpc
143176
global yr
@@ -152,7 +185,7 @@
152185
E=sqrt(OmegaM.*(1+zvec).^3+OmegaL); %Hogg, astro-ph/9905116, Eq. 14
153186
Dc=Dh*dz*cumsum(1./E); %Hogg, Eq. 15
154187
Dm=Dc; %Hogg, Eq. 16, k=0;
155-
Dl=(1+zvec).*Dm; %Hogg, Eq. 20
188+
Dl=(1+zvec).*Dm/Mpc; %Hogg, Eq. 20
156189
%see also Eq. (1.5.46) in Weinberg, "Cosmology", 2008
157190
dVc=4*pi*Dh^3*(OmegaM*(1+zvec).^3+OmegaL).^(-0.5).*(Dc/Dh).^2*dz/Mpc^3;
158191
Vc=cumsum(dVc);
@@ -193,9 +226,65 @@
193226
end %end of Metallicity
194227

195228

229+
%Compute detection rates per unit observer time and per unit source time
230+
%per unit comoving volume as a function of redshifted total mass and eta
231+
function [Rdetections,DetectableMergerRate]=...
232+
DetectionRate(zlistformation,Mtzlist,etalist,MergerRateByRedshiftByMtzByEta,zlistdetection,freqfile,noisefile,Dl,dVc)
233+
fin=load('~/Work/Rai/LIGOfuture_data/freqVector.txt');
234+
noise=load('~/Work/Rai/LIGOfuture_data/dataEarly_low.txt');
235+
236+
flow=10;
237+
df=1;
238+
f=flow:df:500; %BBH focussed
239+
Sf=interp1(fin, noise.^2, f);
240+
241+
Ntheta=1e6;
242+
psi=rand(1,Ntheta)*pi;
243+
phi=rand(1,Ntheta)*2*pi;
244+
costh=rand(1,Ntheta);
245+
cosiota=rand(1,Ntheta);
246+
sinth=sqrt(1-costh.^2);
247+
siniota=sqrt(1-cosiota.^2);
248+
Fplus=1/2*(1+costh.^2).*cos(2*phi).*cos(2*psi)-costh.*sin(2*phi).*sin(2*psi);
249+
Fcross=1/2*(1+costh.^2).*cos(2*phi).*sin(2*psi)+costh.*sin(2*phi).*cos(2*psi);
250+
Theta=1/2*sqrt(Fplus.^2.*(1+cosiota.^2).^2+4*Fcross.^2.*cosiota.^2);
251+
Thetas=sort(Theta);
252+
253+
254+
SNRat1Mpc=zeros(length(Mtzvec),length(etavec));
255+
for(i=1:length(Mtzvec)),
256+
for(j=1:length(etavec)),
257+
[h,Am,psi]=IMRSAWaveform(f, Mtzvec(i), etavec(j), 0, 0, 0, 1, flow);
258+
integral=sum(4*Am.^2./Sf*df);
259+
SNRat1Mpc(i,j)=sqrt(integral);
260+
end;
261+
end;
262+
263+
SNR=zeros(length(zlistdetection),length(Mtzlist),length(etalist));
264+
for(i=1:length(zlistdetection)), SNR(i,:,:)=1./(Dl(i)./Mpc); end;
265+
SNR=SNR.*SNRat1Mpc;
266+
267+
268+
SNR8pre=0.1:0.1:1000;
269+
theta=1./SNR8pre;
270+
pdetect=1-interp1([0,Thetas,1],[(0:Ntheta)/Ntheta,1],theta);
271+
272+
Rdetections=zeros(length(zlistdetection),length(Mtzlist),length(etalist)); %Detections per unit observer time
273+
DetectableMergerRate=zeros(length(zlistdetection),length(Mtzlist),length(etalist)); %Detections per unit source time per unit Vc
274+
SNR8=SNR/8;
275+
pdetection=zeros(size(Rdetections));
276+
pdetection(SNR8>1)=pdetect(floor(SNR8(SNR8>1 & SNR8<max(SNR8pre))*10));
277+
pdetection(SNR8>max(SNR8pre))=1;
278+
DetectableMergerRate=MergerRateByRedshiftByMtzByEta(i,1:length(zlistdetection)).*pdetection;
279+
Rdetections(i,:)=MergerRateByRedshiftByMtzByEta(i,1:length(zlistdetection)).*pdetection.*dVc(1:maxzdetindex)./(1+zlistdetection);
280+
281+
end %end of DetectionRate
282+
196283
%Make a set of default plots
197-
function MakePlots(M1,M2,Z,Tdelay,zvec,Zlist,SFR,Zweight,...
198-
MergerRateByRedshiftByZ, fignumber)
284+
function MakePlots(M1,M2,Z,Tdelay,zformationlist,Zlist,SFR,Zweight,...
285+
MergerRateByRedshiftByZ, Rdetections, DetectableMergerRate, zdetectionlist, Mtzlist, etazlist, fignumber)
286+
287+
zvec=zformationlist;
199288

200289
figure(fignumber), clf(fignumber); %,colormap jet;
201290
plot(zvec, sum(MergerRateByRedshiftByZ,2)*1e9, 'LineWidth', 3), hold on;

‎compas_matlab_utils/IMRSAWaveform.m‎

Lines changed: 159 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,159 @@
1+
function [h, Aeff, psiEff] = IMRSAWaveform (fVec, totalMass, eta, chi, startPhase, ...
2+
startTime, dEff, fLower)
3+
%
4+
% IMRSAWaveform: Generate the parametrised frequency domain BBH coalescence
5+
% waveforms proposed in the paper P. Ajith et al (2009)
6+
%
7+
% usage: [h, Aeff, psiEff] = IMRSAWaveform (f, totalMass, eta, chi, startPhase,
8+
% startTime, dEff, fLower)
9+
%
10+
% fVec : a vector of frequencies at which the waveform is generated
11+
% totalMass : total mass of the binary (in solar masses)
12+
% eta : symmetric mass ratio
13+
% chi : spin parameter
14+
% startPhase : start phase of the waveform (in radian)
15+
% startTime : start time of the waveform (in seconds)
16+
% dEff : effective distance (in Mpc)
17+
% fLower : low-frequency cutoff (Hz)
18+
%
19+
% h : waveform in Fourier domain (complex vector)
20+
% Aeff : amplitude in the Fourier domain
21+
% psiEff : phase in the Fourier domain
22+
%
23+
% P. Ajith, 30.06.2010
24+
%
25+
% $Id: IMRSAWaveform.m 73 2010-07-01 21:39:26Z ajith $
26+
27+
MSOLAR_TIME = 4.92579497077314e-06;
28+
PARSEC_SEC = 1.0292712503e8;
29+
30+
% parameters
31+
M = totalMass*MSOLAR_TIME;
32+
piM = pi*M;
33+
d = dEff*PARSEC_SEC*1e6;
34+
phi = startPhase;
35+
shft = 2*pi*startTime;
36+
delta = sqrt(1.-4.*eta);
37+
mergPower = -2/3;
38+
39+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
40+
% compute the phenomenological parameters
41+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42+
psi0 = 3./(128*eta);
43+
44+
psi1 = 0.;
45+
46+
psi2 = 3715./756. + ...
47+
-9.2091e+02*eta^1 + 4.9213e+02*eta^1*chi^1 + 1.3503e+02*eta^1*chi^2 + ...
48+
6.7419e+03*eta^2 + -1.0534e+03*eta^2*chi^1 + ...
49+
-1.3397e+04*eta^3 ;
50+
51+
psi3 = -16.*pi + 113.*chi/3. + ...
52+
1.7022e+04*eta^1 + -9.5659e+03*eta^1*chi^1 + -2.1821e+03*eta^1*chi^2 + ...
53+
-1.2137e+05*eta^2 + 2.0752e+04*eta^2*chi^1 + ...
54+
2.3859e+05*eta^3 ;
55+
56+
psi4 = 15293365./508032. - 405.*chi^2/8. + ...
57+
-1.2544e+05*eta^1 + 7.5066e+04*eta^1*chi^1 + 1.3382e+04*eta^1*chi^2 + ...
58+
8.7354e+05*eta^2 + -1.6573e+05*eta^2*chi^1 + ...
59+
-1.6936e+06*eta^3 ;
60+
61+
psi5 = 0.;
62+
63+
psi6 = -8.8977e+05*eta^1 + 6.3102e+05*eta^1*chi^1 + 5.0676e+04*eta^1*chi^2 + ...
64+
5.9808e+06*eta^2 + -1.4148e+06*eta^2*chi^1 + ...
65+
-1.1280e+07*eta^3 ;
66+
67+
psi7 = 8.6960e+05*eta^1 + -6.7098e+05*eta^1*chi^1 + -3.0082e+04*eta^1*chi^2 + ...
68+
-5.8379e+06*eta^2 + 1.5145e+06*eta^2*chi^1 + ...
69+
1.0891e+07*eta^3 ;
70+
71+
psi8 = -3.6600e+05*eta^1 + 3.0670e+05*eta^1*chi^1 + 6.3176e+02*eta^1*chi^2 + ...
72+
2.4265e+06*eta^2 + -7.2180e+05*eta^2*chi^1 + ...
73+
-4.5524e+06*eta^3 ;
74+
75+
f1 = 1. - 4.4547*(1.-chi).^0.217 + 3.521*(1.-chi).^0.26 + ...
76+
6.4365e-01*eta^1 + 8.2696e-01*eta^1*chi^1 + -2.7063e-01*eta^1*chi^2 + ...
77+
-5.8218e-02*eta^2 + -3.9346e+00*eta^2*chi^1 + ...
78+
-7.0916e+00*eta^3 ;
79+
80+
f2 = (1. - 0.63*(1.-chi)^0.3)/2. + ...
81+
1.4690e-01*eta^1 + -1.2281e-01*eta^1*chi^1 + -2.6091e-02*eta^1*chi^2 + ...
82+
-2.4900e-02*eta^2 + 1.7013e-01*eta^2*chi^1 + ...
83+
2.3252e+00*eta^3 ;
84+
85+
sigma = (1. - 0.63*(1.-chi)^0.3)*(1.-chi)^0.45/4. + ...
86+
-4.0979e-01*eta^1 + -3.5226e-02*eta^1*chi^1 + 1.0082e-01*eta^1*chi^2 + ...
87+
1.8286e+00*eta^2 + -2.0169e-02*eta^2*chi^1 + ...
88+
-2.8698e+00*eta^3 ;
89+
90+
f3 = 3.2361e-01 + 4.8935e-02*chi^1 + 1.3463e-02*chi^2 + ...
91+
-1.3313e-01*eta^1 + -8.1719e-02*eta^1*chi^1 + 1.4512e-01*eta^1*chi^2 + ...
92+
-2.7140e-01*eta^2 + 1.2788e-01*eta^2*chi^1 + ...
93+
4.9220e+00*eta^3 ;
94+
95+
f3 = f3/piM;
96+
f1 = f1/piM;
97+
f2 = f2/piM;
98+
sigma = sigma/piM;
99+
100+
%fprintf('\nPhenomenological parameters:\n');
101+
%fprintf('psi0 = %5.4e \t psi1 = %5.4e \t psi2 = %5.4e \t psi3 = %5.4e \t psi4 = %5.4e\n',...
102+
% psi0, psi1, psi2, psi3, psi4);
103+
%fprintf('psi5 = %5.4e \t psi6 = %5.4e \t psi7 = %5.4e \t psi8 = %5.4e\n',...
104+
% psi5, psi6, psi7, psi8);
105+
%fprintf('f1 = %5.4e \t f2 = %5.4e \t f3 = %5.4e \t sigma = %5.4e\n',...
106+
% f1, f2, f3, sigma);
107+
108+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
109+
110+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
111+
% now generate the waveforms
112+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
113+
114+
% PN corrections to the inspiral amplitude
115+
alpha2 = -323./224. + 451.*eta/168.;
116+
alpha3 = (27./8. - 11.*eta/6.)*chi;
117+
118+
% correction to the merger power law
119+
epsilon_1 = 1.4547*chi - 1.8897;
120+
epsilon_2 = -1.8153*chi + 1.6557;
121+
122+
% normalisation constants for the merger and ring down amplitude
123+
vMerg = power(pi*M*f1, 1./3.);
124+
vRing = power(pi*M*f2, 1./3.);
125+
w1 = 1. + alpha2*power(vMerg,2.) + alpha3*power(vMerg,3.);
126+
w1 = w1/(1. + epsilon_1*vMerg + epsilon_2*vMerg*vMerg);
127+
w2 = w1*(pi*sigma/2.)*power(f2/f1, mergPower)*(1. ...
128+
+ epsilon_1*vRing + epsilon_2*vRing*vRing);
129+
130+
% amplitude scale
131+
C = M^(5/6)*(f1)^(-7/6)*sqrt(5*eta/24)/(d*pi^(2/3));
132+
133+
Aeff = zeros(size(fVec));
134+
psiEff = zeros(size(fVec));
135+
136+
inspIdx = find(fVec >= fLower & fVec < f1);
137+
mergIdx = find(fVec >= f1 & fVec < f2);
138+
ringIdx = find(fVec >= f2 & fVec <= f3);
139+
outOfBand = find(fVec < fLower | fVec > f3);
140+
141+
v = power(pi*M*fVec, 1./3.);
142+
143+
% effective amplitude - inspiral, merger and ring down
144+
fNorm = fVec/f1;
145+
pnCorr = 1 + alpha2*power(v(inspIdx),2.) + alpha3*power(v(inspIdx),3.);
146+
147+
Aeff(inspIdx) = power(fNorm(inspIdx), -7./6.).*pnCorr;
148+
Aeff(mergIdx) = w1*power(fNorm(mergIdx), mergPower).*(1. ...
149+
+ epsilon_1*v(mergIdx) + epsilon_2*power(v(mergIdx),2));
150+
151+
Lorentzian = sigma./(power(fVec(ringIdx)-f2, 2.) + sigma.*sigma/4.);
152+
Aeff(ringIdx) = w2*Lorentzian/(2*pi);
153+
154+
Aeff = C*Aeff;
155+
psiEff = shft*fVec + phi + psi0*v.^-5.*(1 + psi2*v.^2 + psi3*v.^3 + ...
156+
psi4*v.^4 + psi5*v.^5 + psi6*v.^6 + psi7*v.^7 + psi8*v.^8);
157+
158+
h = Aeff.*exp(1i*(-psiEff));
159+

0 commit comments

Comments
 (0)