Skip to content

Matlab cosmic integration updates #1480

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jeffriley merged 2 commits into from
May 10, 2026
+126 −58
Merged

Matlab cosmic integration updates #1480

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3b9e632
Matlab cosmic integration updates
May 3, 2026
34efe33
Changed ASD to PSD
May 8, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
184 changes: 126 additions & 58 deletions compas_matlab_utils/CosmicHistoryIntegrator.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,39 +1,55 @@
function [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, ...
MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, Rdetections, DetectableMergerRate]=...
CosmicHistoryIntegrator(filename, zlistformation, zmaxdetection, Msimulated, makeplots)
function [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByBinary, ...
FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, MergerRateByRedshiftByBinary,...
MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, pdetectionByRedshiftByMtByEta, ...
DetectableMergerRateByRedshiftByBinary, DetectableMergerRateByRedshiftByMtByEta, ...
RdetectionsByRedshiftByBinary, RdetectionsByRedshiftByMtByEta, RdetectionsPerfectDetectorByRedshiftByBinary]=...
CosmicHistoryIntegrator(filename, noisefile, zlistformation, zmaxdetection, Msimulated, makeplots, data)
% Integrator for the binary black hole merger rate over cosmic history
% COMPAS (Compact Object Mergers: Population Astrophysics and Statistics)
% software package
%
% USAGE:
% [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, ...
% MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, Rdetections, DetectableMergerRate]]=...
% CosmicHistoryIntegrator(filename, zlistformation, zmaxdetection, Msimulated [,makeplots])
% [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByBinary, ...
% FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, MergerRateByRedshiftByBinary,...
% MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, pdetectionByRedshiftByMtByEta, ...
% DetectableMergerRateByRedshiftByBinary, DetectableMergerRateByRedshiftByMtByEta, ...
% RdetectionsByRedshiftByBinary, RdetectionsByRedshiftByMtByEta, RdetectionsPerfectDetectorByRedshiftByBinary]=...
% CosmicHistoryIntegrator(filename, noisefile, zlistformation, zmaxdetection, Msimulated, Msimulated [,makeplots, data])
%
% INPUTS:
% filename: name of population synthesis input file
% should be in COMPAS output h5 format
% noisefile: file containing noise PSD (first column frequency, second PSD)
% zlistformation: vector of redshifts at which the formation rate is
% computed
% zmaxdetection: maximum redshift to which the detection rate is computed
% Msimulated: total star forming mass represented by the simulation (for
% normalisation)
% makeplots: if set to 1, generates a set of useful plots (default = 0)
% data: if provided, this contains a list of simulated binaries as a matrix,
% with each row containing a binary and the columns containing M1, M2, Z and
% tdelay in that order; the filename is then ignored
%
% OUTPUTS:
% SFR is a vector of size length(zlistformation) containing the star formation rate
% (solar masses per Mpc^3 of comoving volume per year of source time)
% Zlist is a vector of metallicities, taken from the COMPAS run input file
% Mtlist is a list of total mass bins
% etalist is a list of symmetric mass ratio bins
% FormationRateByRedshiftByBinary is a matrix of size length(zformationlist) X length(M1)
% which contains a formation rate of merging double compact objects just
% like this binary, in units of formed DCOs per Mpc^3 of comoving volume per year of source time
% FormationRateByRedshiftByZ is a matrix of size length(zformationlist) X length(Zlist)
% which contains a formation rate of merging double compact objects in the given redshift
% which contains the formation rate of merging double compact objects in the given redshift
% and metallicity bin, in units of formed DCOs per Mpc^3 of comoving volume per
% year of source time
% FormationRateByRedshiftByMtByEta is a matrix of size length(zformationlist)
% X length(Mtlist) X length(etalist) which contains a formation rate of merging double compact objects
% in the given redshift, total mass and eta bin, in units of formed DCOs per Mpc^3
% of comoving volume per year of source time
% MergerRateByRedshiftByBinary is a matrix of size length(zformationlist) X length(M1)
% which contains the merger rate of merging double compact objects just
% like this binary, in units of mergers per Mpc^3 of comoving volume per year of source time
% MergerRateByRedshiftByZ is a matrix of size length(zformationlist) X length(Zlist)
% which contains a merger rate of double compact objects in the given redshift
% and metallicity bin, in units of mergers per Mpc^3 of comoving volume per
Expand All @@ -44,21 +60,39 @@
% of comoving volume per year of source time
% zlistdetection is a vector of redshifts at which detection rates are
% computed (a subset of zlistformation going up to zmaxdetection)
% Rdetection is a matrix of size length(zlistdetection) X length(Mtlist) X
% length(etalist) containing the detection rate per year of observer time
% from a given redshift bin and total mass and symmetric mass ratio pixel
% DetectableMergerRate is a matrix of the same size as Rdetection but
% containing the intrinsic rate of detectable mergers per Mpc^3 of comoving
% pdetectionByRedshiftByMtByEta is a matrix of size length(zlistdetection) X length(Mtlist) X length(etalist)
% containing the probability that a binary of a given mass and mass ratio
% is detectable for a given merger redshift
% DetectableMergerRateByRedshiftByBinary is a matrix of size
% length(zlistdetection) X length(M1), containing the intrinsic rate of mergers of
% binaries just like this one per Mpc^3 of comoving
% volume per year of source time
% DetectableMergerRateByRedshiftByMtByEta is a matrix of size
% length(zlistdetection) X length(Mtlist) X length(etalist) containing the detection rate per year of observer time
% from a given redshift bin and total mass and symmetric mass ratio pixel
% RdetectionsByRedshiftByBinary is a matrix of the same size as
% DetectableMergerRateByRedshiftByBinary containing the detection rate
% for binaries just like this one per year of observer time
% from a given redshift bin
% RdetectionsByRedshiftByMtByEta is a matrix of the same size as
% DetectableMergerRateByRedshiftByMtByEta containing the detection rate per
% year of observer time from a given redshift bin and total mass and symmetric mass ratio pixel
% RdetectionsPerfectDetectorByRedshiftByBinary is a matrix of the same size
% as RdetectionsByRedshiftByBinary containing the detection rate per year
% of observer time from a given redshift bin for binaries just like this one
% assuming an imaginary perfectly sensitive detector

%
% EXAMPLE:
% zlist=0:0.01:10;
% [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, ...
% MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, Rdetections, DetectableMergerRate]=...
% CosmicHistoryIntegrator('~/Work/COMPASresults/runs/Zdistalpha1-031803.h5', zlist, 1.5, 90e6, 1);
% [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByBinary, ...
% FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, MergerRateByRedshiftByBinary,...
% MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, pdetectionByRedshiftByMtByEta, ...
% DetectableMergerRateByRedshiftByBinary, DetectableMergerRateByRedshiftByMtByEta, ...
% RdetectionsByRedshiftByBinary, RdetectionsByRedshiftByMtByEta, RdetectionsPerfectDetectorByRedshiftByBinary]=...
% CosmicHistoryIntegrator('~/Work/COMPASresults/runs/Zdistalpha1-031803.h5', '~/Work/Rai/psd-O4-2023_06_v1-L.txt', zlist, 1.5, 90e6, 1);
% figure(10), semilogy(zlist, sum(MergerRateByRedshiftByZ,2)*1e9,'LineWidth',3), set(gca,'FontSize',20),
% xlabel('Redshift z'), ylabel('Formation rate of merging DCO per Gpc^3 per yr')
% xlabel('Redshift z'), ylabel('Merger rate of DCO per Gpc^3 per yr')
%


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

if (nargin<4)
if (nargin<5)
error('Not enough input arguments.');
end;
if (nargin<5), makeplots=0; end;
if (nargin<7),
%load COMPAS data
[M1,M2,Z,Tdelay]=DataRead(filename);
else
if(size(data,2)~=4),
error('The data must have 4 columns: M1, M2, Z, Tdelay');
end;
M1=data(:,1); M2=data(:,2); Z=data(:,3); Tdelay=data(:,4);
%maxNS=0.0; %pretend all are BBH if reading from file
end;



%cosmology calculator
[tL,Dl,dVc]=Cosmology(zlistformation);
%load COMPAS data
[M1,M2,Z,Tdelay,maxNS]=DataRead(filename);

%metallicity-specific SFR
[SFR,Zlist,Zweight]=Metallicity(zlistformation,min(Z),max(Z));

Expand All @@ -90,44 +135,72 @@
dz=zlistformation(2)-zlistformation(1);
etalist=0.01:0.01:0.25;
Mtlist=1:1:ceil(max(M1+M2));
FormationRateByRedshiftByBinary=zeros(length(zlistformation),length(M1));
FormationRateByRedshiftByZ=zeros(length(zlistformation),length(Zlist));
FormationRateByRedshiftByMtByEta=zeros(length(zlistformation),length(Mtlist),length(etalist));
MergerRateByRedshiftByBinary=zeros(length(zlistformation),length(M1));
MergerRateByRedshiftByZ=zeros(length(zlistformation),length(Zlist));
MergerRateByRedshiftByMtByEta=zeros(length(zlistformation),length(Mtlist),length(etalist));
x=zeros(size(M1));
etaindexByBinary=zeros(size(M1)); MtindexByBinary=zeros(size(M1));
for(i=1:length(M1)),
Zcounter=find(Zlist>=Z(i),1);
eta=M1(i)*M2(i)/(M1(i)+M2(i))^2;
etaindex=ceil(eta*100);
Mtindex=ceil(M1(i)+M2(i));
FormationRateByRedshiftByZ(:,Zcounter)=transpose(SFR).*Zweight(:,Zcounter)/Msimulated;
FormationRateByRedshiftByMtByEta(:,Mtindex,etaindex)=transpose(SFR).*Zweight(:,Zcounter)/Msimulated;
etaindex=ceil(eta*100); etaindexByBinary(i)=etaindex;
Mtindex=ceil(M1(i)+M2(i)); MtindexByBinary(i)=Mtindex;
FormationRateByRedshiftByBinary(:,i)=transpose(SFR).*Zweight(:,Zcounter)/Msimulated;
FormationRateByRedshiftByZ(:,Zcounter)=FormationRateByRedshiftByZ(:,Zcounter)+...
FormationRateByRedshiftByBinary(:,i);
FormationRateByRedshiftByMtByEta(:,Mtindex,etaindex)=FormationRateByRedshiftByMtByEta(:,Mtindex,etaindex)+...
FormationRateByRedshiftByBinary(:,i);
tLform=tL+Tdelay(i); %lookback time of when binary would have to form in order to merge at lookback time tL
firsttooearlyindex=find((tLform)>max(tL),1);
if(isempty(firsttooearlyindex)), firsttooearlyindex=length(tL)+1; end;
zForm=interp1(tL,zlistformation,tLform(1:firsttooearlyindex-1));
zFormindex=ceil((zForm-zlistformation(1))./dz)+1;
if(~isempty(zFormindex))
x(i)=SFR(zFormindex(1))*Zweight(zFormindex(1),Zcounter)/Msimulated;
MergerRateByRedshiftByBinary(1:firsttooearlyindex-1,i)=...
transpose(SFR(zFormindex)).*Zweight(zFormindex,Zcounter)/Msimulated;
MergerRateByRedshiftByZ(1:firsttooearlyindex-1,Zcounter)=...
MergerRateByRedshiftByZ(1:firsttooearlyindex-1,Zcounter)+...
transpose(SFR(zFormindex)).*Zweight(zFormindex,Zcounter)/Msimulated;
MergerRateByRedshiftByBinary(1:firsttooearlyindex-1,i);
MergerRateByRedshiftByMtByEta(1:firsttooearlyindex-1,Mtindex,etaindex) =...
MergerRateByRedshiftByMtByEta(1:firsttooearlyindex-1,Mtindex,etaindex) + ...
transpose(SFR(zFormindex)).*Zweight(zFormindex,Zcounter)/Msimulated;
MergerRateByRedshiftByBinary(1:firsttooearlyindex-1,i);
end;
end;

zlistdetection=zlistformation(1:find(zlistformation<=zmaxdetection,1,"last"));
fin=load('~/Work/Rai/LIGOfuture_data/freqVector.txt');
%noise=load('~/Work/Rai/LIGOfuture_data/dataNomaLIGO.txt');
noise=load('~/Work/Rai/LIGOfuture_data/dataEarly_low.txt');
[Rdetections,DetectableMergerRate]=...
DetectionRate(zlistformation,Mtlist,etalist,MergerRateByRedshiftByMtByEta,zlistdetection,fin,noise,Dl,dVc);

%detection probability
pdetectionByRedshiftByMtByEta=...
DetectionProbability(zlistdetection,Mtlist,etalist,noisefile,Dl);
pdetectionsByZByBinary=zeros(length(zlistdetection),length(M1));
for(i=1:length(M1)),
pdetectionByRedshiftByBinary(:,i)=pdetectionByRedshiftByMtByEta(:,MtindexByBinary(i),etaindexByBinary(i));
end;
%Detections per unit source time per unit Vc in each Mt and eta bin
DetectableMergerRateByRedshiftByMtByEta=...
MergerRateByRedshiftByMtByEta(1:length(zlistdetection),:,:).*pdetectionByRedshiftByMtByEta;
%Detections per unit source time per unit Vc for each binary
DetectableMergerRateByRedshiftByBinary=...
MergerRateByRedshiftByBinary(1:length(zlistdetection),:).*pdetectionByRedshiftByBinary;
%Detections per unit observer time in each Mt and eta bin
RdetectionsByRedshiftByMtByEta=...
DetectableMergerRateByRedshiftByMtByEta.*transpose(dVc(1:length(zlistdetection)))./(1+zlistdetection');
%Detections per unit observer time by binary
RdetectionsByRedshiftByBinary=...
DetectableMergerRateByRedshiftByBinary.*transpose(dVc(1:length(zlistdetection)))./(1+zlistdetection');
%Detections per unit observer time per binary by an imaginary perfect detector
RdetectionsPerfectDetectorByRedshiftByBinary=...
MergerRateByRedshiftByBinary(1:length(zlistdetection),:).*...
transpose(dVc(1:length(zlistdetection)))./(1+zlistdetection');
%total detection rate
disp(['Total detection rate: ', ...
num2str(sum(sum(RdetectionsByRedshiftByBinary))), ' per year']);

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

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

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

fin=load('~/Work/Rai/LIGOfuture_data/freqVector.txt');
noise=load('~/Work/Rai/LIGOfuture_data/dataMid_low.txt');
noise=load(noisefile);

flow=10;
flow=max(10,ceil(min(noise(:,1))));
df=1;
f=flow:df:500; %BBH focussed
Sf=interp1(fin, noise.^2, f);
Sf=interp1(noise(:,1), noise(:,2), f);
%Sf=interp1(noise(:,1), noise(:,2).^2, f); %for ASDs

Ntheta=1e6;
psi=rand(1,Ntheta)*pi;
Expand Down Expand Up @@ -275,26 +350,19 @@

for(i=1:length(zlistdetection)),
for(j=1:length(Mtlist)),
SNR(i,j,:)=SNRat1Mpc(ceil(j*(zlistdetection(i)+1)),:)./Dl(i);
SNR(i,j,:)=SNRat1Mpc(ceil(Mtlist(j)*(zlistdetection(i)+1)),:)./Dl(i);
end;
end;
%for(i=1:length(zlistdetection)), SNR(i,:,:)=SNRat1Mpc./Dl(i); end;

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

Rdetections=zeros(length(zlistdetection),length(Mtlist),length(etalist)); %Detections per unit observer time
DetectableMergerRate=zeros(length(zlistdetection),length(Mtlist),length(etalist)); %Detections per unit source time per unit Vc
SNR8=SNR/8;
pdetection=zeros(size(Rdetections));
pdetection=pdetect(max(min(floor(SNR8*10),length(pdetect)),1));

DetectableMergerRate=MergerRateByRedshiftByMtByEta(1:length(zlistdetection),:,:).*pdetection;
Rdetections=DetectableMergerRate.*transpose(dVc(1:length(zlistdetection)))./(1+zlistdetection');

end %end of DetectionRate
pdetection=zeros(length(zlistdetection),length(Mtlist),length(etalist));
pdetection=pdetect(max(min(floor((SNR8-1)*100),length(pdetect)),1));
end %end of DetectionProbability

%Make a set of default plots
function MakePlots(M1,M2,Z,Tdelay,zlistformation,Zlist,SFR,Zweight,...
Expand Down Expand Up @@ -341,11 +409,11 @@ function MakePlots(M1,M2,Z,Tdelay,zlistformation,Zlist,SFR,Zweight,...


figure(fignumber+4), clf(fignumber+4);
RdetectionsByzMt=sum(Rdetections,3); %sum across eta
semilogy(zlistdetection, cumsum(sum(RdetectionsByzMt,2)), 'LineWidth', 3), hold on;
semilogy(zlistdetection, cumsum(sum(RdetectionsByzMt(:,Mtlist<=5),2)), 'LineWidth', 1);
semilogy(zlistdetection, cumsum(sum(RdetectionsByzMt(:,Mtlist>5 & Mtlist<20),2)), 'LineWidth', 1);
semilogy(zlistdetection, cumsum(sum(RdetectionsByzMt(:,Mtlist>=20),2)), 'LineWidth', 1); hold off;
RdetectionsByRedshiftMt=sum(Rdetections,3); %sum across eta
semilogy(zlistdetection, cumsum(sum(RdetectionsByRedshiftMt,2)), 'LineWidth', 3), hold on;
semilogy(zlistdetection, cumsum(sum(RdetectionsByRedshiftMt(:,Mtlist<=5),2)), 'LineWidth', 1);
semilogy(zlistdetection, cumsum(sum(RdetectionsByRedshiftMt(:,Mtlist>5 & Mtlist<20),2)), 'LineWidth', 1);
semilogy(zlistdetection, cumsum(sum(RdetectionsByRedshiftMt(:,Mtlist>=20),2)), 'LineWidth', 1); hold off;
legend('Total rate', 'From M_t<=5 M_o', 'From 5<M_t/M_o<20', 'From M_t>=20 M_o'),
set(gca, 'FontSize', 20); %for labels
xlabel('z'),
Expand Down
Loading