A GPU based multidimensional amplitude analysis to search for tetraquark candidates
Nairit Sur, Leonardo Cristella, Adriano Di Florio, Vincenzo Mastrapasqua
SSur et al.
RESEARCH
A GPU based multidimensional amplitudeanalysis to search for tetraquark candidates
Nairit Sur † , Leonardo Cristella † , Adriano Di Florio † and Vincenzo Mastrapasqua † * Correspondence:[email protected] Tata Institute of FundamentalResearch, Mumbai, IndiaFull list of author information isavailable at the end of the article † Equal contributor
Abstract
The demand for computational resources is steadily increasing in experimentalhigh energy physics as the current collider experiments continue to accumulatehuge amounts of data while physicists indulge in more complex and ambitiousanalysis strategies. This is especially true in the fields of hadron spectroscopy andflavour physics where the analyses often depend on complex multidimensionalunbinned maximum-likelihood fits, with several dozens of free parameters, withthe aim to study the quark structure of hadrons.Graphics processing units (GPUs) represent one of the most sophisticated andversatile parallel computing architectures that are becoming popular toolkits forhigh energy physicists to meet their computational demands. GooFit is anupcoming open-source tool interfacing ROOT/RooFit to the CUDA platform onNVIDIA GPUs that acts as a bridge between the MINUIT minimization algorithmand a parallel processor, allowing probability density functions to be estimated onmultiple cores simultaneously.In this article, a full-fledged amplitude analysis framework developed usingGooFit is tested for its speed and reliability. The four-dimensional fitterframework, one of the firsts of its kind to be built on GooFit, is geared towardsthe search for exotic tetraquark states in the B → J/ψKπ decays that can alsobe seamlessly adapted for other similar analyses. The GooFit fitter running onGPUs shows a remarkable speed-up in the computing performance whencompared to a ROOT/RooFit implementation of the same, running on multicoreCPU clusters. Furthermore, it shows sensitivity to components with smallcontributions to the overall fit. It has the potential to be a powerful tool forsensitive and computationally intensive physics analyses.
Keywords:
High Energy Physics; Flavour Physics; Exotic Hadron Spectroscopy;GPU; CUDA; GooFit; Unbinned Maximum Likelihood
Introduction
Quantum chromodynamics (QCD) is the quantum field theoretical description ofstrong interaction between the quarks and gluons that make up composite hadronsobserved in nature. Due to the complex nature of this non-abelian gauge theorywith peculiar features like “colour confinement” and “asymptotic freedom” [1–4],it is very hard to study the nature of this interaction analytically, especially at lowenergy regimes. In the last 15 years, experimental evidences for a large variety ofnew multiquark bound states, allowed by QCD, have been found that do not fit theexpectations for the conventional quark model (qq (cid:48) q (cid:48)(cid:48) baryons or q¯q (cid:48) mesons). Theexact nature of many of these states still remains a puzzle; even though some ofthem are confirmed by multiple experiments, all the quantum numbers of each state a r X i v : . [ h e p - e x ] J u l ur et al. Page 2 of 16 are still not experimentally determined. Spectroscopic studies of such heavy-flavorstates can provide a deeper understanding of the underlying dynamics of quarksand gluons at the hadron mass scales as well as a valuable insight into various QCDinspired phenomenological models [5, 6].The manifestly exotic charged charmonium-like Z states, which are strong candi-dates for tetraquark states with a possible quark content of | c ¯ cd ¯ u (cid:105) , can be studiedin ongoing collider experiments, namely ATLAS, Belle II, BESIII, CMS, and LHCb.To ascertain the presence of such states with a high degree of significance in three-body decays B → ψ ( nS ) K + π − , complex multidimensional unbinned maximum-likelihood (UML) fits on tens of thousands of data points, with several dozens of freeparameters, must be performed, thus requiring a considerable amount of computa-tional resources. The traditional high-energy physics (HEP) analysis tools such asROOT [7] and its RooFit package [8] may require very long processing times (evendays, depending on the complexity of the model under investigation) even whenthey are run on big clusters comprising multiple CPUs.In this article, we explore the scope of an advanced GPU-accelerated comput-ing setup to reduce the processing times of such complex multidimensional fitsfrequently occurring in the HEP field. An amplitude analysis of the three body decay B → J/ψKπ
The rare exotic Z states can appear as J/ψπ resonances in the quasi 2-body decay B → Z − K + → J/ψπ − K + , where the J/ψ decays into a µ + µ − pair (inclusionof the charge conjugate mode ¯ B → Z + K − → J/ψπ + K − is always implied).However, the decay process is dominated by the intermediate K ∗ ( → Kπ ) resonancesin the quasi 2-body decay B → J/ψK ∗ [9]. These ten kinematically allowed kaonicresonances can interfere with one another as well as with the Z states.Three-body decays with intermediate resonant states, such as P → D + D res , D res → D + D , are generally analysed using a technique pioneered by Dalitz [10].Here, P is the parent particle, D is one of its daughters, D res is the other daughterwhich, being an intermediate resonance, decays into D and D . A two-dimensionalscatter plot of m D D vs. m D D (invariant mass squared of any two daughters),known as Dalitz plot, shows a nonuniform distribution due to the interfering in-termediate resonances, thus to the decay dynamics. If at least one of the threedaughters in the decay is a vector state instead of a pseudoscalar, the traditionalDalitz plot approach becomes insufficient as the angular variables are implicitlyintegrated over, leading to a loss of information about angular correlations amongthe decay products. The K ∗ -only model The kinematics of the process B → J/ψKπ , J/ψ → µ + µ − can be completelydescribed by a four-dimensional variable spaceΦ ≡ (cid:0) m Kπ , m J/ψπ , θ
J/ψ , ϕ (cid:1) . (1)The two angles, θ J/ψ and ϕ are illustrated in Fig. 1.The total decay amplitude of B → J/ψKπ is represented by a coherent sum ofthe Breit-Wigner (BW) contributions associated with all the kinematically allowed ur et al. Page 3 of 16
Figure 1
A sketch illustrating the definition of two independent angular variables, θ J/ψ and ϕ , forthe amplitude analysis of B → J/ψKπ decays. intermediate resonant states. The angle-independent part of the decay amplitudefor each resonance R is given by [11]: A R (cid:0) m R (cid:1) = F ( L B ) B (cid:16) p B M B (cid:17) L B F ( L R ) R (cid:16) p R m R (cid:17) L R M R − m R − iM R Γ( m R ) , (2)where the mass-dependent width of the resonance R is:Γ( m R ) = Γ (cid:18) p R p R (cid:19) L R +1 (cid:18) M R m R (cid:19) F R (3)and • m R is the running invariant mass of the two daughters of R (e.g., m R = m Kπ for a K ∗ ); • M B is the B meson mass; • M R is the nominal mass of R ; • L B ( L R ) is the orbital angular momentum in the B ( R ) decay; • p B is the B daughter momentum (i.e R momentum) in the B rest frame; • F ( L B ) B and F ( L R ) R are the Blatt-Weisskopf form factors [12] for B and R decay,respectively, with the superscript denoting the orbital angular momentum ofthe (sub-)decay; • Γ is the nominal width of R ; • p R and p R are the momenta of R daughters in the former’s rest frame, cal-culated from the running and pole mass of R , respectively.The angle-dependent part of the amplitude is obtained using the helicity formal-ism [13]. For each K ∗ resonance, it is given by A K ∗ λξ (Φ) = H K ∗ λ A K ∗ (cid:0) m Kπ (cid:1) d J ( K ∗ ) λ ( θ K ∗ ) e iλϕ d λξ ( θ J/ψ ) (4)where, A K ∗ (cid:0) m Kπ (cid:1) , defined in Eq. (2), is explicitly written for R ≡ K ∗ and ur et al. Page 4 of 16 • J ( K ∗ ) is the spin of the considered K ∗ resonance; • λ is the helicity of the J/ψ (the quantisation axis being parallel to the K ∗ momentum in the J/ψ rest frame). In general, λ can take the values −
1, 0and 1. For K ∗ s with zero spin, only λ = 0 is allowed; • ξ is the helicity of the µ + µ − system; • H K ∗ λ is the complex helicity amplitude for the decay via the intermediate K ∗ ; • d J ( K ∗ ) λ ( θ K ∗ ) and d λξ ( θ J/ψ ) are the Wigner small-d functions; • θ K ∗ is the K ∗ helicity angle, i.e. the angle between K momentum the in K ∗ rest frame and the K ∗ momentum in the B rest frame (Fig. 1); • θ J/ψ is the
J/ψ helicity angle, i.e. the angle between µ + momentum in the J/ψ rest frame and
J/ψ momentum in the B rest frame; and • ϕ is the angle between the J/ψ → µ + µ − and K ∗ → Kπ decay planes.The signal density function, to be used in the UML fit, is obtained after appro-priately summing over the helicity states and is given by S (Φ) = (cid:88) ξ =1 , − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) K ∗ (cid:88) λ = − , , A K ∗ λξ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (5)The sum over K ∗ includes all kinematically allowed resonance states up to m Kπ =2 .
183 GeV, namely K ∗ (800), K ∗ (892), K ∗ (1410), K ∗ (1430), K ∗ (1430), K ∗ (1680), K ∗ (1780), K ∗ (1950), K ∗ (1980), and K ∗ (2045). As the expression in Eq. (5) issensitive only to the relative phases and amplitudes, we have the freedom to fix oneoverall phase and amplitude in the fit. The helicity amplitude of the K ∗ (892), thedominant resonance, is chosen to be fixed, for λ = 0: (cid:12)(cid:12)(cid:12) H K ∗ (892)0 (cid:12)(cid:12)(cid:12) = 1 , arg (cid:16) H K ∗ (892)0 (cid:17) = 0 (6)The masses and widths of all the resonances are fixed to their world-average val-ues [14]. The LASS parametrization
Generally, P- and D-waves states are considered to be well described by narrowresonance approximations. For the Kπ system, the low mass S-wave K ∗ (800) ap-pears as a broad peak calling for a more careful treatment. The LASS experiment atSLAC used an effective range expansion to model the low-energy behaviour of such Kπ S-wave [15]. We use a similar parametrization where the angle-independent partof the amplitude is a nonresonant contribution interfering with the scalar K ∗ (1430)BW amplitude: A LASS = m Kπ q Kπ sin θ B e iθ B + 2e iθ B (cid:16) m K ∗ (1430) /q K ∗ (1430) (cid:17) Γ K ∗ (1430) M K ∗ (1430) − m Kπ − iM K ∗ (1430) Γ( m Kπ ) , (7)with cot θ B = 1 a q Kπ + 12 b q Kπ and, a = 1 .
95 GeV − , b = 1 .
76 GeV − , (8)where ur et al. Page 5 of 16 • m Kπ is the running mass of the Kπ system; • q Kπ is the momentum of one of the K ∗ daughters in the K ∗ rest frame; • Γ ( m Kπ ) is the running resonance width.Therefore, the signal density with the LASS parametrization for the low-mass Kπ S-wave becomes S (Φ) = (cid:88) ξ =1 , − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) H LASS0 A LASS0 ξ + (cid:88) K ∗(cid:48) (cid:88) λ = − , , A K ∗(cid:48) λξ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (9) Model including exotic Z resonances For the decay B → KZ ( → J/ψπ ), J/ψ → µ + µ − where the Z can either be a Z (4200), and/or a Z (4430), or any other exotic charged state, the angle-dependentamplitude is given as A Zλ (cid:48) ξ (Φ) = H Zλ (cid:48) A Z (cid:16) m J/ψπ + (cid:17) d J ( Z )0 λ (cid:48) ( θ Z ) e iλ (cid:48) ˜ ϕ d λ (cid:48) ξ (˜ θ J/ψ ) e iξα (10)where, • J ( Z ) is the spin of the Z resonance, we consider only 1 + spin-parity of the Z s as per Belle’s result [9]; • λ (cid:48) is the helicity of the J/ψ (quantisation axis parallel to the π momentum inthe J/ψ rest frame); • ξ is the helicity of the µ + µ − system; • H Zλ (cid:48) is the complex helicity amplitude for the decay via the intermediate Z ; • d J ( Z )0 λ (cid:48) ( θ Z ) and d λ (cid:48) ξ (˜ θ J/ψ ) are the Wigner small-d functions; • θ Z is the Z helicity angle, i.e. the angle between K and π momenta in the Z rest frame; • ˜ θ J/ψ is the
J/ψ helicity angle, i.e. the angle between µ and π momenta in the J/ψ rest frame; • ˜ ϕ is the angle between the ( µ + , µ − ) and ( K, π ) planes in the
J/ψ rest frame; • α is the angle between the ( µ + , π ) and ( µ + , Kπ ) planes in the J/ψ rest frame.The amplitudes for different λ (cid:48) values are related by parity conservation: H Zλ (cid:48) = − P ( Z )( − J ( Z ) H Z − λ (cid:48) (11)After inclusion of the Z component, the signal density function of Eq. (5) becomes S (Φ) = (cid:88) ξ =1 , − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) K ∗ (cid:88) λ = − , , A K ∗ λξ + (cid:88) Z (cid:88) λ (cid:48) = − , , A Zλ (cid:48) ξ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (12)The signal density function of the charge conjugate decay, identified through thecharge of the K (or π ) differs only in the sign of ϕ . The implementation of thismodel takes into account this switching of sign and also allows for a possible flavourmis-tagging (typically a few %). For the full fit model with ten K ∗ s and two Z sas well as considering the floating masses and widths for some of the resonances,the total number of free parameters in the 4D probability density function (PDF) ur et al. Page 6 of 16 can exceed 60. The large number of free parameters coupled with a complex PDF,which requires many internal mathematical operations to be executed at each stepof the UML fit, poses a real computational challenge.
Fit model implementation in GooFit
The amplitude analysis fit model described in the previous section is implementedusing the novel GPU based GooFit [16, 17] framework because the standard RooFittool designed to perform such UML fits requires excessive processing time to con-verge. GooFit is an under-development open source analysis tool, used in the HEPapplications for parameters estimation, which interfaces ROOT to the CUDA par-allel computing platform on NVIDIA GPUs [18]. GPU-accelerated computing en-hances application performances by offloading a sequence of elementary but com-putationally intensive operations to the GPU to be processed in parallel, while theremaining code still runs on the CPUs. GooFit acts as an interface between theMINUIT [19] minimization algorithm and the GPU, which allows any PDF to beevaluated in parallel over a huge amount of data. Fit parameters are estimated ateach negative-log-likelihood (NLL) minimization step on the host side (CPU) whilethe PDF/NLL is evaluated on the device side (GPU). Since GooFit is still a limitedopen source tool, being developed by the users themselves according to their specificneeds, significant sections needed for our fit implementation have been either newlyencoded or adapted starting from the existing classes and methods.
Timing comparison
The computing capabilities of GPUs versus CPUs are tested by generating andfitting three sets, each comprising 10,000 Monte Carlo (MC) events (pseudo-experiments) of increasing complexity (number of K ∗ s) of the fit model previouslydescribed. The fitter implemented in ROOT/RooFit is run on an Intel Xeon clusterwith 24 CPUs whereas the GooFit version is run on NVIDIA Tesla K40 GPU with2880 CUDA cores. As the timing test models are for the demonstration purposeonly, they are much less complex than the full model required for the analysis, andprocess a smaller number of events than that expected from a collider experiment.As shown in Fig. 2, it becomes almost impossible to run the fitter on CPUs withinany reasonable timescale when the number of fit parameters is increased. The GPU-based GooFit application provides a striking speed-up in performance with respectto the CPU-based RooFit application. The latter gets so time-expensive that it canbecome unreliable once the full number of parameters is adopted in the fit model. Fit validation
The validation is performed by generating, through MC technique, a distributionaccording to the fit model and then checking, as a result of fitting the same dis-tribution, that the best estimates of parameters returned by the fit are consistentwith their input values.
Validation with the K ∗ -only model A pseudo-data sample of one million events is generated with the ten K ∗ s mentionedin Section 2 with their masses and widths fixed to the nominal values. The helicity ur et al. Page 7 of 16
Figure 2
Comparison of time required by RooFit (CPU-based) and GooFit (GPU+CPU based)fitter frameworks to fit three data sets of 10,000 pseudo-experiments, each generated and fittedaccording to models of increasing complexity in terms of the number of K ∗ components. Figure 3
Projection of m Kπ spectrum of the 4D dataset (black points with error bars) generatedaccording to an ideal signal model. The fit result (red points with error bars) is superimposedalong with the individual signal components corresponding to the different K ∗ s. The post-fitvalues of the helicity amplitude parameters and fit fractions for each component are also displayed.ur et al. Page 8 of 16
Figure 4
Projection of (from top to bottom) m J/ψπ , cos θ J/ψπ , and ϕ of the 4D datasetgenerated according to an ideal signal model (black points with error bars). The fit result (redpoints with error bars) is superimposed along with the individual fit components corresponding todifferent K ∗ s. amplitude parameters for each of these resonances are fixed to the values obtainedby Belle [9].The m Kπ projection of the fit to the generated dataset is shown in Fig. 3 andthe other projections are in Fig. 4. The consistency of the post-fit values of the freeparameters is checked by comparing the pull distributions (normalised residuals)with their generated values as shown in Figs. 5 and 6. ur et al. Page 9 of 16 a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ _ a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m A m p li t u d e s Amplitude Values
Gen valuesPost-fit values a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ _ a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m = p o s t f i t g e n e rr o r Amplitude Pull Values
Figure 5
Comparison of generated and post-fit values of the amplitude parameters (above) andthe corresponding pull distribution (below) obtained from a fit to events generated with all theten K ∗ resonances. The green lines define a ± σ band. As the exact contribution of each resonance to the total signal cannot be preciselyevaluated due to interference effects, an approximate measure is provided throughthe fit fractions. The fit fraction of the j -th resonance R j is given by F F j = (cid:82) Ω | A R j (Φ) | d(cid:126)x (cid:82) Ω | S (Φ) | d(cid:126)x (13)where Ω is the four dimensional domain for the set of variables Φ (Eq. (1)) and S (Φ)is the signal function defined in Eq. (12). The numerator of Eq. (13) is obtainedby setting to zero all the other helicity amplitudes at post-fit level. The sum of allthe fit fractions is not constrained to 100% as a consequence of the non-unitarity ofthe model which stems from the constructive and destructive effects of interferencebetween the resonances. Sensitivity of the fitter to Z contributions Validation exercises are performed for a) the K ∗ -only model but with the LASS line-shape used for the S-wave, and b) model with all ten K ∗ s together with Z (4200) and Z (4430) resonances. The mass, width, and helicity amplitudes of the Z resonancesare fixed to the values obtained by Belle [9]. It is found that the post-fit valuesof parameters are consistent with the ones used for generation in both cases. Thefit fractions of the Z -components are found to be small (about a few percent) as ur et al. Page 10 of 16 b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ _ b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m P h a s e s [ r a d ] Phase Values
Gen valuesPost-fit values b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ _ b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m = p o s t f i t g e n e rr o r Phase Pull Values
Figure 6
Comparison of generated and post-fit values of the phase parameters (above) and thecorresponding pull distribution (below) obtained from a fit to events generated with all the ten K ∗ resonances. The green lines define a ± σ band. expected from the Belle results. This confirms that the fitter is capable of correctlydetecting Z contributions even if they are relatively small.Since the Z contributions are expected to be small, we need to ensure that thefitter does not artificially generate Z peaks owing to statistical fluctuations or al-ternative parametrizations of K ∗ signals such as the LASS lineshape. Pseudo-datais generated with only ten K ∗ s and fitted with a [ten K ∗ s + Z (4200) + Z (4430)]model. The fit fraction for both Z (4200) and Z (4430) are found to be 0.01%. FromFigs. 7 and 8, it can be seen that the post-fit helicity amplitude values for the K ∗ sare close to their generated values indicating that the contribution of the Z s areindeed consistent with zero.Similarly, another set of pseudo-data was generated with all K ∗ s (with LASS forthe S-wave) and fitted with an “all K ∗ s (with LASS) + Z (4200) + Z (4430)” model.The fit fractions for Z (4200) and Z (4430) are found to be 0.002% and 0.003%,respectively. Similar to the earlier test, the post-fit helicity amplitude values for the K ∗ s are found to be close to their generated values signifying that the contributionof the Z s are again consistent with zero. Applicability to real-life use cases
An accurate representation of real data from collider experiments would requirethe inclusion of detection efficiency and background contamination. Keeping that ur et al. Page 11 of 16 a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ _ a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ Z _ _ a _ Z _ _ p a _ Z _ _ a _ Z _ _ p A m p li t u d e s Amplitude Values
Gen valuesPost-fit values a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ _ a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ Z _ _ a _ Z _ _ p a _ Z _ _ a _ Z _ _ p = p o s t f i t g e n e rr o r Amplitude Pull Values
Figure 7
Comparison of generated and post-fit values of the amplitude parameters (above) andthe corresponding pull distribution (below) obtained when a dataset generated with ten K ∗ s isfitted with the [ten K ∗ s + Z (4200) + Z (4430) ] model. The green lines define a ± σ band. Thepulls for the Z components are not defined because the Z s are not present in the generationmodel. in mind, the fit framework is developed in such a way that the efficiency and back-ground models of suitable dimensions can be easily included in the form of analyticalfunctions or binned templates. Generic shapes for efficiency (Fig. 9) and background(Fig. 10) in the form of 2D Bernstein polynomials are adopted to test the effective-ness of the fitter with efficiency and background included. Each of the 4D efficiencyand background shapes is passed into the fitter as 2D (mass variables) ×
2D (angu-lar variables) histograms since the masses and the angles are expected to be fully(or at least mostly) uncorrelated.Typically the background levels found in dedicated B-physics experiments (likeBelle and LHCb) are of the order of a few percent [9]. For this test, the fraction is setto a higher value keeping in mind general purpose detectors like CMS and ATLASthat may record signals with slightly less purity due to the absence of dedicatedhadron identification systems. One million pseudo-data are generated and fittedwith a model including all ten K ∗ s, two Z s as well as the relative efficiency andbackground parametrizations. The relative efficiencies are used to weight the signalmodel, whereas the background is added with a fixed coefficient (equal to [1-signalpurity] which in this study is assumed to be 15%).From Figs. 11 and 12, it can be seen that the post-fit helicity amplitude valuesfor the K ∗ s and the Z s are close to their generated values. The fit fractions of Z (4200) (3.49%) and Z (4200) (1.17%) are found to be a few % as expected from ur et al. Page 12 of 16 b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ _ b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ Z _ _ b _ Z _ _ p b _ Z _ _ b _ Z _ _ p P h a s e s [ r a d ] Phase Values
Gen valuesPost-fit values b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ _ b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ Z _ _ b _ Z _ _ p b _ Z _ _ b _ Z _ _ p = p o s t f i t g e n e rr o r Phase Pull Values
Figure 8
Comparison of generated and post-fit values of the phase parameters (above) and thecorresponding pull distribution (below) obtained when a dataset generated with ten K ∗ s is fittedwith the [ten K ∗ s + Z (4200) + Z (4430) ] model. The green lines define a ± σ band. The pullsfor the Z components are not defined because the Z s are not present in the generation model. p m(K3.23.43.63.844.24.44.64.8 ) [ G e V ] py m ( J / - - y J/ q cos(3 - - - f Figure 9
Simulated template of the relative reconstruction efficiency for the scatter plots of thetwo angular variables (right) and the two mass variables (left) of the decay. In the latter, the 2Dkinematic boundary reflects the decay kinematics. the Belle result [9]. The almost identical post-fit values of the parameters from boththe ideal-world case as well as the real-life case indicate that the fitter can producereliable results while taking into account the detection efficiency and backgroundcontributions.Lastly, the fitter though designed for a 4D amplitude analysis of a pseudoscalardecaying into a vector and two pseudoscalars can be easily adapted to other typesof decays with higher or lower dimensions, occurring in flavour physics studies. ur et al. Page 13 of 16
Figure 10
Simulated background template for the scatter plots of the two angular variables (right) and the two mass variables (left) of the decay. In the latter, the 2D kinematic boundaryreflects the decay kinematics. The z -axis values are arbitrary. a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ _ a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ Z _ _ a _ Z _ _ p a _ Z _ _ a _ Z _ _ p A m p li t u d e s Amplitude Values
Gen valuesPost-fit values a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ a _ K _ _ p a _ K _ _ m a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ _ a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ K _ _ _ a _ K _ _ _ p a _ K _ _ _ m a _ Z _ _ a _ Z _ _ p a _ Z _ _ a _ Z _ _ p = p o s t f i t g e n e rr o r Amplitude Pull Values
Figure 11
Comparison of generated and post-fit values of the amplitude parameters (above) andthe corresponding pull distribution (below) obtained from a fit to events generated with ten K ∗ s+ Z (4200) + Z (4430) model including efficiency and 15% background contribution. The greenlines define a ± σ band.ur et al. Page 14 of 16 b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ _ b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ Z _ _ b _ Z _ _ p b _ Z _ _ b _ Z _ _ p P h a s e s [ r a d ] Phase Values
Gen valuesPost-fit values b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ b _ K _ _ p b _ K _ _ m b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ _ b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ K _ _ _ b _ K _ _ _ p b _ K _ _ _ m b _ Z _ _ b _ Z _ _ p b _ Z _ _ b _ Z _ _ p = p o s t f i t g e n e rr o r Phase Pull Values
Figure 12
Comparison of generated and post-fit values of the phase parameters (above) and thecorresponding pull distribution (below) obtained from a fit to events generated with ten K ∗ s + Z (4200) + Z (4430) model including relative efficiency and 15% background contribution. Thegreen lines define a ± σ band.ur et al. Page 15 of 16
Summary
Searches for exotic multiquark states in collider experiments require complex mul-tidimensional analyses involving several (or even hundreds of) thousands of eventsthat demand considerable computational resources. Conventional CPU-based tech-niques may fall short to meet these ever increasing demands. In this study, usingthe helicity formalism, a four-dimensional amplitude analysis framework for an un-binned maximum-likelihood fit has been implemented. The fitting framework hasbeen developed using the novel GPU based GooFit, which is an under-developmentopen source analysis tool used in the HEP applications for parameters’ estimation,interfacing ROOT to the CUDA parallel computing platform on NVIDIA GPUs.It has been shown that the choice to use GooFit and accelerated performancesprovided by GPUs is crucial to carry out these extreme fits.The fit model has been validated by generating and fitting Monte Carlo pseudo-experiments under different model conditions and assumptions with the known setof K ∗ resonances. Since the low mass S-wave K ∗ (800) is not yet satisfactorily de-scribed by a Breit-Wigner amplitude, the alternative LASS parametrization hasbeen implemented on GooFit and thoroughly tested. The fitter is equipped withthe capability of handling relative detection efficiency and background contamina-tion. Finally, the possible contribution of the exotic Z states has been calculatedand incorporated within the fitter framework with reasonable robustness to allowfor testing with any combination of their spin and parity ( J P C ) values as well aswithout constraints. It is hoped that this kind of fitter implemented within theGooFit framework will considerably augment the capabilities of collider experi-ments in searches and measurements in the field of exotic hadron spectroscopy andbeyond.
List of abbreviations
GPU: graphics processing unit; QCD: quantum chromodynamics; UML: unbinned maximum likelihood; HEP:high-energy physics; BW: Breit-Wigner; PDF: probability density function; NLL: negative log likelihood; MC: MonteCarlo.
Availability of data and materials
The datasets generated and analysed during the current study are available from the corresponding author onreasonable request.
Competing interests
The authors declare that they have no competing interests.
Author’s contributions
All authors have equal contribution. They have read and approved the manuscript.
Acknowledgements
We acknowledge Prof. Alexis Pompili (Universit´a degli Studi di Bari and I.N.F.N. - Sezione di Bari) and Prof. GaganMohanty (Tata Institute of Fundamental Research, Mumbai) for their guidance and constant support as well ashelpful comments on this article. We are grateful to the ReCas-Bari [20] Data Center for the valuable supportreceived while using its HPC resources.
Author details Tata Institute of Fundamental Research, Mumbai, India. Universit´a degli Studi di Bari and I.N.F.N. - Sezione diBari, Italy.
References
1. Gross, D.J., Wilczek, F.: Ultraviolet Behavior of Nonabelian Gauge Theories. Phys. Rev. Lett. , 1343–1346(1973)2. Politzer, H.D.: Reliable Perturbative Results for Strong Interactions? Phys. Rev. Lett. , 1346–1349 (1973)3. Gross, D.J., Wilczek, F.: Asymptotically Free Gauge Theories - I. Phys. Rev. D8 , 3633–3652 (1973)4. Politzer, H.D.: Asymptotic Freedom: An Approach to Strong Interactions. Phys. Rept. , 129–180 (1974) ur et al. Page 16 of 16
5. Olsen, S.L., Skwarnicki, T., Zieminska, D.: Nonstandard heavy mesons and baryons: Experimental evidence.Rev. Mod. Phys. , 015003 (2018)6. Ali, A., Maiani, L., Polosa, A.D.: Multiquark Hadrons. Cambridge University Press, Cambridge (2019)7. Brun, R., Rademakers, F.: ROOT: An object oriented data analysis framework. Nucl. Instrum. Meth. A ,81–86 (1997)8. Verkerke, W., Kirkby, D.P.: The RooFit toolkit for data modeling. eConf C0303241 , 007 (2003)9. Chilikin, K., et al. : Observation of a new charged charmoniumlike state in B → J/ψK − π + decays. Phys.Rev. D , 112009 (2014)10. Dalitz, R.H.: On the analysis of τ -meson data and the nature of the τ -meson. Phil. Mag. Ser.7, 1068–1080(1953)11. Chilikin, K., et al. : Experimental constraints on the spin and parity of the Z (4430) + . Phys. Rev. D (7),074026 (2013)12. Blatt, J.M., Weisskopf, V.F.: Theoretical Nuclear Physics. Springer, New York, NY (1979)13. Richman, J.D.: An Experimenter’s Guide to the Helicity Formalism. CALT-68-1148 (1984)14. Tanabashi, M., et al. : Review of particle physics. Phys. Rev. D , 030001 (2018)15. Aston, D., et al. : A Study of K − π + Scattering in the Reaction K − p → K − π + n at 11 GeV/c. Nucl. Phys. B , 493 (1988)16. R. Andreassen, B. T. Meadows, M. de Silva, M. D. Sokoloff and K. Tomko: Goofit: A library for massivelyparallelising maximum-likelihood fits. Journal of Physics: Conference Series (5), 052003 (2014)17. Schreiner, H., et al. : GooFit 2.0. J. Phys. Conf. Ser. (4), 042014 (2018)18. Nickolls, J., Buck, I., Garland, M., Skadron, K.: Scalable parallel programming with cuda. ACM Queue (2008)19. James, F., Roos, M.: Minuit – a system for function minimization and analysis of the parameter errors andcorrelations. Comput. Phys. Commun. , 343 (1975)20. ReCas Is a project financed by the Italian MIUR., 343 (1975)20. ReCas Is a project financed by the Italian MIUR.