Search for heavy resonances decaying to a photon and a hadronically decaying Z = W = H boson in pp collisions at ffiffi s p = 13 TeV with the ATLAS detector

Many extensions of the Standard Model predict new resonances decaying to a Z,W, or Higgs boson and a photon. This paper presents a search for such resonances produced in pp collisions at ffiffi s p 1⁄4 13 TeV using a data set with an integrated luminosity of 36.1 fb−1 collected by the ATLAS detector at the LHC. The Z=W=H bosons are identified through their decays to hadrons. The data are found to be consistent with the Standard Model expectation in the entire investigated mass range. Upper limits are set on the production cross section times branching fraction for resonance decays to Z=W þ γ in the mass range from 1.0 to 6.8 TeV and for the first time into H þ γ in the mass range from 1.0 to 3.0 TeV.


Introduction
Many proposals for physics beyond the Standard Model (SM) include the prediction of new massive bosons. Examples are Technicolor [1] or Little Higgs [2], as well as extensions to the SM Higgs sector such as including an additional electroweak singlet scalar [3]. Decay modes of these new bosons include final states with a Z or a W boson and a photon. In addition, decays of heavy spin-1 bosons to the 125 GeV Higgs boson and a photon present an interesting search channel [4]. This paper describes a search for massive neutral and charged bosons decaying to a photon and a Z, W, or Higgs boson with subsequent hadronic decay of these bosons. The search uses 36.1 fb −1 of proton-proton (pp) collision data at a center-of-mass energy √ s = 13 TeV collected with the ATLAS detector in 2015 and 2016.
The selection of events collected for this search is based on the presence of high transverse energy photons. The identification of Z, W, and Higgs bosons exploits properties of the highly boosted bosons with merged dijet energy clusters reconstructed as a large-radius jet. The advantage of this final state is that a large fraction of events from the heavy resonance decay is detected since the branching fraction of Z and W bosons into hadrons is approximately 70%. The Higgs boson also decays mainly hadronically, dominated by the decay to a bb quark pair with a branching fraction of 58%. The measurements use the mass of the large-radius jet and other substructure information to identify the Z, W and Higgs bosons. In addition, for the Z and Higgs bosons, the sensitivity is increased by identify the long-lived decays of bottom hadrons within jets.
Previous searches for high-mass resonances decaying to a Z boson and a photon were conducted by Tevatron's D0 experiment [5] as well as at the Large Hadron Collider (LHC). The ATLAS experiment performed searches for heavy resonances decaying to a Z or a W boson and a photon, with subsequent leptonic decays of the Z and W bosons, using data collected at √ s = 7 and 8 TeV [6,7], and for heavy resonances decaying to a Z boson and a photon using 3.2 fb −1 of √ s = 13 TeV data [8]. In the latter search, both the leptonic and hadronic decays of the Z boson were used. The ATLAS experiment also used a 36.1 fb −1 dataset collected at √ s = 13 TeV to search for resonances decaying to a Z boson and a photon with the decays Z → e + e − , µ + µ − [9] and νν [10]. The CMS experiment performed searches for a heavy resonance decaying to a photon and a hadronically or leptonically decaying Z boson using datasets collected at √ s = 7, 8, and 13 TeV [11][12][13][14]. None of the searches revealed the presence of a new resonance. microstrip and transition radiation tracking detectors, covers the pseudorapidity1 range |η| < 2.5. The energies of photons and jets are measured primarily by the calorimeter system, which covers the pseudorapidity range |η| < 4.9. The electromagnetic calorimeter is a high-granularity liquid-argon (LAr) sampling calorimeter with lead absorber plates and is located just outside the solenoid. It spans the region |η| < 3.2 with barrel and endcap sections, segmented into three layers longitudinal in shower-depth in the region |η| < 2.5, with ∆η × ∆φ readout granularity in the second layer of 0.025 × 0.025. Beyond the cryostat of the electromagnetic calorimeter a steel/scintillator tile hadronic calorimeter covers the region |η| < 1.7. It is also segmented into three layers longitudinal in shower-depth with lateral readout granularity of 0.1 × 0.1 in ∆η × ∆φ. Two copper/LAr hadronic endcap calorimeters with similar granularity cover the region 1.5 < |η| < 3.2. In the forward region, electromagnetic and hadronic energy measurements are provided by copper/LAr and tungsten/LAr modules, respectively. A two-level trigger system [18] is used to select the events to be recorded. The first level of the trigger is implemented in hardware using a subset of the detector information to reduce the event rate to at most 100 kHz from the beam bunch crossing rate of 40 MHz. The final data selection is done with a software-based trigger that reduces the event rate to an average of 1 kHz.

Data and Monte Carlo samples
The data used in this search were collected in 2015 and 2016 from LHC pp collisions with a 13 TeV center-of-mass energy and 25 ns bunch spacing. In these collisions, the average number of inelastic pp collisions in each bunch crossing (referred to as pileup) was approximately 25. Events were recorded using a single-photon trigger that imposed a transverse energy threshold of 140 GeV and loose photon identification requirements based on cluster shower-shape variables [19]. The photon trigger is fully efficient for the events selected for this analysis. After requiring all ATLAS subdetectors to be operational, the resulting integrated luminosity is 36.1 fb −1 .
Simulated signal events are used to optimize the event identification and estimate the efficiency of the event reconstruction and selection. SM background processes were simulated to test the parameterization of the jet-photon invariant mass spectra, which is used in the data-driven estimation of the background. All simulated signal and background event samples were generated with Monte Carlo (MC) techniques as described below.
The production and decay of spin-0 and spin-2 Zγ resonances, spin-1 W γ resonances, and spin-1 Hγ resonances were modeled assuming a narrow-width approximation. The decay width of each resonance was set to 4 MeV, which is much smaller than the experimental resolution, and interference between these resonant processes and non-resonant SM production of the corresponding final states was neglected.
The Zγ scalar resonances produced via gluon-gluon fusion, gg → X → Zγ, were modeled by P -B [20, 21], with the CT10 parton distribution function (PDF) set [22], interfaced with P 8.186 [23] for the underlying event, parton showering and hadronization with the CTEQ6L1 PDF set [24] and a set of tuned underlying-event parameters called the AZNLO tune [25]. This model is the same as that used for SM H → Zγ production with the resonance mass varied but the width held fixed.
The spin-2 gg → X → Zγ and qq → X → Zγ, and the spin-1 qq → X → W γ and qq → X → Hγ resonant processes were modeled via effective couplings implemented in M G 5_aMC@NLO v2.3.3 [26], interfaced to the P 8.186 parton shower model with the A14 parameter tune [27] and the NNPDF2.3LO PDF set [28]. The Zγ models produce transversely polarized Z bosons and the W γ models produce longitudinally polarized W bosons. In these samples, the W and Z bosons are forced to decay hadronically, and the Higgs boson is forced to decay to bb.
The dominant SM backgrounds are prompt photons produced in association with jets, hadronically decaying W or Z bosons produced in association with a photon, and tt + γ events. The samples of events containing a photon with associated jets were simulated using the S 2.  [37]. In all simulated signal and background samples, the effects of overlapping inelastic pp collisions were included by overlaying multiple events simulated with P 8.186 using the A3 set of tuned parameters [38] and the MSTW2008LO PDF set [39]. The simulated events were weighted so that the distribution of the number of pileup interactions in the simulation matched the one in the data. The simulated events were then passed through the same event reconstruction algorithms used for the data, including corrections for known differences between data and simulation in the efficiencies of photon reconstruction and selection, in the photon and jet energy scale and resolution, and in the tagging efficiency of heavy-flavor jets.

Reconstruction of photons and jets
The photon reconstruction is described in Ref. [40]. A discriminant, based on lateral and longitudinal shower profiles, is constructed to distinguish prompt photons from hadrons as well as photons from decays of mesons inside jets. In this analysis, two levels of selections are applied: the Loose selection criterion defined in Ref.
[40] at the trigger level, and the Tight selection criterion for the final analysis selection. In the final selection, photon candidates are required to have transverse energy above 250 GeV and pseudorapidity |η| < 1.37. These criteria are applied to reduce the contribution from the SM production of prompt photons with associated jets. The efficiency of a photon within that region to pass the Tight selection criterion is above 90%. The transverse energy E T,iso , deposited within a cone of size ∆R ≡ (∆η) 2 + (∆φ) 2 = 0.4 around the photon cluster, is corrected for the photon energy deposited outside of the cluster, underlying event, and multiple pp interactions; it is required to satisfy E T,iso < 2.45 GeV + 0.022 E γ T , where E γ T is the transverse energy of the photon cluster. The requirement is applied to reduce contamination from hadrons misidentified as photons and is 98% efficient for prompt photons passing the Tight criterion. Corrections mitigating the differences between the calorimeter response to photons in data and simulation are derived using collision data. The photon energy calibration is derived using samples of simulated events followed by data-driven corrections [41].
Large-radius (large-R) jets are reconstructed from topological energy clusters (topocluster) [42] in the calorimeter, using the anti-k t algorithm [43] with a radius parameter R = 1.0. To reduce contributions to the jet transverse momentum (p T ) and mass arising from pileup, a trimming procedure [44] is applied where subjets with R = 0.2 and carrying less than 5% of the original p T of the jet are removed. Jets are calibrated to the level of stable final-state particles using simulation [45]. Differences between data and simulation in jet energy scale and resolution are corrected using in situ methods [46]. The jet energy resolution for jets with p T of 1 TeV is approximately 5%. Jet candidates are required to have p T above 200 GeV and |η| < 2.0 to ensure tracking detector coverage within the jet cone. The jet candidates are required to be separated from any photon candidate by ∆R > 1.0.
The jet mass is computed as a weighted combination of calorimetric mass and track-assisted mass [47]. The calorimetric mass is computed from the massless topocluster four-momenta. The track-assisted mass incorporates information from the calorimeter and the four-momenta of tracks, which are matched to the jets using ghost association [48] and are matched to the primary vertex, which is defined as the vertex with largest sum of squared momenta of associated tracks in the event [49]. The relative weighting of the calorimetric and track-assisted mass in the final mass calculation is based on the expected resolutions of the two mass variables. The mass resolution in the peak of the jet mass distribution for jets originating from Z, W or Higgs bosons ranges from 7% to 15% for jets with transverse momenta of 500 GeV to 2500 GeV. Reconstructed jet mass distributions for various signal hypotheses are shown in Figure 1. Peaks centered at the Z, W, and Higgs boson masses are clearly visible. In the mass distributions of jets arising from the Z and Higgs boson decays, a feature at low jet invariant mass is also prominent for decays of resonances with masses of 2 TeV or lower. This is a result of energy flow outside of the jet cone. For the W γ spin-1 resonance decays the effect is mitigated due to the longitudinal polarization of the W boson, which enhances the collimation of the decay products. The lower mass side of the cores of the Z boson and Higgs boson mass peaks are enhanced from loss of neutrino energy in the jet when at least one of the b-hadrons decays semileptonically.
The jet mass and the D (β=1) 2 jet substructure discriminant are employed to distinguish between jets originating from hadronically decaying Z or W bosons and jets originating from quarks and gluons. The variable D (β=1) 2 is defined as the ratio of two-point and three-point energy correlation functions [50,51], which are based on the energies and pairwise angular distances of particles within a jet. The performance of this discriminant has been studied in MC simulations and data [52,53]. Upper and lower bounds on the jet mass and upper bound on D (β=1) 2 are tuned to achieve 50% efficiency for jets with p T in the range of 300 GeV to 2500 GeV from decays of a Z or W boson. The fraction of jets originating from a quark or a gluon passing this selection varies between approximately 2.2% and 1.3% in this p T range. The D (β=1) 2 discriminant is not used for the H → bb selection.
Track-jets are reconstructed using the anti-k t algorithm with radius parameter R = 0.2 and tracks matched to the large-R jets using ghost association.Track-jets that contain b-hadrons are identified using the MV2c10 tagging algorithm, which exploits the lifetime of b-hadrons and the kinematic properties of their charged decay products [54,55]. The efficiency of the algorithm is 70% when applied to b-jets in simulated tt events. The fraction of track-jets originating from a light quark or a gluon, tagged as originating from a b-hadron, is approximately 0.8% in simulated tt events.

Event selection and categorization
Events considered for analysis are triggered by the presence of a photon candidate with p T greater than 140 GeV. A small fraction of events where adverse instrumental effects were identified are removed. The baseline selection identifies events with one photon candidate and one large-R jet candidate. If more than one photon or jet candidate are found, only the highest-p T objects are considered further. The efficiency of the selection for a resonance with a mass of 3 TeV varies from approximately 60% to 80% depending on the signal hypothesis, as shown in Figure 2. The decrease of the baseline selection efficiency for resonances with lower mass is due to the kinematic thresholds of 200 GeV and 250 GeV required for the jet and photon p T respectively. Differences in the baseline efficiencies between the resonance types are due to different angular distributions of the produced photon-jet system (depending on the spin hypothesis and production mode) and therefore different probabilities to pass the photon and jet p T and |η| requirements.
Following the baseline selection, events are classified into four or fewer subsamples to improve the expected signal sensitivity. The categorization is made in order of decreasing background rejection to achieve high sensitivity throughout the resonance mass (m X ) range. For resonance masses below 3 TeV it is desirable to maximally suppress the SM background, while for very high m X values, due to the steeply falling jet-photon invariant mass (m Jγ ) distribution of background processes, a loose selection is appropriate. For the Zγ search the categories, defined in the next paragraph, are: BTAG, D2, VMASS, and ELSE. In the W γ search only the D2, VMASS, and ELSE categories are used. The Hγ search employs only the BTAG category. The process of event categorization proceeds sequentially starting from the events found in the baseline selection as shown in the diagrams in Figure 3.
The first subsample in the Zγ and Hγ searches, the BTAG category, captures events where the two leading track-jets associated with the selected large-R jet candidate are b-tagged, exploiting the decays of Z and Simulation ATLAS = 13 TeV s Figure 2: Efficiencies for signal events to pass the baseline selection as a function of the resonance mass m X for different signal models, production modes, and spin hypotheses. The uncertainty bars shown are statistical only. The lines represent polynomial fits to the simulated data points. The line corresponding to the Hγ resonance baseline selection efficiency is discontinued above m X of 3 TeV since the search for Hγ does not cover this resonance mass region.
H bosons to a bb quark pair together with strong background suppression. A window requirement on the jet mass is also applied. In the Zγ search the mass interval grows from 80−106 GeV for jets with p T of 500 GeV to 70−110 GeV for jets with p T of 2.5 TeV. These mass intervals are varied such that an approximately constant signal efficiency across the whole jet p T spectrum is maintained, accounting for the jet mass resolution increase as a function of jet p T . In the Hγ search the mass interval is 93−134 GeV independently of the jet p T . In the Zγ and Hγ searches the jet must also have fewer than 30 associated tracks (n trk ) originating from the primary vertex for the event to be accepted in this category. This requirement is made to reject gluon-initiated jets that mimic a two-subjet structure due to gluon splitting [56]. Relative to the baseline selection, the efficiency of selecting in the BTAG category an event originating from a Zγ resonance with a mass of 1 TeV, in the hadronic Z boson decay mode, is 3% to 4% depending on the spin hypothesis and production mode, while only 0.02% of background events enter this category for the same mass value. For the Hγ resonances with a mass of 1 TeV, in the H → bb decay mode, the selection efficiency is 25%. As shown in Figure 4, the BTAG category becomes ineffective for capturing signal events originating from resonances with masses higher than approximately 3 TeV due to b-tagging inefficiency for highly boosted jets. In the Hγ search the categorization process is stopped at this stage, with the remaining events in the baseline selection being rejected from further analysis.
Events not entering the BTAG category in the Zγ search and events from the baseline set in the W γ search are placed in the D2 category if the selected jet satisfies the combined jet mass and D (β=1) 2 discriminant requirements. In the W γ search the n trk < 30 requirement must also be satisfied. The mass window requirement in the Zγ search is the same as in the BTAG category. In the W γ search the mass interval grows from 71−95 GeV for jets with p T of 500 GeV to 60−100 GeV for jets with p T of 2.5 TeV. Approximately 20% to 25% (depending on the spin hypothesis and the production mode) of 1 TeV Zγ resonance decays passing the baseline selection enter this category with the fraction increasing to 22%−28% for 4 TeV resonances. For the W γ decays the fraction is approximately 40% for m X below 4 TeV. The difference in the D2 categorization efficiency between the Zγ and W γ processes is explained by differences between the angular distributions of the Z and W boson decay products due to the longitudinal polarization of the W boson and transverse polarization of the Z boson. For both the Zγ and W γ resonances the D2 categorization efficiency decreases for m X above 4 TeV, as shown in Figure 4, due to the falling discriminating power of the D (β=1) 2 variable. Only approximately 1% (0.5%) of background events passing the baseline selection for resonance masses of 1 TeV (4 TeV) enter this category.
Events that fail to enter either of the first two categories and only pass the jet mass selection enter the VMASS category. The efficiency to enter this category, relative to the baseline selection, is 24% to 27% for Zγ events from a resonance at 1 TeV, growing to 35−36% for a resonance at 4 TeV, and approximately 30% for W γ events. Approximately 9% of the background events enter this category after passing the baseline selection.
Finally, the remaining events passing the baseline selection are assigned to the ELSE category. The Zγ resonance events enter this category with an efficiency of 40% to 50%, and W γ events enter with approximately 30% efficiency for resonances below 4 TeV, growing to 50% for a resonance at 7 TeV.

Signal and background models
The final discrimination between signal and background events in the selected samples is achieved with a fit of a signal+background model to the m Jγ distribution of the selected data events. The fit relies on the parameterization of signal and background m Jγ distributions with a functional form.

Signal model
The shape of the m Jγ distribution is modeled by the sum of a Crystal Ball function [57] representing the core populated by well-reconstructed events, and a Gaussian function modeling the tails populated by poorly reconstructed events as described by Eq. (1): where C is the Crystal Ball function defined by Eq. (2): and G is a Gaussian function with mean µ and standard deviation σ G . The normalization factor N C depends on the shape parameters of the Crystal Ball function. The relative strength of the two distributions, f C , is  a parameter of the signal description model. Since the modeling of the tails of the mass peak is addressed by the Gaussian function, the parameter n C , controlling the shape of the tail of the Crystal Ball function, is fixed to 1. Similarly, the parameter α C , designating the onset of the power law tail of the Crystal Ball function is constrained to be between 0 and 4. Additional free parameters σ C and σ G describe the widths of the Crystal Ball and the Gaussian distributions corresponding to the width of the core of the m Jγ distribution and the width of the tails.
The model described above is fitted to the m Jγ distribution of simulated events of each signal type considered, for m X ranging from 850 GeV to 7 TeV, for each category. To obtain a model varying continuously as a function of the resonance mass, the parameters are interpolated using polynomials of up to the third degree. For all signal hypotheses, σ C is about 20 GeV for a 1 TeV resonance, and grows linearly by 15 GeV per 1 TeV increase of the resonance mass. Simulation ATLAS Simulation ATLAS Simulation ATLAS Figure 4: Efficiencies relative to the baseline selection for signal events to pass the category selection for (a) the Zγ (in the spin 0 production hypothesis), (b) the W γ, and (c) the Hγ resonances as a function of the resonance mass m X . The lines represent polynomial fits to the simulated data points. In the Zγ search the categorization efficiencies have only a small dependence on the production and spin hypotheses. Statistical uncertainties are smaller than the marker size.

Background model
Several SM processes contribute to the predicted event yield with different proportions. In the BTAG category, the dominant SM process is photon production in association with a b-flavored hadron, whereas in the other categories photon production in association with a light or c-flavored hadron dominates. The production of a photon in association with a Z or W boson, or the pair production of top quarks also contribute to the total background. For events with m Jγ > 1 TeV the contribution in the BTAG category from γ + Z production is approximately 32% (13%) for the Zγ (Hγ) search, while in the D2 category the contribution from γ + W/Z production is approximately 15% for both the Zγ and W γ searches. For other categories the contributions from SM γ + W/Z production are below 5%.
Samples of simulated events arising from the processes described above are used to develop the functional modeling of the background and to test the applicability of the functional form, number of parameters, and the range of the fit. Multijet production, where one jet is reconstructed as a photon, also contributes to the event samples. The contribution of this type of events is estimated using a data-driven method [8] and shown to be about 10% of the events passing the baseline selection and to not affect the m Jγ distribution. Multijet production is therefore accounted for in the background model through the γ+jet production.
The family of functions from Ref.
[58], as described by Eq. (3), with up to five parameters, is used for the overall background model: where x is m Jγ divided by the collision energy, and p ≡ (p 1 , p 2 , p 3 , p 4 , p 5 ) is the vector of shape parameters.
In the VMASS and D2 categories, the fit spans the m Jγ range of 800 GeV to 7 TeV; in the BTAG category, it spans the range from 800 GeV to 3.2 TeV; and in the ELSE category, it spans the range from 2.5 TeV to 7 TeV.

Systematic uncertainties
Uncertainties from systematic effects are due to the background estimation as well as the detector modeling, which affect the shape and normalization of the signal m Jγ distributions. These effects are estimated as relative uncertainties for the signal efficiency and the position and width of the signal peak for resonance masses ranging from 1 TeV to 7 TeV in 500 GeV steps. The impact of these effects on the m Jγ distribution is evaluated using simulation. Third-order polynomial interpolations are used to obtain relative variations of the signal efficiency and shape parameters due to each systematic uncertainty, for an arbitrary value of the resonance mass.
Uncertainty in background estimate: A systematic uncertainty associated with the background description arises from a potential bias in the estimated number of signal events due to the functional form chosen for the background paramaterization described by Eq. (3). This effect, referred to as the spurious signal (N spur ), is studied using simulated background events (including γ+jet and W/Z+ jet). The bias is estimated by fitting the signal-plus-background model to simulated background m Jγ distributions in each event category, for each signal hypothesis, with the sample's statistical uncertainties as expected in the data. The absolute number of fitted signal events at a given m X hypothesis defines the number of spurious signal events N spur (m X ). The impact of uncertainties in the background composition has been studied by varying the fraction of the W/Z+ jet and γ+b/c jet backgrounds by 50%, and found to be negligible in the spurious signal estimate. The impact of the spurious signal uncertainty on the exclusion limits is discussed in Section 8.
Luminosity: An integrated luminosity uncertainty of 2.1% is derived, following a methodology similar to that detailed in Ref. [60], from a calibration using beam separation scans performed in August 2015 and May 2016.

Jet energy scale and resolution:
The uncertainties in the jet energy scale and resolution are estimated using γ+jet and dijet events in the data [46].
The impact of the systematic uncertainty in the jet energy scale is a shift of the peak position of the signal m Jγ distribution by 1−3%. The signal mass resolution varies by 5% in the low-mass region (m X < 2.5 TeV) and by 15% in the high-mass region due to the systematic uncertainty in the jet energy resolution. The impact on the signal efficiency from the jet energy uncertainty is 2−6%.
Photon energy scale and resolution: The uncertainties in the photon energy scale and resolution are estimated using electron data samples with Z → ee events and high-purity photon samples with radiative Z → eeγ events [41]. The impact of the systematic uncertainty in the photon energy scale is a shift of 0.5% in the peak position of the signal m Jγ distribution. The signal mass resolution varies by 1% in the low-mass region and by 3% in the high-mass region due to the systematic uncertainty in the photon energy resolution.

Photon identification, isolation, and trigger efficiency:
The uncertainties in the reconstruction, identification, isolation, and trigger efficiency for photons are determined from data samples of Z → γ, Z → ee, and inclusive photon events, using the methods described in Ref. [19]. The impact on the signal efficiency from the photon identification, isolation, and trigger systematic uncertainties is found to be less than 1.5%, 0.5%, and 0.1%, respectively.
Heavy-flavor jet identification: Uncertainties in the b-tagging efficiency for track-jets are derived from the uncertainties measured in dedicated heavy-flavor-enriched data samples, following the methodology described in Ref. [54]. The uncertainties are measured as a function of b-jet p T and range between 2% and 8% for track-jets with p T < 250 GeV. For track-jets with p T > 250 GeV the uncertainty in the tagging efficiencies is extrapolated using simulation [54], and is approximately 9% for track-jets with p T > 400 GeV. The impact of these uncertainties on the signal efficiency is 10−20%.

Number of primary-vertex-originated tracks associated with the jet:
The requirement on the number of tracks originating from the primary vertex and associated with the jet induces a 6% systematic uncertainty in the signal efficiency, as estimated from the comparison of data control samples and samples of simulated events [56]. scale and resolution uncertainties on the signal efficiency and the core width σ C of the signal peak is found to be less than 1%.

Pileup modeling:
The pileup weighting of the simulated signal events, described in Section 3, is varied to cover the uncertainty in the ratio of the predicted and measured inelastic pp cross sections [62]. The pileup uncertainty affects the signal efficiency by 1−2%.

PDF choice:
The uncertainty due to the PDF modeling is evaluated by comparing the signal acceptances for alternative PDF sets to that for the nominal set. The total uncertainty in the acceptance is derived as the standard deviation of the eigen-variations according to the method described in Ref. [63]. The uncertainties in acceptances for signal processes produced via qq annihilation vary from 5% to 2% with increasing resonance mass, while for signal processes produced via gluon-gluon fusion the uncertainties vary from 12% to 2%.

Parton shower:
The uncertainty due to the parton shower modeling is evaluated by comparing the signal acceptances for alternative parton shower models to the acceptance for the nominal model. The alternative parton shower models are defined as the eigen-variations of the P A14 tune [27]. The total uncertainty is derived as the sum in quadrature of all the eigen-variation effects and affects the acceptance by about 2%.
The effects of systematic shifts on signal normalization, position of the signal peak, and the core width σ C of the signal peak are summarized in Table 1 for the Zγ, W γ, and Hγ searches. The effects on the signal acceptances from theoretical uncertainties are also given.  Table 1: Effect of systematic uncertainties from various sources on signal normalization and efficiency, position of the signal peak, and the core width σ C of the signal peak. The last two rows show the theoretical uncertainty effects on the signal acceptances.

Statistical procedure
The data are scrutinized with statistical methods to quantify the presence of a hypothetical resonance and to set a limit on its production. In both cases an unbinned extended maximum likelihood estimator is used to model the data. The parameter of interest is σB had , the signal cross section times the branching fraction of the resonance decaying to (Z/W/H)γ with subsequent hadronic (bb) decays of Z/W (Higgs) bosons. The impact of systematic uncertainties on the signal is modeled with a vector of nuisance parameters, θ, where each component, θ k , is constrained with corresponding Gaussian probability density functions G k (θ k ). The likelihood model, L, for the sample of data events is described by Eq. (4): where j represents the event category, n j is the observed number of events in that category, is the event index, and N j (σB had , θ) is the expected total event yield in category j. The total probability density function f tot, j in category j depends on the photon-jet invariant mass m Jγ, and is a function of the parameter of interest σB had , the parameters of the background modeling function p, the nuisance parameters θ k , as well as the mass m X of the hypothetical resonance. The functional form of f tot, j is given in Eq. (5): where S j and B j are the signal and background probability density functions in category j, described in Section 5. The parameter θ spur, j is an element of the θ vector corresponding to the spurious signal nuisance parameter. The expected yield of signal events N sig, j is given by the product of σB had , integrated luminosity, acceptance and efficiency for a given category j. The expected number of background events N bkg, j is a parameter of the fit. The total expected event yield in a given category, N j , is the sum of the expected signal, background, and spurious signal event yields.
The p-values are computed to examine the compatibility of the data and the background-only hypothesis. First, the local p-value is calculated for the particular value of m X under consideration. The local pvalue is defined as the probability of the background to produce a signal-like excess whose estimated σB had is larger than that found in the fit to the data in all categories simultaneously. This procedure utilizes the ratio of the likelihood value where the most likely value of the parameter of interest σB had is found, to the likelihood value where no signal is allowed (σB had = 0) [64]. The global p-value is defined as the probability of finding, at any value of m X , a signal-like fluctuation more significant than the most significant excess found in the data, in all categories combined. It is calculated approximately by discounting the local p-value by the effective number of search trials possible within the m Jγ range examined.
The modified frequentist (CL s ) method [65,66] is used to set upper limits on the signal σB had at 95% confidence level (CL). To obtain limits on σB, which is the signal cross section times the branching fraction of resonance decays to (Z/W/H)γ with all allowed decays of Z, W, and Higgs bosons, the resulting σB had is divided by the measured hadronic branching fraction of the Z (69.91% [67]) or W (67.41% [67]) bosons or the theoretically calculated branching fraction of the Higgs boson to bb (58.24% [68]).
Closed-form asymptotic formulae [64] are used to calculate the limits. The p-value calculation and σB limit calculations are performed for m X in range of 1 TeV to 6.8 TeV in steps of 20 GeV. The step size is chosen to be much smaller than the experimental m Jγ resolution and the interval spans the interpolation range of the shape parameters. The VMASS and D2 categories are used in the entire range, while the BTAG category is used only for 1 TeV < m X < 3 TeV, and the ELSE category is used only for 3 TeV < m X < 6.8 TeV. The fit interval described in Section 5.2 is selected such that a peak in the m Jγ distribution, resulting from resonance with m X considered, is fully contained.
Due to the small number of events for large m X values, the results are checked with ensemble tests. The limits derived with asymptotic formulae agree well with those obtained with ensemble tests for small m X values, and are underestimated by as much as 30% for large m X values.

Results
The data event yields in the baseline selection as well as in all search channels and in the categories within each channel are shown in Table 2.
The observed m Jγ distributions as well as the background distributions obtained from a global backgroundonly fit with various hypothetical signal curves overlaid are shown in Figures. 5, 6, and 7 for the Zγ, W γ, and Hγ searches, respectively. The event with the highest m Jγ , at a value of 6.3 TeV, has a jet of mass 63 GeV and it is found in the ELSE (VMASS) category of the Zγ (W γ) search.
The smallest local p-value, corresponding to a significance of 2.7 σ, is found in the W γ search at m Jγ = 2.5 TeV utilizing data from all categories simultaneously. This local p-value corresponds to a Hγ search 60,237 59 · · · · · · · · · global significance of less than 1σ. No significant excess is observed in any of the categories and analysis channels. Limits are placed on specific models.
The σB limits on Zγ production, evaluated for resonance masses between 1.0 and 6.8 TeV, are shown in Figure 8. Spin and production hypotheses comprise a spin-0 resonance, produced by gluon-gluon fusion and a spin-2 resonance, produced by gluon-gluon fusion and qq annihilation. The σB limit on W γ resonances evaluated for m X between 1 and 6.8 TeV are shown in Figure 9. This is the first evaluation of such a limit utilizing hadronic W boson decays. The Zγ and W γ limits decrease from approximately 10 fb for a resonance mass of 1 TeV to 0.1 fb for a resonance mass of 6.8 TeV. The Hγ search presented here is the first search for a heavy resonance with this decay mode. The σB limits on the resonances decaying to a Hγ final state are shown in Figure 10. The limit is evaluated for resonance masses between 1 and 3 TeV and varies between 10 fb and 4 fb depending on m X . The limit weakens for m X > 2 TeV due to a decrease in signal efficiency as shown in Figure 4, stemming from the decrease in the b-tagging efficiency for high-momentum jets.
The sensitivity of the resonance search and the strength of the resonance production cross section limit are primarily determined by the available data sample size. Among all the systematic uncertainties, the spurious-signal uncertainty on the background estimation has the largest impact on the limit, in particular in the low-mass region, weakening the limit by up to 20% (1%) at m X = 1 TeV (6.8 TeV). Another important systematic uncertainty is that in the heavy-flavor jet identification efficiency. It weakens the limit by up to 13% (20%) at m X = 1 TeV (3 TeV) in the Hγ analysis, while it has little impact on the limits in the Zγ analysis since the BTAG category is just one of the four categories in this analysis. Events / 40 GeV   Events / 40 GeV Events / 40 GeV

Conclusion
Results are presented from a search for heavy resonances decaying to Zγ, W γ or Hγ final states using 36.1 fb −1 of √ s = 13 TeV pp collision data collected by the ATLAS experiment during the 2015 and 2016 at the LHC run periods. The search sensitivity is enhanced for high-mass resonances by selecting hadronic decays of the Z and W bosons, and the bb decay of the Higgs boson, identified as large-radius jets. Distributions of the invariant mass of photon-jet pairs are used to search for resonances in the mass range from 1.0 to 6.8 TeV for decays to Zγ and W γ, and between 1.0 and 3.0 TeV for decays to Hγ. No evidence for new resonances is found and limits are set based on assumptions on the spin and production model of the resonance.
The 95% confidence level upper limits on the resonance production cross section times decay branching fraction to Zγ and W γ final states vary from about 10 fb to 0.1 fb for masses from 1.0 to 6.8 TeV, with individual studies carried out for resonances with spin 0, 1, and 2 produced via gluon-gluon fusion and qq annihilation. For the spin-1 X → Hγ search the limits vary from about 10 fb to 4 fb for resonance masses between 1.0 and 3.0 TeV. These results set the first limits on the production of Hγ resonances, and this search covers a wider mass range and has a broader scope than previous searches for heavy resonances decaying to Zγ and W γ final states.
(Taiwan), RAL (UK) and BNL (USA), the Tier-2 facilities worldwide and large non-WLCG resource providers. Major contributors of computing resources are listed in Ref. [70].