\section{Introduction}

Applications of atmospheric cosmic rays (CR) have grown in numbers in the last decade in the field of muography. Measurements of the muons flux attenuation or deviation have been successfully applied to the imaging or monitoring of large geophysical, archaeological or industrial structures. Among all charged particles reaching ground level, muons are the most numerous. Their energy loss or scattering when crossing a given length of matter is exploited to reconstruct the density of the medium. The measurement of the muon energy is usually not possible in small and compact field detectors, which are standard trackers using particle physics techniques (scintillators, resistive plate chambers, micro-megas etc).
The simulation of the muon flux follows in general two different approaches :
\begin{itemize}
    \item \textit{Analytical models} which provide semi-empirical formula adjusted on experimental datasets. See for example~: \textbf{Tang} \cite{tang_muon_2006}, \textbf{Shukla} \citep{shukla2018energy}, \textbf{Honda} \citep{Honda_2006qj}, \textbf{Gaisser} \citep{gaisser1990cosmic} ...
    \item \textit{Cosmic rays shower generators} which simulate extended air showers (EAS) from the primary cosmic rays down to the particles at ground. The most popular are : \textbf{PARMA} \citep{sato2015analytical}, \textbf{CRY} \citep{hagmann2007monte}, \textbf{MUPAGE} \citep{carminati2008atmospheric}, \textbf{MCEq}  \citep{fedynitch_mceq_2016}, \textbf{CORSIKA}  \citep{heck1998corsika} ...
\end{itemize}

The sensitivity of the technique relies on the model's accuracy, which should take into account the experimental conditions in the most realistic possible way. For instance, it has been shown in Jourde et al. \citep{jourde2016monitoring} and Tramontini et al. \citep{tramontini2019middle}, that atmospheric conditions (pressure and temperature) are strongly correlated with the muon flux. It is also clear that the geomagnetic field may play a significant role by deflecting charged particles towards the poles \citep{clay1928penetrating}, which leads to a decrease of the flux at the equator and an increase at high latitudes.\\

 This article details the methodology of a flux simulation based on CORSIKA, a very flexible generator which allows to change the atmosphere and geomagnetic field profiles in a given place and time. We configured the hadronic interaction models and necessary options to obtain an accurate flux fitting our experimental data. We also show how the the obtained model compares to various analytical models. The last sections present numerical results on the effects of geomagnetic field and atmospheric conditions (pressure, temperature ..) on the muon flux. A global conclusion is drawn on the obtained tool performance.

\section{Simulation strategy and validation}

\subsection{CORSIKA's parameters}
CORSIKA (COsmic Ray SImulations for KAscade) is a Monte Carlo code for simulating atmospheric particle showers initiated by high energy cosmic ray particles. Primary particles (protons, light nuclei ...) are tracked in the atmosphere until they interact, decay or are absorbed. All secondary particles are explicitly followed along their trajectories. Their parameters are stored when they reach an observational level. For more details on the physics involved in atmospheric showers processes, see reference \cite{heck1998corsika}. In this study, we have selected CORSIKA's options that are believed to be essential to simulate the muon flux in a specific way and we detail some of them below.

\subsubsection{Hadronic interaction models and primary particles trajectories}
The primary CR flux  is composed of several types of particles (H, He, C, O, Fe ...). When a primary particle reaches the top of the atmosphere, it undergoes hadronic interactions leading to the production of secondary CR. Among those particles, pions and kaons decay into muons. Different types of primary particle interaction models are available on CORSIKA. For the hadronic interaction models, we chose FLUKA for the low-energy interactions and QGSJET-II-4 for high energies, the best candidates in their energy domains.\\

The energy range of the primary particle is chosen to match the muon energy measurable in our tomography experiments~: $10^E$ with E=[1, 7]. This total range is split into several parts with a defined step. The number of simulated primary particles is weighted in each bin with an empiric law for the number of showers~: $N=10^{9-E}$. This choice is also motivated by issues on computational time which increases considerably with the energy of the primary particle. The energy spectrum of primary CR follows a roughly exponential law : $E^{-\gamma}$. $\gamma$ is calculated for each intermediate energy range with a fit defined by an analytical model describing the primary spectrum. We select the Papini et al. model \citep{papini1996estimate} based on real data fits. This reference allows to take into account the solar modulation on the expected primary CR. The total flux is the sum of the fluxes of the various constituents of primary CR. It is possible to run simulations for each primary particle type or to apply the \textit{superposition model} explained by Spurio \citep{spurio2015cosmic} on the muon fluxes (results presented here). On top of these energy and momentum parameters, the trajectories of particles are defined by their zenith and azimuth angles, ranging from 0° to 90 ° and from -180° to 180°, respectively. For zenith angles higher than 60° the curvature of the atmosphere is taken into account, it cannot be neglected.

\subsubsection{Earth's magnetic field}
The Earth is protected by a magnetic shield created by the Earth's magnetosphere, which reduces the intensity of the high-energy flux reaching the ground. The geomagnetic field (B) modifies the spectrum of particles bombarding our atmosphere. This presents a low-energy cutoff, where the Earth's magnetic field is able to deflect primary CR below 10 GeV near the equator and close to 1 GeV at higher latitudes. The primary CR intensity also varies with longitude because of the asymmetry of the geomagnetic axis with respect to the Earth's rotation axis \citep{pethuraj2019measurement}. Those "East-West" fluxes show differences in energy intensity up to 100 GeV. This difference is more marked at high altitude than on the ground. Finally, there are significant local variations of the geomagnetic field, which affect the intensity of CR, the most famous being the South Atlantic Anomaly (SAA). All these effects can be taken into account in our simulation.\\

For each place, we declare the horizontal ($B_x$) and vertical ($B_z$) components of the Earth's magnetic field (in~$\mu$T), $B_x$ being the magnetic north. They are generated by \textit{NOAA geomagnetic calculator} according to the reference \citep{tapia_first_nodate}. CORSIKA computes the total magnetic field and its inclination from these two components. The main magnetic field effects are mostly latitude-dependent and therefore related to $B_x$, as we detail in section \ref{sectionB}.

\subsubsection{Atmosphere properties}

The atmosphere can be divided into five layers: the troposphere, stratosphere, mesosphere, thermosphere and ionosphere. The atmospheric profiles presented in section \ref{section4} stop at the end of the stratosphere ( $\sim$ 50 km). The troposphere is the part of the Earth's atmosphere located between the surface and an altitude of about 8 to 15 kilometers, depending on latitude and season. It is thicker at the equator than at the poles. This layer concentrates three quarters of the atmospheric mass and the temperature decreases rapidly with altitude. The stratosphere extends, on average, between 12 and 50 km. It is characterized by an increase in temperature with altitude. The stratosphere begins at low altitude near the poles, because the temperature is lower there. The distribution of atmospheric density is therefore different at opposite latitudes.\\

Muons are produced at typical 10-15 km altitude (troposphere/stratosphere boundary). Their abundance is affected by the density differences in the atmosphere either by direct  re-interaction or by modification of their parent mesons survival probabilities before decay \citep{gaisser2016cosmic, grashorn2010atmospheric}. The effect is more important for high-energy muons, which result from high-energy mesons with larger lifetime due to time dilation and therefore with longer paths in the atmosphere. Thus high-energy muons are more sensitive to temperature changes. An input to CORSIKA is therefore the atmosphere's state in which the CR ray showers are generated. The state of the atmosphere is described by the density of the air at each altitude level. This one is calculated by converting the relative humidity into saturation vapor pressure with the Magnus formula \citep{abreu2012description}. We computed the parameters and altitudes of the layer boundaries from ERA5 data, the latest climate reanalysis produced by the ECMWF which combines large amounts of meteorological observations with estimates made from advanced modeling and data assimilation systems. Some atmospheric density profiles built are represented on Figure \ref{subplot_tout_T_density}.

\subsection{Differential fluxes from analytical models Vs CORSIKA}
The fluxes simulated with CORSIKA, along the procedures described before, are compared to analytical fluxes to check their relevance. For this purpose, we plot differential fluxes as a function of muon energy only, for a muon zenith angle equal to 0°, on Figure \ref{corsikadiff}. At low energies they disagree : the models fit on several orders of magnitude to end up separated at high energies. But we know that low-energy muons are important for calibration and high-energy muons for conducting tomography experiments. Analytical models are known to be poorly adapted to small and large energies, because few measurements are available for their fitting equations. The CORSIKA model instead is probably more reliable over the whole energy range. Furthermore, analytical models are not extrapolated for all zenith angles and they don’t take into account geodesics parameters, a limitation overcome by the CORSIKA approach.
 \begin{figure}[h]
    \centering
  \includegraphics[scale=0.32,trim=1.4cm 0cm 1.5cm 0cm,clip=true]{compflux_0.png}
    \caption{Differential fluxes as function of muon energy for a given zenith angle $\theta$=0°. Comparison between CORSIKA simulation and analytical models. }
    \label{corsikadiff}
\end{figure}

\subsection{Comparison of a CORSIKA flux with real data  \label{corsikavsdata}}
The best validation cross-checks for any simulation is the comparison to real data. Data presented here were taken in Lyon (France, latitude ~45° and close to sea level), in open sky conditions, with a 3 planes muon tracker (so-called muon telescope). We tilted the telescope progressively by step of $\ang{15}$ from the zenith ($\theta=\ang{0}$) to the horizontal ($\theta=\ang{90}$) directions. The muon flux is simulated in Lyon (France) respecting the geodesic constraints. Figure \ref{crossvalidation} displays the data/simulation comparison. Experimental points are the small crosses (1 tilt = 1 color) and data from CORSIKA's simulations are represented by yellow circles. Different fits are made on these two datasets with a simple $cos^2(\theta)$ at first order. Despite the still pending disagreement at large zenith angles, we observe a real improvement in the data/model comparison with respect to the usual $cos^2(\theta)$ (or analytical) fit. The simulated flux under-evaluates slightly the measured one for high zenith angles. To overcome this, it is necessary to constrain the simulated flux to the data to be as close as possible to reality.
\begin{figure}[h]
    \centering
  \includegraphics[scale=0.26,trim=1.4cm 0cm 1.5cm 0cm,clip=true]{compdata_corsikaPROCEEDING_Lyon010819_109.png}
    \caption{Integrated muon fluxes plotted as a function of the zenith angle $\theta$. The colored crosses represent the measured data for different inclinations of the detector. The pink crosses are a fit in $cos^2(\theta)$ on the total measured flux. The circles represent the simulated flux with the CORSIKA. A $cos^2(\theta)$ fit on it is drawn with a solid black line.}
    \label{crossvalidation}
\end{figure}

\section{Numerical Results: effects of the geomagnetic field}\label{sectionB}
In this section, we used our CORSIKA model, validated on a large dataset to simulate muon fluxes (E, $\theta$) for various
Earth magnetic field. We wanted to quantify those effects on the measurable muon flux in open sky. For this purpose, the atmospheric and altitude parameters have been fixed to constant values.

\subsection{$B_x$ influence}
We have set the vertical component as a constant ($B_z$= 20 $\mu$T) and the horizontal component ($B_x$ = 15 $\mu$T) for a simulation of reference. We then compute another simulation with parameters set at $B_x$ = 45 $\mu$T, $B_z$ = 20 $\mu$T. We simulate the muon flux with those two configurations and compute their ratio for four different energy ranges and different zenith angles. Figure \ref{$B_x$15BZ20$B_x$45BZ20_1GeV} (top), on left panel, presents the normalized intensity distribution of the flux ratio over the zenith angle, for each individual energy range.
It shows that when considering a larger energy range, the intensity ratios distribution tends to a narrow peak centered at 0.94834 corresponding to an effect of about 5$\%$. As expected the ratio tends to one for energy ranges greater than 10 GeV but affects the low-energy particles which are more deflected towards the poles. Right panel shows the flux ratio for different zenith angles $\theta$ for the same energy ranges. At low energies (1 to 10 GeV in purple) the flux is higher for $B_x$= 45 $\mu$T until 60° and between 85° and 90°, and lower between 60 and 85°. This probably arises from the fact that high-angles particles cross a larger section of the atmosphere and may be less deflected. For higher energies, the flux ratio remains more constant at each angle (in light blue, green and yellow). The total integrated flux ratio (in black) follows the behaviour of the "1 to 10 GeV" integrated flux.\\
\begin{figure}[h]
    \centering
\includegraphics[scale=0.33,trim=1.5cm 0cm 1.5cm 0cm,clip=true]{BX15BZ20BX45BZ20_12.png}
\includegraphics[scale=0.33,trim=1.5cm 0cm 1.5cm 0cm,clip=true]{BX15BZ20BX20BZ15_12.png}
    \caption{Comparison of two simulations, with same atmosphere and altitude, and different magnetic fields. (Top) The magnetic parameters are (1) $B_x$ = 15 $\mu$T, $B_z$ = 20 $\mu$T and (2) $B_x$ = 45  $\mu$T, $B_z$ = 20 $\mu$T. (Bottom) The magnetic parameters are : (1) $B_x$ = 15  $\mu$T  $B_z$ = 20  $\mu$T and (2) $B_x$ = 20  $\mu$T  $B_z$ = 15  $\mu$T.
    (Left panel) Normalized intensity distribution for different energy ranges from 1 to 10$^4$ GeV (in color) and for the whole energy range (in black). (Right panel) Flux ratio dependence on the zenith angle $\theta$ between two different magnetic field fluxes for the same energy ranges as in the left panel.}
    \label{$B_x$15BZ20$B_x$45BZ20_1GeV}
\end{figure}

\subsection{$B_z$ influence}

The same procedure is followed this time with a fixed horizontal component ($B_x$ = 30  $\mu$T) and various horizontal components ($B_z$~=~20 and 45  $\mu$T). We simulate the muon flux with these two configurations and we compute the ratio as a function of the energy and of the zenith angle. The results are not shown in this issue because only a slight difference at high zenith angles is observed. The $B_z$ component does not seem to affect the muon flux. \\

\subsection{$B$ intensity influence}

Here we have kept the previous simulation of reference, and compute another simulation with where the $B_x$ and $B_z$ components values have been switched. The intensity for the total geomagnetic field has been kept constant for both simulations (Figure \ref{$B_x$15BZ20$B_x$45BZ20_1GeV} (bottom)).
In this case the conclusions are the same as for Figure  \ref{$B_x$15BZ20$B_x$45BZ20_1GeV}. Low-energy muons are more deflected when $B_x$ increases. It is precisely this component that influences the muon flux because in this subsection the magnetic field intensity remained constant in both fluxes compared. \\

These numerical results validate that the magnetic field particularly affects low-energy muons which are deflected towards the poles. The very energetic mesons are less affected by the geomagnetic field. Low-zenith angle muons are less affected because the trajectory is short, at high zenith angle it is less clear because they are less and less energetic as they progress.\\

All the results presented in this section are subject to significant uncertainty. It is statistical and increases when the energy increases. Indeed, we simulate much less high-energy muons and extreme-angles muons. We have limits imposed by computational time and linked to the randomness of CORSIKA.

\section{Numerical results : Atmospheric parameters influence} \label{section4}

We desire to quantify the impact of atmospheric density on the flux, which is ultimately controlled by the temperature and pressure state of the atmosphere. Hence, at the same location, seasonal variations affect the muon flux over time.
For our tests we chose the city of Lyon in France (lat = 45.75, long = 4.75, $B_x$ = 22.71  $\mu$T, $B_z$ = 40.96  $\mu$T), place where we performed the open sky measurements shown in Figure \ref{crossvalidation}
The atmosphere density parameters used to simulate muon fluxes with CORSIKA are determined with "era5tool" and ERA5 datasets, for two different dates.

\subsection{Atmosphere models}
Temperature and density profiles in Lyon, during winter (12/30/20) and summer (08/01/19), are displayed on Figure \ref{subplot_tout_T_density}. These temperature and density profiles highlight which part of the atmosphere may affect the muon flux production and filtering. We observe that the density of the atmosphere decreases with altitude, and that colder atmospheres are denser especially at lower altitude.

\begin{figure}[H]
    \centering
  \includegraphics[scale=0.33,trim=1.5cm 0cm 1.5cm 0cm,clip=true]{suplot_lyon_T_densitu.png}
    \caption{(Top) Temperature and (Bottom) density profiles for altitudes from 0 to 50 km in Lyon during winter (30/12/20) in blue and in summer (08/01/19) in red.}
    \label{subplot_tout_T_density}
\end{figure}

\subsection{Atmosphere flux comparison : Seasonal effects}

Figure \ref{lyon010819-301220_intflux} presents the normalized intensity distribution for different energy ranges : 1-10 GeV, 10-100 GeV,  100-10$^3$ GeV, 10$^3$-10$^4$ GeV and the distribution for the whole energy range. It shows that when considering a larger energy range, the intensity distribution tends to a narrow peak centered at 0.91974 which means that the flux is 10 $\%$ higher in winter in Lyon (France). This effect is quite sizeable and must be properly accounted for, when precise muography is required over for low opacity targets. This effect depends of course on the particles energy. Right panel shows the flux ratio for different zenith angles $\theta$ in the range 0 to 90°, when considering the same energy ranges. It shows that in summer the flux is higher for high energies (100 to 10$^4$ GeV, in green and yellow), and lower in winter. For lower energies, the flux ratio is higher in winter (1 to 100 GeV, in light and dark blue), and lower in summer. These effects increase with the zenith angles $\theta$. Statistical errors of fluxes simulations are present on Figure \ref{lyon010819-301220_intflux} as in Section \ref{sectionB}. \\

 The decrease in atmospheric density increases the muon flux. We have seen before that the effect is more important for high-energy muons, which result from high-energy mesons with larger lifetime due to time dilation and therefore with longer paths in the atmosphere. We observe this effect on Figure \ref{lyon010819-301220_intflux} where the flux increases by 1$\%$ during summer.

 \begin{figure}[H]
    \centering
  \includegraphics[scale=0.33,trim=1.5cm 0cm 1.5cm 0cm,clip=true]{lyon010819Lyon__0m_301220_12.png}
    \caption{(Left) Normalized intensity distribution for several energies range from 1 to 10 000 GeV (in color) and for the whole energy range (in black). (Right) Flux ratio as a function of the zenith angle theta in summer and winter conditions in Lyon, France for the same energy ranges as left panel. The different atmosphere conditions are (1) 08/01/2019, (2) 12/30/20, with constant geomagnetic field and altitude.}
    \label{lyon010819-301220_intflux}
\end{figure}

\section{Conclusion}

In this paper, we present a muon flux simulation workflow accounting for muon-atmosphere interactions, based on the CORSIKA framework. We detailed our simulation strategy and the various relevant inputs from the hadronic interactions models to the atmosphere conditions. In particular, we used meteorological ERA5 pressure and temperature datasets to compute the required atmospheric density profiles. The workflow has been cross-validated against experimental evidences and standard semi-empiric models found in the literature. Simulations prove themselves to be a powerful tools to study and make predictions on tiny effects induced by the geomagnetic field or the atmospheric seasonal variations. Those effects are of increasing importance when one wants to produce muon imagery, and/or, long-term internal state surveys, on both ends of the opacity spectrum. Low-opacity targets imagery is controlled by low-energy muons filtered out by the density of the atmosphere. On the other side, high opacity targets imagery is largely affected by the process at stake for high-energy muon production. This study opens the gate to develop semi-empiric formulas predicting the evolution of the muon energy spectrum for each zenith angle, in relation with the atmospheric state. These formulas will be useful to correct recorded muon fluxes on the fly, when direct open-sky measurements are not available or not sufficiently refined in terms of energy description.

\section*{Acknowledgements}
This work was the subject of a CIFRE agreement between the ArcelorMittal Maizières Research SA and IP2I (Lyon).\\

\vspace{2cm}
\color{black}
\noindent\rule[0.25\baselineskip]{\textwidth}{1pt}


\section{INTRODUCTION}
\dl{The Gravitational wave (GW) catalogs of signals from the compact binary coalescence (CBC), present in the first half of the third observing run (O3a), had been recently released \cite{gwtc2,gwtc2_1} \footnote{The third observing run (O3) of the Advanced LIGO \cite{LIGO-det} and Advanced Virgo \cite{Virgo-det} (April 1, 2019, to March 27, 2020) detectors was divided into two segments as O3a run  (April 1, 2019 - October 1, 2019,)  and O3b (November 1, 2019 - March 27, 2020). } This includes intermediate black hole (IMBH) system GW190521 \cite{GW190521}, the merger of two compact objects with the unequal mass ratio GW190814 \cite{GW190814}, the inspiral of a binary neutron star (BNS) system, GW190425 \cite{GW190425}. Further,  the first confident observation of GW signal from a neutron star-black hole (NSBH) binary coalescence was made during O3b run \cite{NSBH,gwtc3}.}

Besides CBC signals, the current ground-based detectors are also sensitive to detecting GW short-duration bursts, such as signals from the core-collapse supernova (CCSN), cosmic strings, non-linear memory effects and isolated Neutron stars \cite{ccsn_ligo,cosmic_string,memory,magnetar,ns_glitch}. The LIGO-Virgo-Kagra collaboration has been actively searching for GW bursts signals since the fifth science run.  The all-sky burst search for short-duration transient  includes the burst signal duration of up to a few seconds \cite{allskyfith,allskysixth,allakyo1,alskyo2,allskyo3}. This search is identified as the unmodeled  search, sensitive to detect GW signals of different morphologies, including GW signals from the CBC system.

\dl{The morphology independent coherent WaveBurst (cWB) short-duration algorithm identifies triggers with short-duration excess power in the time-frequency domain \cite{cwb_2008}. The cWB made a significant contribution in the detection of the first GW signal from the binary black hole (BBH) merger GW150914  \cite{gw150914_cwb}  and the intermediate-mass black hole (IMBH) merger event \cite{gw190425_cwb}.  In cWB,  each trigger is associated with multiple attributes. The cWB  adopts an empirical classification of noise-based glitches based on various threshold values applied to different attributes.  The threshold values, though \textit{ad-hoc},  are chosen based on the simulation exercise.  These threshold choices vary based on the detector sensitivity during each run, the network combination, and the type of the transient signal.}

Machine learning (ML) techniques are the possible option to address such problems as they offer powerful tools for classification between signal and noisy transients \cite{noise_cl1,noise_cl2,noise_cl3,noise_cl4,noise_cl5,noise_cl6,noise_cl7,noise_cl8,noise_cl9,noise_cl10}. More generically, ML ideas have the notable potential to improve the detection sensitivity of unmodeled GW signals \cite{gmm,cwb_ml1,cwb_ml2,Mishra2021, Cuoco_2020}. ML techniques have also been used to enhance searches and classification of GW signals from CCSN \cite{ml_ccsn}.

\dl{In an earlier work \cite{gmm},  the authors had proposed an alternative ML based approach to address the \textit{ad-hoc} thresholding on the attribute set. In this Gaussian mixture model (GMM)-based
supervised machine learning approach, we modeled the cWB trigger attributes in the multidimensional attribute space. Thus, it provides an alternative approach to  thresholds applied on multidimensional space to the detection problem under the scalar log-likelihood ratio. In this work, we continue to use the Gaussian mixture models in obtaining the model for the attributes in the multidimensional parameter set.  Here, we explore this approach by carefully choosing the attribute set  and investigating the dependence on the search sensitivity with the appropriate usage of the re-parametrization approach. We present the search sensitivity for the generic short-duration gravitational-wave transients using the data from the O3a run with the low-frequency (16–1024 Hz) analysis.}

\dl{The paper is organized as follows. Section \ref{Sec:methodology} describes the unmodeled search algorithm cWB and Gaussian mixture modeling for cWB triggers. Section \ref{Sec:analysis} discusses the dataset used to study the GMM approach to all-sky search, Gaussian mixture model generation, and sensitivity improvement to generic signal morphologies and core-collapse supernova waveforms. Section  \ref{Sec:results} discusses the GMM method results on cWB triggers of coincident events from the O3a run of Advanced LIGO detectors. Finally, in Section \ref{Sec:discussions}, we summarize the significance of the GMM method on an all-sky search for short-duration transient signals using a minimal model approach.}

\section{METHODOLOGY}

\label{Sec:methodology}

\subsection{Coherent WaveBurst algorithm}

The cWB algorithm is a morphology-independent,  multidetector GW signal detection algorithm that coherently maps the multidetector's data into the multiresolution time-frequency scale domain using the Wilson-Daubechiers-Meyer wavelet transformation \cite{wdm}. It is based on the constraint maximum-likelihood ratio approach applied to the strain data in the time-frequency domain. The clusters of excess energy pixels are selected above the noise level of the detector and labeled as triggers if they exceed the thresholds on coherent energy ($E_c$) and network correlation coefficient ($c_c$) \cite{cwb_Klimenko,cwb_2021}. For each trigger, the cWB computes a collection of attributes that characterize signal as well as noisy transients properties.

We consider the following attribute sets for this study,  which are generic for the short-duration transient signals. The network coherent energy $E_{c}$, the effective correlated signal-to-noise $\eta_{c}$, the network correlation coefficients  $c_{c0}$ and $c_{c2}$, network energy short-duration $N_{ED}$, the ratio between the reconstructed energy and the total energy  $N_{\mathrm{norm}}$, the residual noise energy measure  $\chi^{2}$ and attributes pertaining to the noise vetoes like $Q_{\mathrm{veto}}$ and  $L_{\mathrm{veto}}$.  We list this attribute set  in Table \ref{tab:attributes} \cite{cwb_page}.  In standard cWB-based analysis, the threshold on the multidimensional attribute space is placed to distinguish the noise based triggers from the astrophysical GW event.

\begin{table*}[htp]
\begin{tabular}{c c
\hline
  \hline
 Attribute & Definition \\ [1ex]
 \hline

$E_{c}$ &   Cross correlation of reconstructed waveforms  between the detector pairs.\\
$\eta_{c}$ &   $\eta_{c} = \sqrt{\frac{E_{c}}{(K-1)}}$, where $K$ is the number of detector.\\
$c_{c0}$ & $c_{c0}=\frac{E_{c}}{(|E_{c}|+E_{n})}$, where $E_{n}$  is the energy of the residual noise.\\
$c_{c2}$ & $c_{c2}=\frac{E_{c} \times c_{c0}}{(|E_{c}|+E_{n})}$\\
$N_{ED}$ & Energy disbalance of the trigger between the detectors.\\
$N_{\mathrm{norm}}$ & Ratio between the reconstructed energy and the total energy. \\
$\chi^{2}$ &  $\chi^{2} = \frac{E_{n}}{N}$,  where N is the number of independent wavelet amplitude of the event. \\
$Q_{\mathrm{veto0}}$ & Energy distribution of an event over different time segments. \\
$Q_{\mathrm{veto1}}$ & An estimate of quality factor assuming the signal to be a CosGaussian  (Q factor).\\
$L_{\mathrm{veto0}}$ & Central frequency of the reconstructed signal, to identify narrow band glitches.\\
$L_{\mathrm{veto1}}$ & Root mean square frequency of the reconstructed signal.\\
$L_{\mathrm{veto2}}$ & Energy ratio between pixel energy and total energy of the event.\\[1ex]
\hline

\end{tabular}
\caption{Summary of the selected attributes associated with the short-duration burst transients in cWB and used for the GMM-based analysis.}
\label{tab:attributes}
\end{table*}

\subsection{Gaussian mixture modeling in postproduction }
 \dl{In \cite{gmm}, we proposed the Gaussian mixture modeling approach \cite{bishop:2006} to construct two distinct GMM models in the multidimensional trigger attribute set for astrophysical GW signals and noise glitches. We  used the Bayesian information criterion (BIC) method to obtain the optimum number of  Gaussians and the GMM model parameters (mean, covariance and weights).  For each model (signal and noise), we had obtained the maximum log-likelihood statistics $W = \ln(\hat{\mathcal{L}} )|_{\hat{K}}$, where $\hat{\mathcal{L}}$ is the maximum value of the likelihood function with the optimum number of Gaussians $\hat{K}$ \cite{gmm,bic_78,bic}. Then, we defined the detection statistic, $T$,  as
\begin{align}
\label{T}
T = W_{\text{s}}- W_{\text{n}}\,.
\end{align}
where the subscript ${\mathrm s}$ and ${\mathrm n}$ stand for signal and noise. We had used this quantity as a ranking statistic for each cWB trigger instead of applying threshold values on the multidimensional  attribute set.  We suggest the reader to follow \cite{gmm} for more details.}

\dl{A well behaved multidimensional Gaussian distribution is preferred to obtain GMM. When the distribution of a selected subset of attributes is not well behaved (highly peaked), the GMM overestimates the required number of Gaussians. Thus, we either precondition the data by reparametrizing some of the cWB attributes or combining two attributes in a single attribute. }

We reparameterize the subset of attributes  $\eta_{c}$, $c_{c0}$, $c_{c2}$, $N_{ED}$, $E_{c}$, $Q_{\mathrm{veto}}$ and $L_{\mathrm{veto}}$, which shows the non-Gaussian distribution by either applying a logarithmic or inverse sigmoid transformation. Additionally, we define a new parameter $L_{\mathrm{ratio}}: L_{\mathrm{veto1}}/L_{\mathrm{veto0}}$ instead of using $L_{\mathrm{veto0}}$ and $L_{\mathrm{veto1}}$ as two different attributes. The distribution of attributes after the reparameterization is well behaved and better conditioned than before reparametrization.  This shows a direct implication on the number of Gaussians in the GMM model.  For example,  for the data set considered here,  the number of optimum Gaussians in the original attribute set for the signal and noise model is 113 and 115, whereas the optimum number of Gaussians is reduced to 90 and 82 for the reparametrized attribute set.  More details about the dataset is in Sec.  \ref{Sec:analysis}.
More details on reparametrization can be found in Appendix.

\section{O3\lowercase{a} ANALYSIS WITH GMM}

\label{Sec:analysis}

\dl{This study uses publicly available O3a data \cite{GWOSC_O3a,GWOSC} from the LIGO Hanford and LIGO Livingston (HL) network in the low-frequency range (16-1024 Hz)\footnote{As we focus on maximization of the detection efficiency,  we consider the HL network for this study rather than HLV because of the similar reasoning mentioned in \cite{allskyo3}.}. The total amount of coincidence data between the two detectors used for this analysis is 104.9 days. The distribution of accidental triggers is calculated by time-shifting the data of one LIGO detector with respect to the other LIGO detector by an amount that breaks any correlation between the detectors for an actual signal. For the cWB O3a analysis, 1000 years of background data were generated.}

\begin{table}[htp]
\begin{tabular}{c c c c
\hline
  \hline
 Attribute & LF1 & LF2 & LF3 \\ [1ex]
 \hline
$c_{c0}$ & > 0.8 &> 0.8 &> 0.8\\
$c_{c2}$ & > 0.8 &> 0.8 & > 0.8 \\

$N_{\mathrm{norm}}$ & > 2.5 & > 2.5 & > 2.5\\

$Q_{\mathrm{veto0}}$ & = 0 & $\neq  0$ & \\
$Q_{\mathrm{veto1}}$ & $\leq 3$ & $\leq 3$ & > 3 \\
\hline

\end{tabular}
\caption{Definitions of the search bins used for the O3a  low-frequency short-duration burst analysis. The different thresholds are applied to attributes $c_{c0}$, $c_{c2}$, $N_{\mathrm{norm}}$, $Q_{\mathrm{veto0}}$ and  $Q_{\mathrm{veto1}}$ to classify the triggers into different bins based on background trigger distribution.}
\label{tab:bin_thershold}
\end{table}

\subsection{The standard O3a cWB all-sky analysis}

\dl{In standard cWB O3a analysis, triggers are required to pass frequencies above 24 Hz and high network correlation with $c_c$ above 0.8 for the HL detectors. These triggers are classified into three different bins  as LF1, LF2, and LF3 based on background trigger morphologies to isolate the triggers to a small part of the parameter space  \cite{allskyo3}.  LF1 contains a population of glitches dominated by very short and loud blip-type glitches, with negligible energy outside the single oscillation \cite{blip_glitch}. LF2 contains triggers resembling  blip glitches with Q factor $\leq$ 3, and LF3 contains remaining high Q factor triggers. Table \ref{tab:bin_thershold} lists the thresholds applied to cWB attributes for classification into different bins used in O3a analysis. After applying thresholds on the cWB attributes, the events are ranked based on inverse false alarm rate (iFAR). A trials factor of three is applied to iFAR values corresponding to three different background bins, and a threshold of 100 years on iFAR is used for the significant detection \cite{allskyo3}. By contrast, GMM based postprocessing does not require the creation of multiple analysis bins (see below).}

\subsection{GMM model generation}

For O3a GMM analysis, we choose most of the attributes required to characterize the generic short-duration transient signal from standard cWB (Table \ref{tab:attributes}) to generate the GMM model.
We exclude the attributes reconstructed signal's central frequency, duration, and strain from this analysis, as these  attributes are strongly affected by the priors on the injection parameters for the foreground model  \cite{gmm,cwb_page}. Since we do not want the GMM to develop a bias towards the choice of the signal population used in the simulation.  For GMM analysis, we consider those triggers with $\eta_{c} > 5.5$ and $c_{c}>0.5$  for the HL network.  The cWB uses different thresholds on attributes to generate gravitational wave candidates and veto out the noise triggers in its standard configuration.

We consider $70\%$ of noise background as well as simulated burst signals for training and the rest for testing the GMM method. We use the equal percentage of generic waveforms (Gaussian Pulse, sine-Gaussian wavelets, and white noise burst) as the training set aimed to improve the detection efficiency of all injected waveforms. We construct the GMM for the noise and signal triggers using this simulation \cite{gmm}. More specifics of the models can be found in Appendix.

\subsection{Data simulation}

\dl{A set of generic short-duration burst waveforms (typically used in the all-sky short-duration bursts  search \cite{allskyo3}) are injected in the HL O3a noise to estimate the search sensitivity using the cWB plus GMM method. They are sine-Gaussian wavelets (SG), Gaussian pulses (GA), and band-limited white-noise bursts (WNB) signals injected over a range of amplitude expressed in terms of the root-sum-squared strain amplitude ($h_{\mathrm{rss}}$). The SG waveforms are characterized by mean frequency $f_{0}$ and quality factor $Q$,  the GA waveforms are described by the duration $\tau$, and the WNB signals are characterized by lower frequency bound $f_{\mathrm{low}}$,  frequency bandwidth $\Delta{f}$, and duration $\tau$.  The SG and GA signals are injected as $h_{\mathrm{rss}} = (\sqrt{3})^{N}5 \times 10^{-23} Hz^{-1/2}$ over a grid of maximum strain values with N ranges from 0 to 8,  and  WNB waveforms are injected uniformly as the square of the signal's distance. Linearly polarized signals are used  for GA, and sine-Gaussian wavelets use both elliptical (SGE) and linearly (SG) polarized waveforms \cite{allskyo3, allskysixth}. Elliptical waveforms are uniform in the cosine of the source inclination angle, which is defined as the angle between total angular momentum and the line of sight. The WNB represents isotropic emission at the source and carries an equal amplitude from both polarization of GW strain at the detector \cite{Sutton_wnb}. Table \ref{tab:waveforms} lists the simulated signals and their characteristic parameter values used in this analysis.}

\begin{table}[htb]
\begin{tabular}{c c c c}
\hline
\hline
\multicolumn{3}{c}{Sine-Gaussian Burst (SGW)}\\
\hline
No. &$f_0$ (Hz) & $Q$ &  \\
\hline
1 & 70 & 3 &  \\
2 &70 & 9 &  \\
3 &70 & 100 &  \\
4 &100 & 9 &  \\
5 &153 & 9 & \\
6 &235 & 3 &  \\
7 &235 & 9 &  \\
8 &235 & 100 &  \\
9 &361 & 9 &  \\
10 &554 & 9 & \\
11 &849 & 3 & \\
12 &849 & 9 & \\
13 & 849 & 100 & \\
\hline
\multicolumn{3}{c}{White-Noise Burst (WNB)}\\
\hline
 & $f_{\text{low}}$ (Hz) & $\Delta f$ (Hz) & $\tau$ (s)\\
\hline
14 & 150 & 100 & 0.1 \\
15 & 300 & 100 & 0.1 \\
16 & 750 & 100 & 0.1 \\

\hline
\multicolumn{3}{c}{Gaussian Pulse (GP)}\\
\hline
 &  & & $\tau$ (s) \\
\hline
17 & & & 0.1 \\
18 &  & & 1 \\
19 & & & 2.5 \\
20 & & & 4 \\
\hline
\hline
\end{tabular}
\caption{List of generic burst waveforms and their characteristic parameters.}\label{tab:waveforms}
\end{table}

\dl{The low-frequency,  generic all-sky short-duration burst search  targets the frequency range of GW burst from most of the core-collapse supernovae (CCSN). To demonstrate the  robustness of the cWB plus GMM approach against different signal morphologies, here,  we analyze the search sensitivity to CCSN  waveforms used in all-sky short-duration transient analyses during O3a  \cite{allskyo3}. The  waveforms are generated  from five different  3D CCSN simulations  models like $\mathrm{s18}$  \cite{s18},  $\mathrm{m20}$ \cite{m20},  $\mathrm{s9}$ \cite{s9}, $\mathrm{m39}$ \cite{m39}, and $\mathrm{35OC}$ \cite{35OC}. The waveforms are distributed uniformly in distance  with maximum distance for $\mathrm{s18}$, $\mathrm{m20}$, $\mathrm{s9}$, $\mathrm{m39}$, and  $\mathrm{35OC}$ CCSN models set to 25 kpc, 5kpc, 5kpc, 70kpc, and 70 kpc, respectively \cite{allskyo3}.}

\subsection{ Sensitivity improvement for generic signal morphologies}

\begin{figure*}[htp]

    {\includegraphics[width=1\linewidth,height=0.55\textheight]{roc_commonlabel.pdf}}
\caption{ROC curves for the cWB plus GMM search (solid lines)
compared to the standard cWB search (dotted lines) for the simulated waveforms. Panel (a) shows GA class waveforms with duration 0.1, 1, 2.5 and, 4 s. (b) shows elliptically polarized SG waveforms with a quality factor $Q=3$ and frequency at 70, 235, and 849 Hz. Panel (c) and (d) show the SG waveforms with $Q=9$ and 70-900 Hz frequency range. (e) shows SG with $Q=100$ and frequency at 70, 235, and 849 Hz. (e) WNB is characterized by $f_{\mathrm{low}}$ at 150, 300, and 750 Hz with frequency bandwidth and duration at 100 and 0.1 s.}
\label{fig:roc_curve}
\end{figure*}

We use the receiver-operating-characteristic (ROC) curve to characterize the sensitivity improvement. The ROC curve of the standard cWB and the cWB plus GMM relates the detection efficiency against the false alarm rate (FAR) is shown in Fig. \ref{fig:roc_curve}.  The detection efficiency is the fraction of the recovered injections by a given algorithm.  The standard cWB uses $\eta_c$, and cWB plus GMM uses T as a detection statistic.  We show that the cWB plus GMM algorithm significantly enhances the detection efficiency for all the 20 waveforms compared to the standard cWB.

In order to understand this improvement,  we show the two-dimensional distribution of attributes $c_{c0}$, $N_{\mathrm{norm}}$, $Q_{\mathrm{veto0}}$ and $Q_{\mathrm{veto1}}$ of the signal triggers with iFAR $\geq$ 100 years in cWB plus GMM (by purple color) and standard cWB (by gray color) analysis in Fig. \ref{Fig.2dhist}. From panels (a) and (c), it is evident that with GMM, we have improved the signal/noise classification ability that we can now recover events below $c_{c0}$ of 0.8.

The most significant improvement after using cWB plus GMM is shown in the recovery of GA morphology (see Fig. \ref{fig:roc_curve}(a)).  The GA signal mostly resembles blip glitches and lies mainly in bin LF1.  During the O3 run, the search sensitivity of standard cWB to the GA waveforms has worsened compared to its sensitivity O2.  This is primarily due to an abundance of blip glitches  that resemble GA signals. However, with GMM based postproduction approach appears to have modeled them well, which is also evident in the recovery of GA with $Q_{\mathrm{veto1}} <3 $ (in Fig \ref{Fig.2dhist} panels (b) and (d)) those otherwise would get classified as bin LF1 events in the standard cWB.  In standard cWB analysis, detection efficiency for triggers falls in each bin estimates separately. However, the signal waveforms in bin LF1 are mostly buried in background glitches.  As a result, cWB plus GMM has shown an enhanced sensitivity for Gaussian pulses with an improvement of greater than $20\%$ at an iFAR of 100 years. The signal waveforms lying in bin LF1 show comparatively lower significance to other morphologies in both cWB plus GMM and standard cWB method at the given signal strength because of excess glitches.

The cWB plus GMM approach marginally improves the SG class of injections  (Fig. \ref{fig:roc_curve}(b)-(e)). This indicates that our approach can recover other burst morphologies considerably better than standard cWB. For WNB,  the cWB plus GMM approach consistently improves the search sensitivity with respect to standard cWB (Fig. \ref{fig:roc_curve}(f)). WNB falls in the  bin LF3,  the cleanest part of the background distribution and cWB plus GMM has shown improvement of around $10\%$ throughout  the iFAR threshold.

It should also be noted that categorizing the events into various bins in standard cWB analysis iFAR is penalized by a trials factor of three. Whereas in cWB plus GMM based analysis, there is no such bin classification.  This can also be attributed to some improvement in detection efficiency.

\begin{figure*}[t!]
  \centering
    \scalebox{0.235}[0.19]{\includegraphics{2Dhist_sep_zoom.pdf}}
  \caption{The two-dimensional distribution of signal trigger attributes $c_{c0}$ \textit{vs} $N_{\mathrm{norm}}$ ((a) and (c)) and $Q_{\mathrm{veto1}}$ \textit{vs} $Q_{\mathrm{veto0}}$ ((b) and (d)) for triggers recovered by cWB plus GMM  (purple) and standard cWB (gray) at iFAR $\geq$ 100 years.  The inset shows a zoom-in of the 2D distribution of $Q_{\mathrm{veto1}}$ \textit{vs} $Q_{\mathrm{veto0}}$. The black dashed lines show the postproduction threshold applied to these attributes for the standard cWB all-sky search, which classify the events into various bins as shown in Table \ref{tab:bin_thershold}.}
  \label{Fig.2dhist}
\end{figure*}

\subsection{Robustness test to CCSN}

To test the robustness against the different morphologies of waveforms and distribution, we applied the GMM model on the core-collapse supernovae (CCSN) injections used in  \cite{allskyo3}  that are not included in the training data set. Fig. \ref{fig:roc_sn} shows the ROC curves for injections with different CCSN waveforms. The  models $\mathrm{35OC}$ and $\mathrm{m39}$  describe the explosion driven by rapid rotation of  massive progenitor star with  zero age main sequence (ZAMS) mass of $35 M_{\odot}$ and $39 M_{\odot}$, which creates strong GW emission \cite{allskyo3}. Both standard cWB and cWB plus GMM show comparable high efficiency to CCSN model $\mathrm{35OC}$ at low FAR because GW signal occurs frequency above 100 Hz, which is optimally sensitive to the cWB low-frequency search method. However, the CCSN model $\mathrm{m39}$,  cWB plus GMM enhances the detection efficiency compared to the standard cWB. Since  GW amplitude peaks at a frequency $\sim$ 750 Hz, at the edge of the postproduction threshold  applied for standard cWB. The  standard cWB and cWB plus GMM give comparatively low efficiency for all the other models such as $\mathrm{m20}$, $\mathrm{s18}$ and $\mathrm{s9}$, which describes GW emission from nonrotating progenitors with ZAMS mass ranges from $9 M_{\odot}$ to $20 M_{\odot}$. Among them,  cWB plus GMM performs better than standard cWB at all FAR threshold values.  The ML based cWM plus GMM method does not focus on one class but shows good improvement for a broad class of waveforms at detectable iFAR levels.
\begin{figure}[H]

    {\includegraphics[width=1\linewidth,height=0.193\textheight]{SN_roc.pdf}}
\caption{ROC curves for CCSN waveforms. As before, the cWB plus GMM search ROC curves are represented by solid lines, while the standard cWB search ROC curves are plotted as dotted lines.}
\label{fig:roc_sn}
\end{figure}

\section{SEARCH RESULTS}
\label{Sec:results}

\dl{Following the sensitivity study,  we applied the cWB plus GMM on the O3a data from LIGO Hanford and LIGO Livingston.
The maximum background for cWB plus GMM analyses is 200 years since we used 20\% of the total O3a background as test data.
The results of this search from the O3a low-frequency region are shown in Fig. \ref{fig:zerolag}.  We plot the cumulative number of recovered events by standard cWB (by green color) and cWB plus GMM (by purple color) in terms of their iFAR values.}

\dl{The cWB plus GMM method recovers all CBC events identified by the targeted cWB search for CBC in GWTC-2 and GWTC-2.1~\cite{gwtc2, gwtc2_1,3-OGC}. No new candidates were found with high statistical significance in this search. The most significant non-CBC event observed at UTC 2019-09-30 23:46:52 with an iFAR of 0.3 years in the cWB plus GMM method shows 0.01 years of iFAR in the standard cWB search. The second most significant event at UTC 2019-05-11 04:12:15 was with an iFAR of 0.15 years in cWB plus GMM, which had shown iFAR of 0.002 years in standard cWB.  We note that the two loudest non-CBC events in standard cWB all-sky search are at UTC 2019-09-28 02:11:45 and UTC 2019-08-04 08:35:43, with an iFAR of 0.53 years, and 0.19 years~\cite{allskyo3}, showed less significance in the cWB plus GMM search with an iFAR of 0.006 and 0.05, respectively.  After excluding the CBC signals, the distribution of coincident events is close to the predicted background for the analyzed time.}

Table~\ref{tab:cbc} lists the CBC events from GWTC-2 and GWTC-2.1 that are also observed by standard cWB all-sky burst search \cite{allskyo3}. We list their attribute values,  $T$ statistic value, total source mass,  and the corresponding iFAR values for both standard all-sky burst cWB search~\cite{allskyo3} and cWB plus GMM search~\footnote{The significance values for standard cWB shown in this table are different from those listed in GWTC-2 and GWTC-2.1 because the significance estimates in these publications are based on a version of cWB optimized for CBC events. Since the cWB plus GMM approach is applied in the all-sky search context, we find it fairer to compare the event significance obtained for all-sky searches.}. We observe that the cWB plus GMM method is less sensitive to low mass BBH systems (with total source mass < 60 $M_{\odot}$) since these events are more significant in the standard cWB all-sky search. The exception is GW190412,  which has significantly asymmetric component masses and shows the contribution from higher modes\footnote{The component masses for GW190412 system are $m_{1} = 30.1^{+4.6}_{-5.3} M_{\odot}$ and $m_{2}=8.3^{+1.6}_{-0.9} M_{\odot}$.} \cite{GW190412}.  The IMBH system GW190521 \cite{GW190521} of total mass $163 M_{\odot}$ was observed with high significance in cWB plus GMM with an iFAR of greater than 200 years, which shows an iFAR for 65 years in all-sky burst cWB search \cite{allskyo3} \footnote{GW190521 reported with an iFAR greater than 4900 years using targeted CBC searches \cite{GW190521}.}.  A similar trend  is observed for other events in Table~\ref{tab:cbc}, with total source mass $M > 100 M_{\odot}$. The CBC event count acts as a verification step just to review the sensitivity of cWB plus GMM search against standard cWB all-sky search.

 \begin{figure}[H]

    {\includegraphics[width=1.0\linewidth,height=0.2\textheight]{zerolag_zoom.pdf}}
\caption{The cWB plus GMM (purple) and standard cWB  (green) observed the cumulative number of events versus iFAR in the low-frequency range. The iFAR distribution for all coincident data is plotted using triangular markers, while the line with circular markers shows the distribution after the removal of known CBC events. The expected mean value of background shows as a solid gray line with 1, 2 and 3$\sigma$ Poisson uncertainty (shaded regions).}
\label{fig:zerolag}
\end{figure}

\begin{table*}[htp]
\centering
\scriptsize
\begin{tabular}{c c c c c c c c c c
\hline
 \hline
 & &  &   &  &  & & & \multicolumn{2}{c}{iFAR in years} \\
 \hline
 Event & $\eta_{c}$ & $Q_{\mathrm{veto}}$ & $c_{c0}$ & $\chi^{2}$ & $N_{\mathrm{norm}}$ & $T$ & $M(M_{\odot})$ & cWB+GMM & cWB\\ [1ex]
 \hline\hline
 GW190408\_181802 & 8.59 & 0.92 & 0.96 &0.13 & 5.09& -0.41& $43.0^{+4.2}_{-3.0}$  &0.30 &25.14 \\
 GW190412 & 11.69 & 4.16 & 0.95 &0.06 &5.4 & 13.21& $38.4^{+3.8}_{-3.7}$ &15.62&14.86  \\
 GW190421\_213856 & 6.46 & 0.31 & 0.97 & -0.07 & 4.41 &-0.38 & $72.9^{+13.4}_{-9.2}$ & 0.30&0.04 \\
 GW190426\_190642 & 5.52 & 0.45 & 0.88 & 0.08 & 4.07 & -4.85 & $184.4^{+41.7}_{-36.6}$ & 0.02 & 0.01 \\
 GW190503\_185404 & 7.34 & 0.34 & 0.93 &-0.02 &4.76 &1.65 & $71.7^{+9.4}_{-8.3}$ & 0.84& 0.70\\
 GW190513\_205428 & 7.05 & 1.67 & 0.86 & 0.15 & 3.77 &-2.99 & $53.9^{+8.6}_{-5.9}$ & 0.07& 0.28\\
 GW190517\_055101 & 6.08 & 0.19 & 0.88 & -0.15 & 3.05 & -2.79 & $63.5^{+9.6}_{-9.6}$ &0.08 &0.01 \\
 GW190519\_153544 & 10.13 & 0.53 & 0.89 & 0.01 & 7.63 &18.04 & $106.6^{+13.5}_{-14.8}$ & 33.83 &7.78 \\
 GW190521 & 9.24 & 0.60 & 0.92 &-0.16 & 10.53 &32.45 & $163.9^{+39.2}_{-23.5}$ & $>200 $ &65.38 \\
 GW190521\_074359 & 14.19 & 0.56 & 0.96 &-0.08 &8.44 &72.77 & $74.7^{+7.0}_{-4.8}$ & $>200$ & 326.88\\
 GW190602\_175927 & 7.25 & 0.43 & 0.95 & -0.13 &6.5 & 0.73 &  $116.3^{+19.0}_{-15.6}$ &0.54&0.51 \\
 GW190706\_222641 & 9.29 & 0.79 & 0.83 & -0.10 &7.36 &24.93 & $104.1^{+20.2}_{-13.9}$ & $>200$ &65.38 \\
 GW190727\_060333 & 5.86 & 0.35 & 0.96 &  0.17 &4.96 & -2.94 & $67.1^{+11.7}_{-8.0}$ &0.07 & 0.006\\
 GW190728\_064510 & 6.50 & 3.94 & 0.87 &-0.13 & 2.55& -4.93 & $20.6^{+4.5}_{-1.3}$ &0.02 &0.051 \\
 GW190828\_063405 & 10.27 & 0.84 & 0.82 &0.10 & 5.01&8.78 & $58.0^{+7.7}_{-4.8}$ & 7.52&163.44 \\
 GW190915\_235702 & 8.07 & 0.42 & 0.95 &0.06 &4.29 &5.29 & $59.9^{+7.5}_{-6.4}$ & 3.07&5.36 \\
 GW190929\_012149 & 5.97 &  0.22 &  0.85 &0.103 &3.44 &-6.20 & $104.3^{+34.9}_{-25.2}$ & 0.01&0.009 \\[1ex]
 \hline
\end{tabular}
\caption{Table of CBC candidate events identified by cWB in GWTC-2 and GWTC-2.1. The columns show the effective correlated signal-to-noise $\eta_{c}$, energy distribution of the event over different time segments  $Q_{\mathrm{veto0}}$, network correlation coefficients  $C_{c0}$, residual noise energy measure  $\chi^{2}$, the ratio between the reconstructed energy and the total energy  $N_{\mathrm{norm}}$  and source total mass $M$ reported in GWTC-2 \cite{gwtc2} and GWTC-2.1 \cite{gwtc2_1}. The two rightmost columns report event significances for the cWB plus GMM search and the standard cWB all-sky search.}
\label{tab:cbc}
\end{table*}

\section{DISCUSSIONS}
\label{Sec:discussions}

\dl{This paper provides a robust method of direct usage of a supervised machine learning approach at the postproduction stage of a well-established model-independent cWB detection algorithm. The work employs a Gaussian mixture modeling approach to model cWB triggers (signal as well as noise) in a multidimensional attribute set. The GMM modeling helps to classify background glitches from the transient GW signals in a log likelihood-based test statistic. While this approach was initially conceptualized in \cite{gmm} here, we extend it further in conditioning the attribute set by first taking an appropriate subset targeted for a given signal class. Then introduce appropriate reparametrization of the attributes ensuring the well-behaved distributions for the construction of GMM.}

\dl{We use the cWB triggers set from the all-sky short-duration burst search during O3a coincident data from the two Advanced LIGO detectors. We consider the signals from generic morphology and specific CCSN models used in \cite{allskyo3} for sensitivity study. For generic morphology, there is a definite improvement in sensitivity for a fixed value of false alarm rate for all waveforms using the cWB plus GMM method over the standard cWB. We observed notable improvement for the Gaussian pulse, which falls in the bin LF1 containing a very short and loud (blip-type) population of glitches. Sensitivity improvement to CCSN models clearly demonstrates the robustness of the cWB plus GMM method against the variation in morphology as well as distribution of waveforms as CCSN waveforms were not part of the training set. }

\dl{The application of cWB plus GMM on the coincident data search is consistent with the all-sky short-duration burst results using the standard cWB search as published in  \cite{allskyo3}.  No additional, significant events are observed.  The search recovered the GWs signals from the BBH merger as found by the targeted search for CBC signals.  We observed higher significance for massive CBC events compared to that reported in the all-sky burst \cite{allskyo3}. This clearly shows the promise of the cWB plus GMM approach at the postproduction stage of the burst algorithm as it gives a better handle in distinguishing the short-duration bursts from that of the short-duration noisy glitches. We plan to further extend this work to Generalized mixture modeling \cite{generalized}, which allows us to model more complex distributions for the fourth and fifth observing runs of ground-based GW detectors~\cite{o4_o5}.}

\begin{acknowledgments}
The authors thank Shubhanshu Tiwari for valuable discussion and suggestions. D.L. acknowledges support from Swiss National Science Foundation  (SNSF) grant number 200020-182047.  A.P. acknowledges the SERB Matrics grant MTR/2019/001096 and SERB-Power-fellowship grant SPF/2021/000036 of Department of Science and Technology,  India for support. The authors are thankful to LIGO-Virgo-Kagra collaboration for provisions of the cWB attribute data used for the all-sky short-duration O3a analysis. The authors are grateful for the computational resources and data provided by the LIGO Laboratory and supported by National Science Foundation Grants No.PHY-0757058 and No.PHY-0823459. C.M. and I.S.H. are supported by the Science and Technology Research Council (grant No. ST/V005634/1) and the European Cooperation in Science and Technology (COST) action CA17137. The open data is available in the Gravitational Wave Open Science Center \cite{gwopen}, a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. The authors also acknowledge the use of the LDG clusters for computational/numerical work.
\end{acknowledgments}

\renewcommand{\thesubsection}{\Alph{subsection}}

\section*{APPENDIX: REPARAMETRIZATION OF ATTRIBUTES}
\label{app:reparameterization}

\begin{table}[!htp]
\caption{Details of reparametrization.}
\label{tab:repar}
\begin{subtable}{0.45\textwidth}{
\centering
\begin{tabular}{c c}
\hline
\hline
Original  attribute set & Reparametrized  attribute set \\ [1ex]
 \hline
$\eta_{c}$ & $log_{10}(\eta_{c})$\\
$c_{c0}$ & $logit(c_{c0})$\\
$c_{c2}$ & $logit(c_{c2})$\\
$N_{ED}$ & $log_{10}(N_{ED} + 10^{3})$\\
$E_{c}$ & $log_{10}(E_{c})$\\
 $N_{\mathrm{norm}}$ &  $N_{\mathrm{norm}}$\\
 $\chi^{2}$ &  $\chi^{2}$\\
 $Q_{\mathrm{veto0}}$ & $log_{10}(Q_{\mathrm{veto0}} + 1)$\\
 $Q_{\mathrm{veto1}}$ & $log_{10}(Q_{\mathrm{veto1}})$\\
 $L_{\mathrm{ratio}}$ & $logit(L_{\mathrm{ratio}})$\\
$L_{\mathrm{veto2}}$ & $logit(L_{\mathrm{veto2}} \times 0.99)$\\[1ex]

\end{tabular}}
\end{subtable}
\vspace{0.4in}
\begin{subtable}{0.45\textwidth}{
\centering
\begin{tabular}{c c c
  \hline
  \hline
  & \multicolumn{2}{c}{Number of optimum Gaussians in GMM} \\
  \hline
Model & Original attribute set & Reparametrized attribute set \\  [1ex]
 \hline
 Signal & 113 & 90\\
 Noise & 115 & 82 \\[1ex]
\hline
\end{tabular}}
\end{subtable}
\end{table}

Given the cWB trigger attributes $E_{c}$, $\eta_{c}$, $c_{c0}$, $c_{c2}$, $N_{ED}$,  $N_{\mathrm{norm}}$,  $\chi^{2}$, $Q_{\mathrm{veto0}}$ ,$Q_{\mathrm{veto1}}$, $L_{\mathrm{veto0}}$, $L_{\mathrm{veto1}}$, and  $L_{\mathrm{veto2}}$. We reparametrize  some of the attributes which do not follow a well behaved Gaussian distribution. In Table \ref{tab:repar}, we show the reparametrized attributes along with their original form. We further show the reduction in the number of Gaussians after the reparameterization of the attributes. The distribution of reparametrized cWB trigger attributes for the signal initially chosen for GMM analysis and after the reparameterization is shown in Fig. \ref{fig:corner_plot}.

\begin{figure*}
\centering
\includegraphics[width=0.9\textwidth,height=0.68\linewidth]{corner_worepar_signal.pdf}
\hspace{\fill}
\includegraphics[width=0.9\textwidth,height=0.68\linewidth]{corner_repar_signal.pdf}
\caption{The corner plot showing the one and two dimensional distribution of the reparametrized cWB trigger attribute for the signals before (top) and after (bottom) the reparameterization.}
\label{fig:corner_plot}
\end{figure*}

\clearpage

\bibliographystyle{apsrev4-1}

\subsection{Footnotes}%

\footnote{Automatically placing footnotes into the bibliography requires using BibTeX to compile the bibliography.};
author-year citation styles place the footnote at the bottom of the text column.
Note: due to the method used to place footnotes in the bibliography,
\emph{you must re-run Bib\TeX\ every time you change any of your document's footnotes}.

\begin{table}[b]
\caption{\label{tab:table4}%
Numbers in columns Three--Five are aligned with the ``d'' column specifier
(requires the \texttt{dcolumn} package).
Non-numeric entries (those entries without a ``.'') in a ``d'' column are aligned on the decimal point.
Use the ``D'' specifier for more complex layouts. }
\begin{ruledtabular}
\begin{tabular}{ccddd}
One&Two&
\multicolumn{1}{c}{\textrm{Three}}&
\multicolumn{1}{c}{\textrm{Four}}&
\multicolumn{1}{c}{\textrm{Five}}\\
\hline
one&two&\mbox{three}&\mbox{four}&\mbox{five}\\
He&2& 2.77234 & 45672. & 0.69 \\
C\footnote{Some tables require footnotes.}
  &C\footnote{Some tables need more than one footnote.}
  & 12537.64 & 37.66345 & 86.37 \\
\end{tabular}
\end{ruledtabular}
\end{table}