Skip to content

Commit 3b9e632

Browse files
author
Ilya Mandel
committed
Matlab cosmic integration updates
1 parent effdbcd commit 3b9e632

1 file changed

Lines changed: 125 additions & 58 deletions

File tree

‎compas_matlab_utils/CosmicHistoryIntegrator.m‎

Lines changed: 125 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -1,39 +1,55 @@
1-
function [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, ...
2-
MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, Rdetections, DetectableMergerRate]=...
3-
CosmicHistoryIntegrator(filename, zlistformation, zmaxdetection, Msimulated, makeplots)
1+
function [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByBinary, ...
2+
FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, MergerRateByRedshiftByBinary,...
3+
MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, pdetectionByRedshiftByMtByEta, ...
4+
DetectableMergerRateByRedshiftByBinary, DetectableMergerRateByRedshiftByMtByEta, ...
5+
RdetectionsByRedshiftByBinary, RdetectionsByRedshiftByMtByEta, RdetectionsPerfectDetectorByRedshiftByBinary]=...
6+
CosmicHistoryIntegrator(filename, noisefile, zlistformation, zmaxdetection, Msimulated, makeplots, data)
47
% Integrator for the binary black hole merger rate over cosmic history
58
% COMPAS (Compact Object Mergers: Population Astrophysics and Statistics)
69
% software package
710
%
811
% USAGE:
9-
% [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, ...
10-
% MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, Rdetections, DetectableMergerRate]]=...
11-
% CosmicHistoryIntegrator(filename, zlistformation, zmaxdetection, Msimulated [,makeplots])
12+
% [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByBinary, ...
13+
% FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, MergerRateByRedshiftByBinary,...
14+
% MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, pdetectionByRedshiftByMtByEta, ...
15+
% DetectableMergerRateByRedshiftByBinary, DetectableMergerRateByRedshiftByMtByEta, ...
16+
% RdetectionsByRedshiftByBinary, RdetectionsByRedshiftByMtByEta, RdetectionsPerfectDetectorByRedshiftByBinary]=...
17+
% CosmicHistoryIntegrator(filename, noisefile, zlistformation, zmaxdetection, Msimulated, Msimulated [,makeplots, data])
1218
%
1319
% INPUTS:
1420
% filename: name of population synthesis input file
1521
% should be in COMPAS output h5 format
22+
% noisefile: file containing noise ASD (first column frequency, second ASD)
1623
% zlistformation: vector of redshifts at which the formation rate is
1724
% computed
1825
% zmaxdetection: maximum redshift to which the detection rate is computed
1926
% Msimulated: total star forming mass represented by the simulation (for
2027
% normalisation)
2128
% makeplots: if set to 1, generates a set of useful plots (default = 0)
29+
% data: if provided, this contains a list of simulated binaries as a matrix,
30+
% with each row containing a binary and the columns containing M1, M2, Z and
31+
% tdelay in that order; the filename is then ignored
2232
%
2333
% OUTPUTS:
2434
% SFR is a vector of size length(zlistformation) containing the star formation rate
2535
% (solar masses per Mpc^3 of comoving volume per year of source time)
2636
% Zlist is a vector of metallicities, taken from the COMPAS run input file
2737
% Mtlist is a list of total mass bins
2838
% etalist is a list of symmetric mass ratio bins
39+
% FormationRateByRedshiftByBinary is a matrix of size length(zformationlist) X length(M1)
40+
% which contains a formation rate of merging double compact objects just
41+
% like this binary, in units of formed DCOs per Mpc^3 of comoving volume per year of source time
2942
% FormationRateByRedshiftByZ is a matrix of size length(zformationlist) X length(Zlist)
30-
% which contains a formation rate of merging double compact objects in the given redshift
43+
% which contains the formation rate of merging double compact objects in the given redshift
3144
% and metallicity bin, in units of formed DCOs per Mpc^3 of comoving volume per
3245
% year of source time
3346
% FormationRateByRedshiftByMtByEta is a matrix of size length(zformationlist)
3447
% X length(Mtlist) X length(etalist) which contains a formation rate of merging double compact objects
3548
% in the given redshift, total mass and eta bin, in units of formed DCOs per Mpc^3
3649
% of comoving volume per year of source time
50+
% MergerRateByRedshiftByBinary is a matrix of size length(zformationlist) X length(M1)
51+
% which contains the merger rate of merging double compact objects just
52+
% like this binary, in units of mergers per Mpc^3 of comoving volume per year of source time
3753
% MergerRateByRedshiftByZ is a matrix of size length(zformationlist) X length(Zlist)
3854
% which contains a merger rate of double compact objects in the given redshift
3955
% and metallicity bin, in units of mergers per Mpc^3 of comoving volume per
@@ -44,21 +60,39 @@
4460
% of comoving volume per year of source time
4561
% zlistdetection is a vector of redshifts at which detection rates are
4662
% computed (a subset of zlistformation going up to zmaxdetection)
47-
% Rdetection is a matrix of size length(zlistdetection) X length(Mtlist) X
48-
% length(etalist) containing the detection rate per year of observer time
49-
% from a given redshift bin and total mass and symmetric mass ratio pixel
50-
% DetectableMergerRate is a matrix of the same size as Rdetection but
51-
% containing the intrinsic rate of detectable mergers per Mpc^3 of comoving
63+
% pdetectionByRedshiftByMtByEta is a matrix of size length(zlistdetection) X length(Mtlist) X length(etalist)
64+
% containing the probability that a binary of a given mass and mass ratio
65+
% is detectable for a given merger redshift
66+
% DetectableMergerRateByRedshiftByBinary is a matrix of size
67+
% length(zlistdetection) X length(M1), containing the intrinsic rate of mergers of
68+
% binaries just like this one per Mpc^3 of comoving
5269
% volume per year of source time
70+
% DetectableMergerRateByRedshiftByMtByEta is a matrix of size
71+
% length(zlistdetection) X length(Mtlist) X length(etalist) containing the detection rate per year of observer time
72+
% from a given redshift bin and total mass and symmetric mass ratio pixel
73+
% RdetectionsByRedshiftByBinary is a matrix of the same size as
74+
% DetectableMergerRateByRedshiftByBinary containing the detection rate
75+
% for binaries just like this one per year of observer time
76+
% from a given redshift bin
77+
% RdetectionsByRedshiftByMtByEta is a matrix of the same size as
78+
% DetectableMergerRateByRedshiftByMtByEta containing the detection rate per
79+
% year of observer time from a given redshift bin and total mass and symmetric mass ratio pixel
80+
% RdetectionsPerfectDetectorByRedshiftByBinary is a matrix of the same size
81+
% as RdetectionsByRedshiftByBinary containing the detection rate per year
82+
% of observer time from a given redshift bin for binaries just like this one
83+
% assuming an imaginary perfectly sensitive detector
5384

5485
%
5586
% EXAMPLE:
5687
% zlist=0:0.01:10;
57-
% [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, ...
58-
% MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, Rdetections, DetectableMergerRate]=...
59-
% CosmicHistoryIntegrator('~/Work/COMPASresults/runs/Zdistalpha1-031803.h5', zlist, 1.5, 90e6, 1);
88+
% [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByBinary, ...
89+
% FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, MergerRateByRedshiftByBinary,...
90+
% MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, pdetectionByRedshiftByMtByEta, ...
91+
% DetectableMergerRateByRedshiftByBinary, DetectableMergerRateByRedshiftByMtByEta, ...
92+
% RdetectionsByRedshiftByBinary, RdetectionsByRedshiftByMtByEta, RdetectionsPerfectDetectorByRedshiftByBinary]=...
93+
% CosmicHistoryIntegrator('~/Work/COMPASresults/runs/Zdistalpha1-031803.h5', '~/Work/Rai/aligo_O4high.txt', zlist, 1.5, 90e6, 1);
6094
% figure(10), semilogy(zlist, sum(MergerRateByRedshiftByZ,2)*1e9,'LineWidth',3), set(gca,'FontSize',20),
61-
% xlabel('Redshift z'), ylabel('Formation rate of merging DCO per Gpc^3 per yr')
95+
% xlabel('Redshift z'), ylabel('Merger rate of DCO per Gpc^3 per yr')
6296
%
6397

6498

@@ -71,15 +105,26 @@
71105
Mpc=Mpcm/c; %Mpc in seconds
72106
yr=3.15569e7; %year in seconds
73107

74-
if (nargin<4)
108+
if (nargin<5)
75109
error('Not enough input arguments.');
76110
end;
77111
if (nargin<5), makeplots=0; end;
112+
if (nargin<7),
113+
%load COMPAS data
114+
[M1,M2,Z,Tdelay]=DataRead(filename);
115+
else
116+
if(size(data,2)~=4),
117+
error('The data must have 4 columns: M1, M2, Z, Tdelay');
118+
end;
119+
M1=data(:,1); M2=data(:,2); Z=data(:,3); Tdelay=data(:,4);
120+
%maxNS=0.0; %pretend all are BBH if reading from file
121+
end;
122+
123+
78124

79125
%cosmology calculator
80126
[tL,Dl,dVc]=Cosmology(zlistformation);
81-
%load COMPAS data
82-
[M1,M2,Z,Tdelay,maxNS]=DataRead(filename);
127+
83128
%metallicity-specific SFR
84129
[SFR,Zlist,Zweight]=Metallicity(zlistformation,min(Z),max(Z));
85130

@@ -90,44 +135,72 @@
90135
dz=zlistformation(2)-zlistformation(1);
91136
etalist=0.01:0.01:0.25;
92137
Mtlist=1:1:ceil(max(M1+M2));
138+
FormationRateByRedshiftByBinary=zeros(length(zlistformation),length(M1));
93139
FormationRateByRedshiftByZ=zeros(length(zlistformation),length(Zlist));
94140
FormationRateByRedshiftByMtByEta=zeros(length(zlistformation),length(Mtlist),length(etalist));
141+
MergerRateByRedshiftByBinary=zeros(length(zlistformation),length(M1));
95142
MergerRateByRedshiftByZ=zeros(length(zlistformation),length(Zlist));
96143
MergerRateByRedshiftByMtByEta=zeros(length(zlistformation),length(Mtlist),length(etalist));
97-
x=zeros(size(M1));
144+
etaindexByBinary=zeros(size(M1)); MtindexByBinary=zeros(size(M1));
98145
for(i=1:length(M1)),
99146
Zcounter=find(Zlist>=Z(i),1);
100147
eta=M1(i)*M2(i)/(M1(i)+M2(i))^2;
101-
etaindex=ceil(eta*100);
102-
Mtindex=ceil(M1(i)+M2(i));
103-
FormationRateByRedshiftByZ(:,Zcounter)=transpose(SFR).*Zweight(:,Zcounter)/Msimulated;
104-
FormationRateByRedshiftByMtByEta(:,Mtindex,etaindex)=transpose(SFR).*Zweight(:,Zcounter)/Msimulated;
148+
etaindex=ceil(eta*100); etaindexByBinary(i)=etaindex;
149+
Mtindex=ceil(M1(i)+M2(i)); MtindexByBinary(i)=Mtindex;
150+
FormationRateByRedshiftByBinary(:,i)=transpose(SFR).*Zweight(:,Zcounter)/Msimulated;
151+
FormationRateByRedshiftByZ(:,Zcounter)=FormationRateByRedshiftByZ(:,Zcounter)+...
152+
FormationRateByRedshiftByBinary(:,i);
153+
FormationRateByRedshiftByMtByEta(:,Mtindex,etaindex)=FormationRateByRedshiftByMtByEta(:,Mtindex,etaindex)+...
154+
FormationRateByRedshiftByBinary(:,i);
105155
tLform=tL+Tdelay(i); %lookback time of when binary would have to form in order to merge at lookback time tL
106156
firsttooearlyindex=find((tLform)>max(tL),1);
107157
if(isempty(firsttooearlyindex)), firsttooearlyindex=length(tL)+1; end;
108158
zForm=interp1(tL,zlistformation,tLform(1:firsttooearlyindex-1));
109159
zFormindex=ceil((zForm-zlistformation(1))./dz)+1;
110160
if(~isempty(zFormindex))
111-
x(i)=SFR(zFormindex(1))*Zweight(zFormindex(1),Zcounter)/Msimulated;
161+
MergerRateByRedshiftByBinary(1:firsttooearlyindex-1,i)=...
162+
transpose(SFR(zFormindex)).*Zweight(zFormindex,Zcounter)/Msimulated;
112163
MergerRateByRedshiftByZ(1:firsttooearlyindex-1,Zcounter)=...
113164
MergerRateByRedshiftByZ(1:firsttooearlyindex-1,Zcounter)+...
114-
transpose(SFR(zFormindex)).*Zweight(zFormindex,Zcounter)/Msimulated;
165+
MergerRateByRedshiftByBinary(1:firsttooearlyindex-1,i);
115166
MergerRateByRedshiftByMtByEta(1:firsttooearlyindex-1,Mtindex,etaindex) =...
116167
MergerRateByRedshiftByMtByEta(1:firsttooearlyindex-1,Mtindex,etaindex) + ...
117-
transpose(SFR(zFormindex)).*Zweight(zFormindex,Zcounter)/Msimulated;
168+
MergerRateByRedshiftByBinary(1:firsttooearlyindex-1,i);
118169
end;
119170
end;
120171

121172
zlistdetection=zlistformation(1:find(zlistformation<=zmaxdetection,1,"last"));
122-
fin=load('~/Work/Rai/LIGOfuture_data/freqVector.txt');
123-
%noise=load('~/Work/Rai/LIGOfuture_data/dataNomaLIGO.txt');
124-
noise=load('~/Work/Rai/LIGOfuture_data/dataEarly_low.txt');
125-
[Rdetections,DetectableMergerRate]=...
126-
DetectionRate(zlistformation,Mtlist,etalist,MergerRateByRedshiftByMtByEta,zlistdetection,fin,noise,Dl,dVc);
173+
174+
%detection probability
175+
pdetectionByRedshiftByMtByEta=...
176+
DetectionProbability(zlistdetection,Mtlist,etalist,noisefile,Dl);
177+
pdetectionsByZByBinary=zeros(length(zlistdetection),length(M1));
178+
for(i=1:length(M1)),
179+
pdetectionByRedshiftByBinary(:,i)=pdetectionByRedshiftByMtByEta(:,MtindexByBinary(i),etaindexByBinary(i));
180+
end;
181+
%Detections per unit source time per unit Vc in each Mt and eta bin
182+
DetectableMergerRateByRedshiftByMtByEta=...
183+
MergerRateByRedshiftByMtByEta(1:length(zlistdetection),:,:).*pdetectionByRedshiftByMtByEta;
184+
%Detections per unit source time per unit Vc for each binary
185+
DetectableMergerRateByRedshiftByBinary=...
186+
MergerRateByRedshiftByBinary(1:length(zlistdetection),:).*pdetectionByRedshiftByBinary;
187+
%Detections per unit observer time in each Mt and eta bin
188+
RdetectionsByRedshiftByMtByEta=...
189+
DetectableMergerRateByRedshiftByMtByEta.*transpose(dVc(1:length(zlistdetection)))./(1+zlistdetection');
190+
%Detections per unit observer time by binary
191+
RdetectionsByRedshiftByBinary=...
192+
DetectableMergerRateByRedshiftByBinary.*transpose(dVc(1:length(zlistdetection)))./(1+zlistdetection');
193+
%Detections per unit observer time per binary by an imaginary perfect detector
194+
RdetectionsPerfectDetectorByRedshiftByBinary=...
195+
MergerRateByRedshiftByBinary(1:length(zlistdetection),:).*...
196+
transpose(dVc(1:length(zlistdetection)))./(1+zlistdetection');
197+
%total detection rate
198+
disp(['Total detection rate: ', ...
199+
num2str(sum(sum(RdetectionsByRedshiftByBinary))), ' per year']);
127200

128201
if(makeplots==1), %make a set of default plots
129202
MakePlots(M1,M2,Z,Tdelay,zlistformation,Zlist,SFR,Zweight,...
130-
MergerRateByRedshiftByZ, Rdetections, DetectableMergerRate, zlistdetection, Mtlist, etalist, 1);
203+
MergerRateByRedshiftByZ, RdetectionsByRedshiftByMtByEta, DetectableMergerRateByRedshiftByMtByEta, zlistdetection, Mtlist, etalist, 1);
131204
end;
132205

133206
end %end of CosmicHistoryIntegrator
@@ -136,7 +209,7 @@
136209
%Load the data stored in COMPAS .h5 output format from a file
137210
%Select only double compact object mergers of interest, and return the
138211
%component masses, metallicities, and star formation to merger delay times
139-
function [M1,M2,Z,Tdelay, maxNS]=DataRead(file)
212+
function [M1,M2,Z,Tdelay]=DataRead(file)
140213
if(exist(file, 'file')~=2),
141214
error('Input file does not exist');
142215
end;
@@ -157,7 +230,7 @@
157230
%NSBH=(((type1==13) & (type2==14)) | ((type1==14) & (type2==13)));
158231
%mergingDCO=mergingBNS | mergingNSBH | mergingBBH;
159232
%BNScount=sum(mergingBNS); NSBHcount=sum(mergingNSBH); BBHcount=sum(mergingBBH);
160-
maxNS=max(max(mass1(type1==13)), max(mass2(type2==13)));
233+
%maxNS=max(max(mass1(type1==13)), max(mass2(type2==13)));
161234
chirpmass=mass1.^0.6.*mass2.^0.6./(mass1+mass2).^0.2;
162235
q=mass2./mass1;
163236
seedCE=h5read(file,'/BSE_Common_Envelopes/SEED');
@@ -235,16 +308,17 @@
235308

236309
%Compute detection rates per year of observer time and per year of source time
237310
%per Mpc^3 of comoving volume as a function of total mass and eta
238-
function [Rdetections, DetectableMergerRate]=...
239-
DetectionRate(zlistformation,Mtlist,etalist,MergerRateByRedshiftByMtByEta,zlistdetection,freqfile,noisefile,Dl,dVc)
311+
%Compute the detection probability as a function of detection redshift,
312+
%total mass and eta
313+
function [pdetection]=...
314+
DetectionProbability(zlistdetection,Mtlist,etalist,noisefile,Dl)
240315

241-
fin=load('~/Work/Rai/LIGOfuture_data/freqVector.txt');
242-
noise=load('~/Work/Rai/LIGOfuture_data/dataMid_low.txt');
316+
noise=load(noisefile);
243317

244-
flow=10;
318+
flow=max(10,ceil(min(noise(:,1))));
245319
df=1;
246320
f=flow:df:500; %BBH focussed
247-
Sf=interp1(fin, noise.^2, f);
321+
Sf=interp1(noise(:,1), noise(:,2).^2, f);
248322

249323
Ntheta=1e6;
250324
psi=rand(1,Ntheta)*pi;
@@ -275,26 +349,19 @@
275349

276350
for(i=1:length(zlistdetection)),
277351
for(j=1:length(Mtlist)),
278-
SNR(i,j,:)=SNRat1Mpc(ceil(j*(zlistdetection(i)+1)),:)./Dl(i);
352+
SNR(i,j,:)=SNRat1Mpc(ceil(Mtlist(j)*(zlistdetection(i)+1)),:)./Dl(i);
279353
end;
280354
end;
281-
%for(i=1:length(zlistdetection)), SNR(i,:,:)=SNRat1Mpc./Dl(i); end;
282355

283-
SNR8pre=1:0.1:1000;
356+
SNR8pre=1:0.01:100;
284357
theta=1./SNR8pre;
285358
pdetect=1-interp1([0,Thetas,1],[(0:Ntheta)/Ntheta,1],theta);
286359
pdetect(1)=0; %set of measure zero to exceed threshold, but enforce just in case
287360

288-
Rdetections=zeros(length(zlistdetection),length(Mtlist),length(etalist)); %Detections per unit observer time
289-
DetectableMergerRate=zeros(length(zlistdetection),length(Mtlist),length(etalist)); %Detections per unit source time per unit Vc
290361
SNR8=SNR/8;
291-
pdetection=zeros(size(Rdetections));
292-
pdetection=pdetect(max(min(floor(SNR8*10),length(pdetect)),1));
293-
294-
DetectableMergerRate=MergerRateByRedshiftByMtByEta(1:length(zlistdetection),:,:).*pdetection;
295-
Rdetections=DetectableMergerRate.*transpose(dVc(1:length(zlistdetection)))./(1+zlistdetection');
296-
297-
end %end of DetectionRate
362+
pdetection=zeros(length(zlistdetection),length(Mtlist),length(etalist));
363+
pdetection=pdetect(max(min(floor((SNR8-1)*100),length(pdetect)),1));
364+
end %end of DetectionProbability
298365

299366
%Make a set of default plots
300367
function MakePlots(M1,M2,Z,Tdelay,zlistformation,Zlist,SFR,Zweight,...
@@ -341,11 +408,11 @@ function MakePlots(M1,M2,Z,Tdelay,zlistformation,Zlist,SFR,Zweight,...
341408

342409

343410
figure(fignumber+4), clf(fignumber+4);
344-
RdetectionsByzMt=sum(Rdetections,3); %sum across eta
345-
semilogy(zlistdetection, cumsum(sum(RdetectionsByzMt,2)), 'LineWidth', 3), hold on;
346-
semilogy(zlistdetection, cumsum(sum(RdetectionsByzMt(:,Mtlist<=5),2)), 'LineWidth', 1);
347-
semilogy(zlistdetection, cumsum(sum(RdetectionsByzMt(:,Mtlist>5 & Mtlist<20),2)), 'LineWidth', 1);
348-
semilogy(zlistdetection, cumsum(sum(RdetectionsByzMt(:,Mtlist>=20),2)), 'LineWidth', 1); hold off;
411+
RdetectionsByRedshiftMt=sum(Rdetections,3); %sum across eta
412+
semilogy(zlistdetection, cumsum(sum(RdetectionsByRedshiftMt,2)), 'LineWidth', 3), hold on;
413+
semilogy(zlistdetection, cumsum(sum(RdetectionsByRedshiftMt(:,Mtlist<=5),2)), 'LineWidth', 1);
414+
semilogy(zlistdetection, cumsum(sum(RdetectionsByRedshiftMt(:,Mtlist>5 & Mtlist<20),2)), 'LineWidth', 1);
415+
semilogy(zlistdetection, cumsum(sum(RdetectionsByRedshiftMt(:,Mtlist>=20),2)), 'LineWidth', 1); hold off;
349416
legend('Total rate', 'From M_t<=5 M_o', 'From 5<M_t/M_o<20', 'From M_t>=20 M_o'),
350417
set(gca, 'FontSize', 20); %for labels
351418
xlabel('z'),

0 commit comments

Comments
 (0)