Search for Beyond Standard Model Physics at BaBar Citation Li, Yunxuan (2022) Search for Beyond Standard Model Physics at BaBar. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/tz7n-d662. https://resolver.caltech.edu/CaltechTHESIS:05232022-144829107 Abstract This thesis reports searches for beyond Standard Model physics in e + e − collisions in two directions. We report the first search for a dark matter bound state. The existence of dark matter bound states could arise in a simple dark sector model in which a dark photon ( A ′ ) is light enough to generate an attractive force between dark fermions. We report herein a search for a J P C = 1 − − darkonium state, the Υ D , produced in the reaction e + e − → γ Υ D , Υ D → A ′ A ′ A ′ , where the dark photons subsequently decay into pairs of leptons or pions, using data collected with the BaBar detector. No significant signal is observed, and we derive limits on the γ − A ′ kinetic mixing ( ϵ ) as a function of the dark sector coupling constant for 0.001 < m A ′ < 3.16 GeV and 0.05 < m Υ D < 9.5 GeV. Bounds on the mixing strength ϵ down to 5 × 10 − 5 − 10 − 3 are set for a large fraction of the parameter space. We also report a measurement of R ( D ) = B ( B → D̄ τ ν̄ τ )/ B ( B → D̄ ℓ ν̄ ℓ ) and R ( D * ) = B ( B → D̄ * τ ν̄ τ )/ B ( B → D̄ * ℓ ν̄ ℓ ) , where ℓ refers to either an electron or muon. We select samples by reconstructing tag-side B mesons in semileptonic decays and signal-side τ in a purely leptonic decay. Using data collected with the BaBar detector, we measure R ( D ) = 0.316 ± 0.062(stat) ± 0.019(syst) and R ( D * ) = 0.226 ± 0.022(stat) ± 0.012(syst) , which agree with the Standard Model expectations by 0.26 σ and 1.10 σ , respectively. Taken together, the results are in agreement with the Standard Model within 1.51 σ level. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: Beyond Standard Model Physics, Electron-Positron Collider, BaBar Collaboration, High Energy Physics, Dark Matter, B Meson Decays Degree Grantor: California Institute of Technology Division: Physics, Mathematics and Astronomy Major Option: Physics Minor Option: Computer Science Thesis Availability: Public (worldwide access) Research Advisor(s): Hitlin, David G. Thesis Committee: Porter, Frank C. (chair) Hitlin, David G. Patterson, Ryan B. Wise, Mark B. Echenard, Bertrand Defense Date: 20 May 2022 Non-Caltech Author Email: yxli528 (AT) gmail.com Record Number: CaltechTHESIS:05232022-144829107 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:05232022-144829107 DOI: 10.7907/tz7n-d662 ORCID: Author ORCID Li, Yunxuan 0000-0001-9004-2471 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 14598 Collection: CaltechTHESIS Deposited By: Yunxuan Li Deposited On: 27 May 2022 23:21 Last Modified: 04 Aug 2022 23:21 Thesis Files PDF - Final Version See Usage Policy. 16MB Repository Staff Only: item control page

Search for Beyond Standard Model Physics at BABAR Thesis by Yunxuan Li In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy CALIFORNIA INSTITUTE OF TECHNOLOGY Pasadena, California 2022 Defended May 20, 2022 ii © 2022 Yunxuan Li ORCID: 0000-0001-9004-2471 All rights reserved iii ACKNOWLEDGEMENTS I would like to express my deepest and sincerest gratitude to my advisors Prof. David G. Hitlin and Prof. Frank C. Porter. They introduced me to the field of experimental high energy physics, created a great research environment to allow me to try new things, and taught me how to be a great physicist as role models. I want to thank David for all the insightful discussions. During my PhD study, I was always impressed by his physics intuition generating new ideas, understanding complicated results, and getting insights from the numbers. It always enlightened me to think from first principles when I struggled with complex projects. I also want to thank him for sharing lots of past experiences in the high energy physics area. These cherished experiences give me a picture of how this field was developed and a sense of what I was really doing in a historical perspective. I also want to thank Frank for his continuous and detailed support on the thesis. His expertise and academic rigor in physics and statistics encouraged me to question authenticity and validate each hypothesis before taking them for granted. I also want to thank him for allowing me to test my immature and even crazy ideas. I wish to thank Prof. Bertrand Echenard for initiating me into the area of dark matter and spending a great deal of time and effort discussing the project, diagnosing the analysis strategy, debugging the programs, and writing papers. I acknowledge his patience and hand-on support for my growth toward becoming a good researcher, and I really enjoyed the time we worked together. I thank my group members Daniel Chao, Jae Hong Kim, Dexu Lin, James Y. Oyang, Léo Borrel, Sophie Middleton, Jason M. Trevor, and Viktor A. Shcherbatyuk for the help and encouragement during my Ph.D. journey. I thank all the members in the BABAR Collaboration, especially Fabio Anulli, Gerard Bonneaud, and Janis McKenna. Thank you for the great opportunities to let me represent the collaboration and present my research in conferences. I want to thank my good friend Prof. Yiqiu Ma, for not only his immeasurable support on my career development as a senior, but also his open-mindedness to share in-depth thoughts as a friend. I want to thank my good friends Baoyi Chen and Hao Xie. Thank you for your time spent with me on my random and sometimes naive ideas, and your companionship along the journey. I also want to thank Tao Liang, Ruofei Shen, Chenshuo Yue, Zilong Chen, Xingsheng Luan, Xiaolin Mao, iv Junlong Kou, Xiang Li, Boqiang Shen, Yu-An Chen, and Yichuan Song ... for all the fun times together. I am extremely grateful to my parents and sister. Thank you for your love and unconditional support. I really thank my wife Yicheng Wang. Thank your for being in my life and always being there for me. v ABSTRACT This thesis reports searches for beyond Standard Model physics in e+ e− collisions in two directions. We report the first search for a dark matter bound state. The existence of dark matter bound states could arise in a simple dark sector model in which a dark photon (A0) is light enough to generate an attractive force between dark fermions. We report herein a search for a J PC = 1−− darkonium state, the ΥD , produced in the reaction e+ e− → γΥD,ΥD → A0 A0 A0, where the dark photons subsequently decay into pairs of leptons or pions, using data collected with the BABAR detector. No significant signal is observed, and we derive limits on the γ − A0 kinetic mixing () as a function of the dark sector coupling constant for 0.001 < m A0 < 3.16 GeV and 0.05 < mΥD < 9.5 GeV. Bounds on the mixing strength  down to 5 × 10−5 − 10−3 are set for a large fraction of the parameter space. We also report a measurement of R(D) = B(B → D̄τ ν̄τ )/B(B → D̄` ν̄` ) and R(D∗ ) = B(B → D̄∗ τ ν̄τ )/B(B → D̄∗ ` ν̄` ), where ` refers to either an electron or muon. We select samples by reconstructing tag-side B mesons in semileptonic decays and signal-side τ in a purely leptonic decay. Using data collected with the BABAR detector, we measure R(D) = 0.316 ± 0.062(stat) ± 0.019(syst) and R(D∗ ) = 0.226 ± 0.022(stat) ± 0.012(syst), which agree with the Standard Model expectations by 0.26σ and 1.10σ, respectively. Taken together, the results are in agreement with the Standard Model within 1.51σ level. vi TABLE OF CONTENTS Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Illustrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii PART I Introduction 1 Chapter I: The Standard Model of Particle Physics . . . . . . . . . . . . . . . Chapter II: Beyond Standard Model Physics . . . . . . . . . . . . . . . . . . 2 6 PART II PEP-II and BABAR 8 Chapter III: The PEP-II Accelerator . . . . . . . . . . . . . . . . . . . . . . Chapter IV: The BABAR Detector . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Silicon Vertex Tracker (SVT) . . . . . . . . . . . . . . . . . . . . . 4.2 Drift Chamber (DCH) . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Detector of Internally Reflected Cherenkov Light (DIRC) . . . . . . 4.4 Electromagnetic Calorimeter (EMC) . . . . . . . . . . . . . . . . . 4.5 Instrumented Flux Return (IFR) . . . . . . . . . . . . . . . . . . . . Chapter V: Data Record and Luminosity . . . . . . . . . . . . . . . . . . . . 10 11 11 11 12 13 14 16 PART III Search for Dark Matter at BABAR 18 Chapter VI: Search for Darkonium in Electron-Positron Collisions . . . . . . 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Dark Matter Bound States Model . . . . . . . . . . . . . . . . . . . 6.3 Dataset and Signal Monte Carlo Simulations . . . . . . . . . . . . . 6.4 Candidate Reconstruction . . . . . . . . . . . . . . . . . . . . . . . 6.5 Event Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Signal Modeling and Efficiency . . . . . . . . . . . . . . . . . . . . 6.7 Systematic Uncertainties . . . . . . . . . . . . . . . . . . . . . . . . 6.8 Signal Significance Estimation . . . . . . . . . . . . . . . . . . . . . 6.9 Signal Extraction and Upper Limit of Coupling . . . . . . . . . . . . 6.10 Summary and Outlook . . . . . . . . . . . . . . . . . . . . . . . . . 19 19 21 25 25 27 37 43 46 51 54 vii PART IV Precise Measurement of Semileptonic B Meson Decays 57 Chapter VII: Measurement of R(D) and R(D∗ ) Using Semileptonic Tagging and Leptonic τ Decays . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 7.2 Analysis Strategy Overview . . . . . . . . . . . . . . . . . . . . . . 61 7.3 Simulation Samples . . . . . . . . . . . . . . . . . . . . . . . . . . 64 7.4 Event Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . 65 7.5 Signal Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 7.6 Signal Extraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 7.7 Systematic Uncertainties . . . . . . . . . . . . . . . . . . . . . . . . 91 7.8 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 PART V Conclusion 108 Chapter VIII: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 Bibliography 110 Appendix A: Deep Adversarial Network for Systematic Uncertainty Reduction 115 viii LIST OF ILLUSTRATIONS Number Page 1.1 The Standard Model of particle physics. The picture is taken from [1]. 5 3.1 The design of the PEP-II accelerator. . . . . . . . . . . . . . . . . . 10 4.1 The cross-sectional view of the silicon vertex tracker design (upper) along the beam axis; (bottom) orthogonal to the beam axis. . . . . . . 12 4.2 (upper) The side view of drift chamber; (bottom) The cell layout of the drift chamber. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 4.3 The geometry of DIRC. . . . . . . . . . . . . . . . . . . . . . . . . 14 4.4 The side view of the electromagnetic calorimeter barrel and forward endcap. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4.5 The overview of the IFR system. . . . . . . . . . . . . . . . . . . . . 15 5.1 The integrated luminosity over time. . . . . . . . . . . . . . . . . . . 16 6.1 The dark sector scenario under study: The dark sector has an additional U(1)0 group and the corresponding dark photon A0 (γ0 in the picture) plays a role as intermediate messenger between SM and dark sector (hidden sector). . . . . . . . . . . . . . . . . . . . . . . . . . 20 6.2 Constraints on the dark photon parameter space from a re-interpretation of the BABAR dark Higgsstrahlung searches. The solid purple curve corresponds to the current BABAR limit for the parameters αD = 0.5, m χ = 3.5 GeV. The figure is taken from [13]. . . . . . . . . . . . . . 22 6.3 Diagram for ηD and ΥD production and decay at BABAR . The figure is taken from [13]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 6.4 Dark photon branching ratios to specific final states as a function of the dark photon mass. The two peaks around 0.78 GeV and 1.0 GeV for q q̄ correspond to the ω and φ resonances. . . . . . . . . . . . . . 24 6.5 Monte Carlo-generated ΥD and dark photon masses. . . . . . . . . . 26 6.6 Illustration of angle-related physical variables used in the MVA. A_helicity (upper left), A_angle (upper right), A_dihedral (bottom left), and A_poslepton_helicity (bottom right). . . . . . . . . . . . . 31 6.7 Signal and background density distributions of input variables. . . . . 32 6.8 Signal and background density distributions of input variables. . . . . 33 ix 6.9 6.10 6.11 6.12 6.13 6.14 6.15 6.16 6.17 6.18 6.19 6.20 6.21 6.22 Linear correlation matrix of variables used for MVA training for signal MC (upper) and Run3 data (bottom). . . . . . . . . . . . . . . 34 Illustration of the stacked Random Forest model. Figure taken from [20]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 The distribution of the classifier scores for each event category for the data (markers) and signal Monte Carlo (solid lines) samples, for the prompt dark photon decays. The MC simulations are arbitrarily normalized. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 The distribution of the classifier scores for data (markers) and signal MC simulations (solid lines) for dark photon lifetimes corresponding to (top) cτA0 = 0.1mm, (middle) cτA0 = 1mm, and (bottom) cτA0 = 10mm. The MC simulations are arbitrarily normalized. . . . . . . . . 37 Fits to the classification score distributions of data with a function 2 of the form p(x) = ea·x +b·x+c , for C0 (upper), C1 (middle), and C2 (bottom) category of events. . . . . . . . . . . . . . . . . . . . . . . 38 The (mΥD , m A0 ) distribution for events passing all selection criteria for prompt dark photon decays. . . . . . . . . . . . . . . . . . . . . 39 The (mΥD , m A0 ) mass distribution of event candidates passing all selection criteria for the datasets optimized for each dark photon lifetime. 39 Example of CBF fits to the ΥD and A0 mass spectrum (mΥD = 8.0 GeV and m A0 = 0.5 GeV). . . . . . . . . . . . . . . . . . . . . . . . 41 (Left) The ΥD and (right) the dark photon mass resolution as a function of mΥD and m A0 for each category of events. . . . . . . . . . . . . . . 42 The acceptance, selection efficiency, and signal efficiency as a function of mΥD and m A0 for each category of events. . . . . . . . . . . . 43 Signal efficiency (including the branching fractions) as a function of mΥD and m A0 . The two horizontal bands around m A0 = 0.78 GeV and 1.0 GeV correspond to the ω and φ resonances. . . . . . . . . . . . . . . . . 44 The total systematic uncertainties as a function of mΥD and m A0 for combined systematic uncertainties. . . . . . . . . . . . . . . . . . . 46 Distribution of the ml_score for the C2 category for data. The sideband region is shown as a red rectangle. The optimal cut on the classifier score is shown as a solid line. The width of signal region is the δ we set to contain at least 500 events. . . . . . . . . . . . . . . . 47 The distributions of sideband events for the (upper) C0 , (middle) C1 , and (bottom) C2 categories of events. . . . . . . . . . . . . . . . . . 49 x 6.23 6.24 6.25 6.26 6.27 6.28 7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 7.10 7.11 Illustration of the background estimate method for a given mass hypothesis in the mΥD − m A0 space. The dataset displayed here is sideband data for C0 category of events. . . . . . . . . . . . . . . . The distribution of the maximum number of signal events after scanning the parameter space 0 < mΥD < 10 GeV, 0 < m A0 < 3 GeV. . . The 90% CL upper limits on the e+ e− → γΥD cross section for prompt dark photon decays. . . . . . . . . . . . . . . . . . . . . . . The 90% CL upper limits on the e+ e− → γΥD cross section for dark photon lifetimes corresponding to (top left) cτA0 = 0.1 mm, (top right) cτA0 = 1 mm, and (bottom) cτA0 = 10 mm. . . . . . . . . . . . The 90% CL upper limits on the kinetic mixing  2 as a function of the ΥD mass, mΥD , and dark photon mass, m A0 , assuming (top left) αD = 0.1, (top right) αD = 0.3, (middle left) αD = 0.5, (middle right) αD = 0.7, (bottom left) αD = 0.9, and (bottom right) αD = 1.1. The 90% CL upper limits on the kinetic mixing ε for (top) various ΥD masses assuming αD = 0.5 and (bottom) various αD values assuming mΥD = 9 GeV together with current constraints (gray area). . . . . . . Feynman diagram of B → D(∗) τν decay. . . . . . . . . . . . . . . . Visualization of previous measurements of R(D) and R(D∗ ). The plot is from HFLAV [22]. . . . . . . . . . . . . . . . . . . . . . . . The distributions of reconstructed D meson mass for (upper left) D+ decays without π 0 , (upper right) D0 decays without π 0 , (bottom left) D+ decays with π 0 , and (bottom right) D0 decays with π 0 . . . . . . . The distribution of mass difference between reconstructed D∗ meson and its daughter D meson for (left) D∗+ → D0 π + , (middle) D∗+ → D+ π 0 , and (right) D∗0 → D0 π 0 decays. . . . . . . . . . . . . . . . . Histograms of variables used for the C1 classifier. . . . . . . . . . . . Histograms of variables used for the C1 classifier. . . . . . . . . . . . Histograms of variables used for the C1 classifier. . . . . . . . . . . . Area under the ROC curve for BDT classifiers with different numbers of trees. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Importance of each variable for learning the C1 classifier. . . . . . . . z1 distribution for signal, normalization, D∗∗ lν, B B̄ combinatorial, and continuum events. . . . . . . . . . . . . . . . . . . . . . . . . . Histograms of variables used for the C2 classifier. . . . . . . . . . . . 50 51 53 54 55 56 58 61 69 70 74 75 76 77 78 78 79 xi 7.12 7.13 7.14 7.15 7.16 7.17 7.18 7.19 7.20 7.21 7.22 7.23 7.24 7.25 7.26 8.1 The area under the ROC curve for BDT classifiers with different number of trees. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 Importance of each variable for learning the C2 classifier. . . . . . . . 81 z2 distribution for signal, normalization, D∗∗ lν, B B̄ combinatorial, and continuum events. . . . . . . . . . . . . . . . . . . . . . . . . . 81 KDE learned densities for event components in the D+ l subset. . . . 84 KDE learned densities for event components in the D0 l subset. . . . . 85 KDE learned densities for event components in the D∗+ l subset. . . . 86 KDE learned densities for event components in the D∗0 l subset. . . . 87 Comparison of z1 score distributions of the data with the projections of fit results for (upper left) D+ l, (upper right) D0 l, (bottom left) D∗+ l, and (bottom right) D∗0 l subsets. . . . . . . . . . . . . . . . . . 88 Comparison of z2 score distributions of the data with the projections of fit results for (upper left) D+ l, (upper right) D0 l, (bottom left) D∗+ l, and (bottom right) D∗0 l subsets. . . . . . . . . . . . . . . . . . 89 MC/data comparison of z1 and z2 score distributions in normalization enriched region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 MC/data comparison of z1 and z2 score distributions in backgrounds enriched region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 (upper) The current mass ranges for excited D meson states. (bottom left) The R(D∗∗ ) for different excited D meson masses. (bottom right) The Φ(B → D∗∗ lν) for different excited D meson masses. . . . . . . 99 The π 0 efficiency correction values in dependence of the π 0 momentum. The figure is from [48]. . . . . . . . . . . . . . . . . . . . . . . 101 z1 and z2 distribution of sidebands for generic MC and sideband data. The difference between data and histogram indicates the discrepancy between MC and data. . . . . . . . . . . . . . . . . . . . . . . . . . 102 The comparison of preliminary results with previous measurments for R(D(∗) ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Architecture of unsupervised domain adaptation. The figure is taken from [49]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 xii LIST OF TABLES Number Page 5.1 Integrated luminosity of on-peak and off-peak data collected for each run. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 6.1 Number of signal Monte Carlo events generated for each category of final states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 6.2 PID selection rules of different final states. . . . . . . . . . . . . . . 30 6.3 Algorithm to find converged signal efficiency, cross-section upper limits, and kinetic mixing strength upper limits. . . . . . . . . . . . . 53 7.1 Previous results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 7.2 Definition of event types in the B B̄ system. . . . . . . . . . . . . . . 62 7.3 Cross sections used to convert the sizes generic simulated data to the equivalent on-peak dataset. . . . . . . . . . . . . . . . . . . . . . . . 64 7.4 Generic simulated data. Multiplier is the factor by which the size of the corresponding on-peak dataset exceeds that of the given simulated dataset. More simulated datasets are usually generated to better study the decay characteristics. . . . . . . . . . . . . . . . . . . . . . . . . 65 7.5 Simulated signal data. . . . . . . . . . . . . . . . . . . . . . . . . . 66 7.6 Simulated normalization data. . . . . . . . . . . . . . . . . . . . . . 67 7.7 D meson decay modes used in the analysis. . . . . . . . . . . . . . . 68 7.8 Proportions of each component after reconstruction, evaluated using the generic MC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 7.9 Fit results for the yields of all the components in four subsets. . . . . 89 7.10 Evaluated relative systematic uncertainties from the B → Dlν and B → D∗ lν form factors. . . . . . . . . . . . . . . . . . . . . . . . . 92 7.11 Evaluated relative systematic uncertainties from the B → D∗∗ lν form factors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 7.12 Decay modes fluctuated to evaluate the B → D(∗) lν branching fractions on generic MC. The world average values are from [34]. . . . . 94 7.13 Evaluated relative systematic uncertainties from B → D(∗) lν branching fractions on MC. . . . . . . . . . . . . . . . . . . . . . . . . . . 94 xiii 7.14 7.15 7.16 7.17 7.18 7.19 7.20 7.21 7.22 7.23 7.24 7.25 7.26 B → D∗∗ (D(∗) π)(l/τ)ν decays and method used to estimate the systematic uncertainties. "MC" means the uncertainties can be directly estimated using MC sample. "Phase Space" means we do not have the corresponding MC samples, therefore a phase-space-based estimation is applied. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Decay modes fluctuated to evaluate the resonant B → D∗∗ (1P)(l/τ)ν branching fractions on generic MC. The setup values are from [34, 45]. 96 Evaluated relative systematic uncertainties from resonant B → D∗∗ (1P)(l/τ)ν branching fractions. . . . . . . . . . . . . . . . . . . . . . . . . . . 96 MC samples used for assessing the gap between inclusive and the sum of exclusive B → Xc `ν branching fractions. . . . . . . . . . . . 97 Evaluated relative systematic uncertainties from the resonant B → D∗∗ (2S)lν branching fractions. . . . . . . . . . . . . . . . . . . . . . 97 Evaluated relative systematic uncertainties from non-resonant B → D(∗) πlν branching fractions. . . . . . . . . . . . . . . . . . . . . . . 97 Evaluated relative systematic uncertainties from resonant B → D∗∗ (2S)τν and non-resonant B → D(∗) πτν branching fractions. . . . . . . . . . 99 Evaluated relative systematic uncertainties from the D → Kππ branching fractions. . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Evaluated relative systematic uncertainties from B(Υ(4S)) decays to charged or neutral B mesons. . . . . . . . . . . . . . . . . . . . . . . 100 Evaluated relative systematic uncertainties from the calibration factor on the B B̄ background shape difference between generic MC and data. 103 Evaluated relative systematic uncertainties from the calibration factor on the B B̄ background shape difference between generic MC and data. 103 Summary of evaluated uncertainties (preliminary). . . . . . . . . . . 105 Efficiencies for signal and normalization events. . . . . . . . . . . . 106 Part I Introduction 1 2 Chapter 1 THE STANDARD MODEL OF PARTICLE PHYSICS The ultimate goal of particle physics is to understand the fundamental law of matter in the universe, which is made of elementary particles. In other words, elementary particles are the building blocks of the universe and are thought to have no internal structure. According to the idea of reductionism, if we understand the characteristics of elementary particles and how they interact with each other, we are able to fully explain and link together all physical phenomena in the universe. The corresponding theory to model elementary particles and their interactions is often referred as the Theory of Everything (TOE). Developed during the twentieth century, the standard model (SM) is a quantum field theory to describe three of four fundamental forces: electromagnetism, strong, and weak interactions. Although it does not include gravity, the standard model is currently our best attempt to unify the fundamental interactions together, and it has been extensively tested by a series of experiments. Some famous experiments include discovery of Z boson at CERN, W boson at CERN and Fermilab, J/ψ meson at SLAC and BNL, and Higgs boson at CERN. The Lagrangian of the standard model is 1 1 1 LSM = − Bµν B µν − Tr[Wµν W µν ] − Tr[G µν G µν ] 4 4 4 / + h.c. + i ψ̄ Dψ Õ + ψi yi j ψ j φ + h.c. (1.1) i,j=e,µ,τ + |D µ φ| 2 − V(φ). The first term is the gauge term, where Bµν, Wµν and G µν are the field tensors for the U(1), SU(2), and SU(3) groups. Using the weak mixing angle θ w , the A µ, Z µ, W µ± bosons are mixtures of Wµν and G µν components: A µ = W11µ sin θ w + Bµ cos θ w, Bµ = A µ cos θ w − Z µ sin θ w Z µ = W11µ cos θ w − Bµ sin θ w, W11µ = −W22µ = A µ sin θ w + Z µ cos θ w √ √ ∗ W µ+ = W µ−∗ = W12µ / 2, W12µ = W21µ = 2W µ+ . (1.2) 3 The second term involves kinetic terms of fermions and how fermions interact with gauge fields. ψ is a vector including all fermion fields. Explicitly, the Lagrangian can be re-written as / + h.c. LF,F−G =i ψ̄ Dψ ! ν L =(ν̄L , l¯L )σ̃ µiD µ + l¯R σ µiD µ lR + ν̄R σ µiD µ νR + h.c. lL ! u L +(ū L , d¯L )σ̃ µiD µ + ū R σ µiD µ u R + d¯R σ µiD µ dR + h.c., dL (1.3) where l = (e, µ, τ), ν = (νe, ν µ, ντ ), u = (u, c, t), and d = (d, s, b) are all threecomponent generation indices. The σ are the Pauli matrices satisfying " ! ! ! !# 1 0 0 1 0 −i 1 0 σµ = , , , 0 1 1 0 i 0 0 −1 (1.4) σ̃ µ = [σ 0, −σ 1, −σ 2, −σ 3 ]. D µ is covariant derivative combining electromagnetic, weak, and strong interactions in the form of D µ = ∂µ − ig0Y Bµ − igW µaT a − igs G aµ t a where g0, g, gs are the coupling strengths of the hypercharge, weak and strong interaction. Y,T a, t a are the corresponding hypercharge operator, SU(2) generator, and SU(3) generator. The interaction between fermions and gauge fields can be further decomposed to charged and neutral currents. Taking leptons as an example, its charged current part of the Lagrangian is g µ L LCC = − √ ( jW,L W µ + h.c.) 2 2 µ jW,L = ν̄l γ µ (1 − γ 5 )l. (1.5) Its neutral current includes the QED Lagrangian as well as a weak neutral current: γ L LNC = L L + L LZ γ µ µ L L = −e jγ,L A µ, jγ,L = −ēγ µ e g µ µ ¯ µ (cl − cl γ 5 )l. L LZ = − j Z,L Z µ, j Z,L = ν̄l γ µ (cVν − c νA γ 5 )νl + lγ V A 2 cos θ w (1.6) The last two terms are related to the Higgs mechanism. They describe the Higgs field and how the Higgs field gives mass to gauge and fermion fields. For the 4 Higgs-gauge interaction, since the hypercharge Y = +1, we have D µ = ∂µ + igWµ · σ 1 + ig0 B µ · . 2 2 So the Higgs-gauge term can be written as 1 1 LH,G = (∂µ H)(∂ µ H) − λv 2 H 2 − λvH 3 − λH 4 2 4 v 2 g 2 + −µ v 2 g2 + Wµ W + Z Zµ 2 µ 4 8 cos θ w 2 vg + −µ v 2 g2 v 2 g 2 + −µ 2 v 2 g2 µ + Wµ W H + Z Z H + W W H + Z µ Z µ H 2, µ µ 2 4 4 cos θ w 2 8 cos θ w 2 (1.7) from which we can get the mass of the Higgs as well as the mass of W and Z bosons: p mH = 2λv 2 vg mW = (1.8) 2 vg . mZ = 2 cos θ w Since the mass term of fermion fields mψ̄ψ = m(ψ̄L ψR +ψ̄R ψL ) is not gauge invariant, the way to give mass to fermions is via Yukawa coupling. Once diagonalized, it can be written as Õ LY ukawa = ψi yi j ψ j φ i,j=e,µ,τ   yl v ¯ yl ¯ = √ ll − √ hll 2  2 yu yu v + √ ūu − √ hūu 2  2 yd v ¯ yd ¯ + √ dd − √ h dd, 2 2 (1.9) so the mass of leptons and quarks are yi v mi = √l , i = e, µ, τ 2 i y v mi = √u , i = u, c, t 2 i y v mi = √d , i = d, s, b. 2 (1.10) 5 Figure 1.1: The Standard Model of particle physics. The picture is taken from [1]. Fig. 1.1 shows all the elementary particles in the standard model. In conclusion, the standard model has 17 types of elementary particles. 12 of them are fermions, which are the building blocks of the matter in our universe. The two type of fermions, quark and lepton, are distinguished by whether they carry color charge. In addition, we have 4 gauge bosons which act as force carriers: photon for electromagnetism, gluon for strong interaction, and W and Z boson for weak interaction. We also have the Higgs boson to explain the mass of the other particles. 6 Chapter 2 BEYOND STANDARD MODEL PHYSICS Although the standard model has been tested by many experiments, it is still an incomplete theory due to the lack of explanation of the following phenomena and theoretical puzzles: • Gravity: the standard model does not explain gravity, which is currently best described with general relativity and explains the gravitational force due to curved space-time from massive objects. • Dark matter: The existence of dark matter has been supported by many cosmological and astrophysical observations, and it is measured to be 45 times that of ordinary matter in our universe. However, particles and interactions in the SM only explain the physics of ordinary matter, leaving dark matter unexplained. • Dark energy: Dark energy is around 70% of the entire energy in our universe, and it is also not explained in the SM. • Neutrino masses: the neutrinos in the SM are massless, however, neutrino oscillation experiments show that they do have mass. • Hierarchy problem: The SM introduces a Higgs field to explain particle masses, however, the mass parameter in the Higgs field needs to be finetuned to cancel quantum correlations. However, relying on a numerical cancellation to explain the large discrepancy between weak force and gravity is uncomfortable for many physicists. • Strong CP problem: The CP symmetry in the strong interaction sector is allowed to be broken in theory, while no experiments observe such violation. Therefore, one of the main goals of particle physics is to search for and understand the physics beyond the SM. To search for BSM physics, particle physics has three categories of approaches, referred to as the energy frontier, the cosmic frontier, and the intensity frontier. While the energy and cosmic frontiers often provide a direct search for the production of new particles, the intensity frontier provides indirect 7 probes of new physics with intense beam sources and super-sensitive detectors. With much better sensitivity, intensity frontier experiments can probe for new physics at mass scales of hundreds to thousands of TeV. Part II PEP-II and BABAR 8 9 The primary goal of the experiment is to study CP violation in the decays of B mesons. With its high luminosity, a sensitive measurement of the CKM matrix element Vub can be made, which is a fundamental parameter in the Standard Model. This high-luminosity experiment is designed such that electrons and positrons collide at the Υ(4S) resonance, which then decays to BB pairs more than 95% of the time. The B mesons subsequently decay to other particles, which can be detected or reconstructed from specific sub-systems of the BABAR detector. At the hardware level, charged or neutral final-state particles generate hits or clusters in the detector. The optimized particle identification algorithms (PID) or cluster-finding algorithms are applied to identify final-state particles and reconstruct the entire decay process. Besides the Υ(4S) resonance, the BABAR experiment also operated at the Υ(3S) and Υ(2S) resonances, as well as off-resonance to better understand the continuum backgrounds. The BABAR experiment can be also applied to study a range of other physics including, but not limited to, other B meson decays, the physics of charm and τ leptons, and two-photon physics. 10 Chapter 3 THE PEP-II ACCELERATOR Figure 3.1: The design of the PEP-II accelerator. The PEP-II accelerator is designed to deliver the B mesons to the experiment. It is an e+ e− collider operating at the Υ(4S) resonance of 10.58 GeV. The two storage rings are asymmetric with 9 GeV for the electron beam and 3.1 GeV for the positron beam, which is a key design feature of PEP-II, as shown in Fig. 3.1. This design results in B mesons with significant momenta in the lab frame, and its decay length can be measured. The PEP-II system includes the following four major subsystems: Injector, HighEnergy Ring (HER), Low-Energy Ring (LER), and Interaction Region (IR) [2]. The Injector includes the extraction and transport lines from the SLAC two-mile linac. The HER is responsible for the 9 GeV electron beam. It is originated from PEP ring but with new vacuum, RF, feedback, and diagnostics systems. The LER is a new ring mounted above the 9-GeV HER, responsible for the 3.1 GeV positron ring. The IR is designed to focus the beam spots, bring two beams into head-on collision, and separate them cleanly. 11 Chapter 4 THE BABAR DETECTOR In order to achieve the sensitivity required, the BABAR detector needs to have maximum possible acceptance in the center-of-mass (CM) system, excellent vertex resolution to measure B meson decay time, and great discrimination between charged particles to tag B meson decays. The detector mainly includes six components: Silicon Vertex Tracker (SVT), Drift Chamber (DCH), Detector of Internally Reflected Cherenkov light (DIRC), Caesium Iodide Electromagnetic Calorimeter (EMC), a 1.5 T superconducting coil, and Instrumented Flux Return (IFR). The details of design [3] are summarized below. 4.1 Silicon Vertex Tracker (SVT) The main goal of the SVT is to reconstruct the decay vertices of the two B mesons to determine the time between two decays, and measure time-dependent CP asymmetries in B meson decays. The SVT is also responsible for the detection of low-momentum charged particles which do not reach the drift chamber. As shown in Fig. 4.1, the SVT includes five concentric cylindrical layers of doublesided silicon detectors, divided in azimuth into modules. The inner three layers are barrel-style structures with six detector modules. The outer two layers have 16 and 18 detector modules, respectively, and are employed with a new arch structure so that detectors are electrically connected across an angle. The silicon vertex detectors use double-sided silicon strip detectors coupled with polysilicon bias resistors. The p+ and n+ strips are fabricated on around 300-µmthick high-resistivity silicon. The front-end signal processing is performed by ICs mounted on hybrid circuits, fabricated in a radiation-hard CMOS technology. 4.2 Drift Chamber (DCH) The drift chamber is the main tracking system of the BABAR detector used to reconstruct tracks with transverse momentum larger than 100 MeV/c. The performance goal of the DCH includes a spatial resolution better than 140 µm, a PID for low momentum tracks by dE/dx with resolution of 7%, and a resolution of σpt ≈ 0.3% × pt for momentum above 1 GeV/c. 3.3 The Silicon Vertex Tracker 83 12 84 An Introduction to the BABAR Experiment Figure 3-4. Detail of the inner part of the apparatus, showing the cross-sectional view of the silicon vertex tracker in a plane containing the beam axis. y Each module is divided into two electrically separated forward and backward half-modules, the shorter (if they are not exactly symmetrical) always being the forward one. In order to follow the designed structure of the fourteen half-modules, detectors with six different wafer geometries are needed. x The inner sides of the detectors have strips oriented perpendicular to the beam direction to measure the coordinate ( -side), whereas the outer sides, with longitudinal strips, allow the coordinate measurement ( -side). The read-out electronics will be placed outside the active area, the -side strips will be connected to them with flexible Upilex fanout circuits glued to the inner faces of the half-modules. In the two outer modules, the number of strips exceeds the number of readout channels, so some fraction of the strips is “ganged”, i.e., two strips are connected to the same readout channel. The “ganging” is performed by the fanout circuits. The strips are daisy-chained Figure 3-5. Layout of the BAR silicon vertexoftracker. view incapacitance the plane of between detectors, resulting in aBAtotal strip length to 26Cross-sectional and a maximum Figure 4.1: The cross-sectional view of theupsilicon vertex tracker design (upper) orthogonal to the beam axis. 35 pF. Also, for the -side, a short fanout extension is needed to connect the ends of the strips to alongthethe beam axis; (bottom) orthogonal to the beam axis. readout electronics. outer detector contain of external done after theThe wire bonding. between The modules are supported on Half-modules fromlayers two toisfour detectors. connections two adjacent detectors Kevlar mounted to the cones in the and wire backward A carbonon the ribs -side, and with the end fanout on located both sides, areforward made with bonds.directions. The bending of the fiber space which isare mounted on the magnets The side viewframe and supports layout the of assembly, drift chamber shown inbending Fig. 4.2. Theusing driftkinematic chamber is mounts. The space frame is needed to allow some motion of the magnet during the assembly R EPORT BABAR P HYSICS ORKSHOP enda 280-cm-long madeOFofTHE aluminum. TheWforward procedure. cylinder. Its flat endplates are BABAR silicon vertex will have 340region silicon detectors in total, covering an areawith of about plateThe is made thinner in tracker the acceptance of the detector compared the rear 1 , with a total of 150.000 readout channels. Details of the silicon vertex tracker parameters endplate. Thein drift cellshere arethearranged in 10 superlayers of 4forlayers each, scheduled are shown Table 3-4: intrinsic point resolution is calculated tracks incident at 90 , assuming a signal-to-noise ratio, . using the pattern AUVAUVAUVA for axial (A) and stereo (U, V) superlayers. The of Silicon Microstrip Detectors 4.3 3.3.3 Detector Internally Reflected Cherenkov Light (DIRC) As shown in Fig. 4.3, the DIRC is designed to provide excellent kaon identification Silicon vertex detectors use double-sided silicon strip detectors, AC coupled with polysilicon bias 0 → + π − . It resistors [5]. between They have pthestrips on the p-side and nmodes strips on n K strips to distinguish two-body decay B0the→opposite, π + π − n-side. and BThe on the n-side must be interleaved by p implants (p-stops). They are fabricated on (300 30) is also designed to identify muons in the low-momentum range. The DIRC system includes 144OFlong, barsWof synthetic quartz arranged in a 12-sided polygonal R EPORT THE Bstraight ABAR P HYSICS ORKSHOP the design of the detector is optimized to reduce the material in the forward end. The forward endplate is made thinner (12 mm) in the acceptance region of the detector compared to the rear endplate (24 mm), and all the electronics is mounted on the rear endplate. The inner cylinder is made of 1 mm beryllium, which corresponds to 0.28% radiation lengths (X ). The outer cylinder consists of 2 layers of carbon fiber on a Nomex core, corresponding to 1.5% X . 13 324 1015 1749 35 90 202 551 973 IP 68 17.19 469 An Introduction to the236 BABAR Experiment 1618 Figure 3-7. Side view of the BABAR drift chamber. The dimensions are in mm. R EPORT OF THE BABAR P HYSICS WORKSHOP Figure 4.2: (upper) The side view of drift chamber; (bottom) The cell layout of the drift chamber. Figure 3-8. Cell layout in the BABAR drift chamber. barrel. 4.4 R EPORT OF THE BABAR P HYSICS WORKSHOP Electromagnetic Calorimeter (EMC) The primary goal of EMC is to measure the energy of photon and electrons with precise energy and angular resolution. When a photon or an electron travels through the crystal, it interacts with the thallium-doped cesium iodide (CsI(Tl)) and triggers 14 Figure 4.3: The geometry of DIRC. an electromagnetic cascade and deposits its energy to create scintillation light. The photons are then reflected inside the crystal and collected with the PIN diodes. As shown in Fig. 4.4, the calorimeter includes a cylindrical barrel section and a forward conic endcap. It is located asymmetrically about the interaction point with the inner radius extending 112.7 cm in the backward direction and 180.9 cm in the forward direction. The barrel crystals covers angle between −0.775 ≤ cos θ ≤ 0.892 in the laboratory frame, corresponding to −0.916 ≤ cos θ ≤ 0.715 in the c.m.s. frame. 4.5 Instrumented Flux Return (IFR) The IFR shown in Fig. 4.5 is designed to return the flux of the 1.5T magnet and to provide muon identification and neutral hadron detection. The system consists of a central part (Barrel) and two plugs (End Caps) to cover the angle from 300 mrad in the forward direction to 400 mrad in the backward direction. The resistive plate chambers (RPC) are constructed of bakelite sheets separated by thick polycarbonate spacers with cylindrical symmetry and fiberglass frame. 15 98 An Introduction to the BABAR Experiment 112.7 cm 180.9 cm 140.8° e– 26.9° IP 15.8° e+ 91 cm 7–96 8184A1 Figure 4.4: The side view of the electromagnetic calorimeter barrel and forward Figure 3-11. The EMC layout: Side view showing dimensions (in mm) of the calorimeter barrel endcap. and forward endcap. The barrel consists of 5760 CsI(Tl) crystals, arranged in 48 polar-angle ( ) rows of distinct crystal sizes, each having 120 identical crystals in azimuthal angle ( ). The crystals are grouped in 280 modules, each spanning crystals in and , respectively (except for the most-backward module, which has only 6 crystals in ). The modules are made from 300 thick carbon fiber composite (CFC) material and held from the rear by an aluminum strongback. They are mounted in an aluminum support cylinder, which in turn is fixed to the coil cryostat. By supporting the crystals from the rear, minimal material is placed in front of them. The EMC front material consists of two 1 -thick cylinders of aluminum, separated by foam, which provide a gas seal and rf shielding. Additional front material due to the liquid radioactive source calibration system consists of the equivalent of another 3 of aluminum. (All thicknesses are quoted at perpendicular incidence.) Cooling, cables, and services are located at the back of each module and thus do not add to the inactive-materials budget. The total amount of external material in front of the calorimeter barrel (endcap) is about 0.25 (0.20 ) at a polar angle of 90 (20 ), the major contributors being the DIRC (drift chamber endplate) with about 0.14 in both cases. Areas of particular concern are the forward and backward end of the barrel, where the material grows to 0.5 and 0.39 , respectively. R EPORT OF THE BABAR P HYSICS WORKSHOP Figure 4.5: The overview of the IFR system. 16 Chapter 5 DATA RECORD AND LUMINOSITY The BABAR experiment operated from 1999 to 2008. An integrated luminosity of 531 fb−1 was collected mostly at the Υ(4S) resonance, but also at the Υ(3S), Υ(2S), as well as off-resonance. Fig. 5.1 shows the integrated luminosity over time. The luminosity for each run is shown in Table 5.1. Figure 5.1: The integrated luminosity over time. 17 Run period Run1 Run2 Run3 Run4 Run5 Run6 Run7 Run8 Total Luminosity (fb−1 ) on-peak off-peak 20.4 2.6 61.3 6.7 32.3 2.4 99.6 10.0 133.2 14.3 78.3 7.8 27.9 2.6 13.6 1.4 465.7 47.9 Data sample on-peak Υ(4S) Υ(4S) Υ(4S) Υ(4S) Υ(4S) Υ(4S) Υ(3S) Υ(2S) Table 5.1: Integrated luminosity of on-peak and off-peak data collected for each run. Part III Search for Dark Matter at BABAR 18 19 Chapter 6 SEARCH FOR DARKONIUM IN ELECTRON-POSITRON COLLISIONS 6.1 Introduction The existence of dark matter is overwhelmingly supported by astrophysical and cosmological observations, and it is one of the most important tasks of contemporary physics to understand its nature and properties. The concept of dark matter can be traced back to 1880s, when Lord Kelvin estimated the number of dark bodies in the Milky Way from the observed velocity dispersion of the stars orbiting the center of the galaxy. From his results, he concluded that many of our stars may be dark bodies. In 1933, Fritz Zwicky reached a similar conclusion by applying the virial theorem to the Coma Cluster in an attempt to estimate its mass, finding a result 400 times larger than the visible mass. During the last decades, much evidence for dark matter has been accumulated from measurements of galaxy rotation curves, velocity dispersions, galaxy clusters, and gravitational lensing. Many dark matter candidates have been proposed during the last decades. A longleading paradigm is in the form of Weakly Interacting Massive Particles (WIMPs), as the dark matter self-annihilation cross section is roughly consistent with a weakscale massive particle interacting via the electroweak force. Many beyond Standard Model theories, such as Supersymmetry [4] or models containing extra dimensions, include WIMP-like candidates. However, the absence of observation of WIMPs so far has motivated the exploration of alternative models, such as light dark matter and the existence of a dark sector. The dark sector (also known as hidden sectors) is a collection of particles that do not interact directly with SM particles, i.e. the dark sector fields are singlets under the SM gauge groups. They are motivated by many BSM theories, for example, in string theory constructions and type-II compactifications. The dark sector has its own symmetries, which could be arbitrarily complex. If the dark sector contains an extra U(1) symmetry, the dark and visible sectors are indirectly connected via the so-called vector portal, namely they can interact via kinetic mixing between the dark gauge field (A0µ ) and the SM hypercharge field F µν : 20  F µν,Y Fµν,D Fµν,D = ∂µ A0ν − ∂ν A0µ, with  denoting the mixing strength, as shown in Fig. 6.1. Values of mixing strength in the 10−12 − 10−3 range have been predicted in the literature [5, 6, 7]. The kinetic mixing term can be removed by redefining the vector field A0µ → A0µ, A µ → A µ +  A0µ, so that SM fermions pick up a small dark charge ∼  e leading to "milli-charged" interaction between the two sectors. The dark photon can be massive, and its mass can arise via the Higgs or the Stueckelberg mechanisms. Its mass range is predicted to be in the MeV-GeV range [8, 9, 10, 11], though much smaller (sub-eV) masses are also possible [12]. Besides the vector portal, there are a few other indirect interactions that can connect the dark sector to the SM, including the scalar and Fermion portals. In this analysis we focus on the vector portal. Dark sector DM !(1)′ #!" !"!" &′# SM !" 3 ×!" 2 ×"(1) Figure 6.1: The dark sector scenario under study: The dark sector has an additional U(1)0 group and the corresponding dark photon A0 (γ0 in the picture) plays a role as intermediate messenger between SM and dark sector (hidden sector). 21 6.2 Dark Matter Bound States Model The possibility of dark sector bound states was first proposed by An et al.[13]. A specific model contains a Dirac dark matter field ( χ) charged under an additional U(1) gauge group in the dark sector, with the corresponding vector boson acting as a mediator between the dark sector and SM via kinetic mixing. The mass of dark matter particles is O(few GeV). The Lagrangian for this model is L =L SM + χ̄iγ µ (∂µ − igD A0µ ) χ − m χ χ̄ χ 1  1 − A0µν A0µν − Fµν A0µν + m2A0 A0µ A0µ, 4 2 2 (6.1) where A0µ is the vector boson mediator with field strength A0µν , Fµν is the SM g2 hypercharge field,  is the kinetic mixing strength, and αD = 4πD is the strength of the dark electromagnetic interaction. When gD is sufficiently large, the force between the dark fermions mediated becomes attractive, resulting in the formation of dark matter bound states ( χ χ̄). The two lowest (1S) bound states, 1 S0 and 3 S1 , are respectively denoted ηD and ΥD in analogy with the SM. The critical conditions for the existence of stable bound states has been determined numerically [14], 1.68m A0 ≤ αD m χ . If we assume αD = 0.5 and m χ = 3.5 GeV, the dark photon must be lighter than 1 GeV, as shown in Fig. 6.2. The value of αD has to be large enough to produce a bound state; m χ higher than 3.5 GeV has been excluded by other measurements [15, 16]. The quantum numbers of the two lowest (1S) bound states suggest the following production mechanisms at electron-positron colliders: e+ e− → ηD + A0 e+ e− → ΥD + γ. These production processes are mediated by a mixed γ − A0 propagator, as shown in Fig. 6.3. The ηD (ΥD ) decays into 2 (3) pairs of dark photons, leading to multi-lepton and/or multi-quark final states. 22 SM particles prompt for t above decay containing si muons or pio We turn t tial state rad non-relativist matter field e im t [ · p/ + + e e ! ⌘D +V ; e e ! 𝑚 ⌥!D! (GeV) + ; p+p ! ⌥D +X (2) nihilation (cr We use the Figure 6.2: Constraints on the dark photon parameter space from a re-interpretation The last process represents the direct production of ⌥ FIG. 2. Left: Constraint on the dark photon parameter space from the BA BA D of the BABAR dark Higgsstrahlung searches. The solid purple curve corresponds to function [30] from q q̄Band fusion. production are mediated byisolid purple curve production dark bound and ⌥The The D . figure the current ABAR decay limitAll forof the parameters αDprocesses =states 0.5, m χ⌘=D 3.5 GeV. taken from [13]. parameters ↵D =V0.5, m = 3.5 GeV. The dashed curve shows the future a mixed propagator, as shown in Fig. purple 1. on the m mV plane for the SIDM scenario are shown with 2 = 10 7 andh0d $ $ favored for SIDM solving the galactic small-scale structure problems [3] for ↵D + ! $ $ where and "µ⌥D th is e e ! (⌘D V, ⌥!D ) ! 3V channels are shown in thick purple curves, $ # $ ! ! are shown in# thin blue curves. $Allowed regions are in the arrow R direction. A ⌘D (0) is the + " ! $ the e e ! for ↵the " ! ¯ + 2V channel $are shown in dot-dashed black curves account k D = ! $ ! excluded by CDMSlite [37] and$ LUX [38]. The region mV $. 30 MeV is ruled photon, we od for ↵D = 0.3. tween ⌥D an 𝜖 As discussed in the introduction, sufficiently strong dark interaction strength and light dark photon will result in the formation of dark matter particles ( ¯). The two lowest (1S) bound states, 1 S0 (J P C = 0 + ) and 3 S1 (J P C = 1 ), will be called ⌘D and ⌥D , respectively. The condition for their existence has been determined nu2 merically [26] 2 , 1.68mV < ↵D m , with ↵D = gD /(4⇡). Their quantum numbers suggest the following production mechanisms at colliders: # # ! ! # " " # ! ! " ! ! # ! " # ! " " Figure 6.3: Diagram for ηD and ΥD production and decay at BABAR . The figure is FIG. 1. Diagram for ⌘D and ⌥D production and decay at taken from [13]. B-factories. beams, the most important production channel is from analysis is searching forfusion, production. The⌥ e+ e−. emit an initial state radiation theThis quark-anti-quark q q̄ ! Generalizing In order to obtainΥDthe rate for the process in cal(2), D first (ISR) photon to produce a ΥD state. The ΥD then + decays into three dark photons, we calculate of cross e e section ! ¯Vis with ,¯ culations of [42],the theamplitude production given by and these dark photons subsequently decay to lepton or quark-antiquark pairs. In 2 2 having the same2four2 momentum p (with p = m ), and Z 1 2 X this analysis, we will4⇡ only↵ search A0 → X + X − (X =dx e, µ, π). for 2 D apply the projection operator, = Q pp(n)!⌥D q The ΥD production state ⌧radiation s via initial x is calculated in [13] by scross-section q applying ΥD rest frame, h a non-relativistic1expansion ⇣ ⌧ ⌘ to the dark matter field in⇣the ⌘i ⌧ ⇧ = R (0)(6 p + m ) (6 p m ) , [17] (3) ⌘ relation between ⌘ 5 and⇥ using the element the wave function fq/p (x)fq̄/p(n) + fq̄/p (x)fand , (10)to 32⇡m3 the Dmatrix q/p(n) x x where ⌧ = the m2V⌘/s, fq/p(n) and fq̄/p(n) area the relevant to select state [28]. We find leading-order D bound differential cross for section: structure functions this process, and Qq is the quark charge Unlike de2 B-factories, 2 d in+ units of e. 4⇡↵↵ 2 [R (0)]2 (1 only + cosmuonic ✓) e e !⌘D V D ⌘D 3 Le↵ = current unc refrain fromm In the limit 2 periments in ↵D /2. We ob not cover th d Ve+< m M e 200 ! ⌥D d cos ✓ A Outlook. " and gested ⇥ sensitive to (s same time.4m W strong where ✓self-in is th facilitates th the CM fram 23 derive the effective kinetic mixing term between ΥD and the photon s 1 αD µν L e f f = − D FµνΥD , D = RΥD (0). 2 2m3χ When m A0  αD m χ , the following differential cross section is obtained ! 4m2χ dσe+ e− →γΥD 2πα2  2 D2 ≈ 1− d cos θ s s # " 8s2 (s2 + 16m4χ ) sin2 θ −1 , × (s − 4m χ )2 (s + 4me2 − (s − 4me2 ) cos 2θ)2 (6.2) (6.3) where θ is the angle between γ and the initial e− in the center-of-mass frame. The ΥD particle will subsequently decay to three dark photons; the differential decay rate is calculated by generalizing the approach in [18] to massive dark photon case, 3 [R (0)]2 dΓ(ΥD → A0 A0 A0) 2αD ΥD = dx1 dx2 3πm2χ 39x 8 + 4x 6 F6 − 16x 4 F4 + 32x 2 F2 + 256F0 × 2 , (x − 2x1 )2 (x 2 − 2x2 )2 (x 2 + 2(x1 + x2 − 2))2 (6.4) where x1,2 = E1,2 /m χ , x = m A0 /m χ , and F6 = x12 + (x1 + x2 )(x2 − 2) − 30 F4 = (x12 + x1 x2 − 2x1 )(3x2 − 10) − 10x2 (x2 − 2) − 21 F2 = x14 + 2x13 (x2 − 2) + x12 (x2 (3x2 − 22) + 28)+ 2x1 (x2 − 2)(x2 (x2 − 9) + 12) + x2 (x2 − 2)(x2 (x2 − 2) + 24) + 24 (6.5) F0 = x14 + 2x13 (x2 − 2) + x12 (3x2 (x2 − 3) + 7)+ x1 (x2 − 1)(x2 − 2)(2x2 − 3) + (x2 − 1)2 (x2 (x2 − 2) + 2). We assume that the coupling of dark photons to SM fermions is universal and proportional to their charge. Based on the following constraints, • B(A0 → hadrons)/B(A0 → l + l − ) = σ(e+ e− → hadrons)/σ(e+ e− → l + l − ) = R, Õ • B(A0 → l + l − ) + B(A0 → hadrons) = 1, l=e,µ,τ (6.6) we obtain the dark photon branching ratios into specific final states, as shown in Fig. 6.4. A′ branching ratio 24 1.0 e+ e− µ +µ − π +π − qq 0.8 0.6 0.4 0.2 0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 m (GeV) A′ Figure 6.4: Dark photon branching ratios to specific final states as a function of the dark photon mass. The two peaks around 0.78 GeV and 1.0 GeV for q q̄ correspond to the ω and φ resonances. Integrating (6.3) over θ, and combining with (6.4) and the dark photon branching ratios, we obtain the total cross section of the signal channel e+ e− → γΥD,ΥD → A0 A0 A0, A0 → X + X − (X = e, µ, π), which is a function of  and m A0 only. One complication arises from the dark photon lifetime, which can become sufficiently large to produce displaced decay vertices and affect the mass resolution and signal efficiency. The partial decay width of the dark photon in the case of m A0 > 2me is given by Γ ≈ 31 m A0 α 2 , where m A0 is the dark photon mass, α is the strength of SM electromagnetic interaction,  is the kinetic mixing strength. Therefore, dark photon decay length l in the laboratory frame is l = γvcτ = p 3~c · , m2 α 2 (6.7) which becomes significant when the kinetic mixing and the dark photon mass are small. For instance, if m A0 = 10 MeV and  = 10−4 , the flight length of dark photon is around 0.2 m. The effect of the dark photon flight length becomes non-negligible above O(100 µm), a few times the resolution of the beam spot location. 25 6.3 Dataset and Signal Monte Carlo Simulations Signal MC events are generated using MadGraph, which calculates matrix elements from first principles. Our final states consist of three pairs of leptons or pions, leading to 10 different possibilities. We do not consider the 3π + π − combination, as this channel has much more background and much lower signal-to-noise ratio than the other possibilities. These final states can be divided into three categories based on the number of pion pairs: 0π + π − pair (C0 ), 1π + π − pair (C1 ) and 2π + π − pairs (C2 ). We simulate zero-lifetime signal MC events in the range 0 < mΥD ≤ 10 GeV and 0 < m A0 ≤ 3 GeV. The number of events simulated for each final state is listed in Table 6.1. The dark photon and dark Upsilon mass hypotheses we generated are shown in Fig. 6.5. The range of dark Upsilon mass is [0.05, 0.15, 0.3, 0.5, 1, 1.5, 2, 2.25, 3, 4, 4.5, 5, 6, 7, 7.5, 8, 9, 10] GeV. The range of dark photon mass is [0.01, 0.02, 0.05, 0.1, 0.2, 0.3, 0.47, 0.5, 0.71, 0.75, 1, 1.42, 1.5, 2, 2.38, 2.5, 2.85, 3] GeV. The generated events are passed through the detector response simulation based on GEANT4 [19] and reconstructed using the same software chain as the experimental data. The dark photons are reconstructed using the prompt reco sequence. Signal MC events with non-zero lifetimes are also generated using MadGraph and reconstructed using the same software chain as the experimental data using prompt reco sequence. The dark photon flight lengths we simulated are [0.1, 1, 10] mm. The range of dark Upsilon mass is [ 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] GeV. The range of dark photon mass is [0.02, 0.05, 0.1 , 0.2, 0.5] GeV. This signal MC will be used to evaluate the impact of dark photon lifetime on the signal efficiencies. The background is complex due to multiple lepton and pion pairs in the final states and cannot be simulated accurately by current generators. Instead, we use 5% of the data (Run3 data) as a background sample to optimize our selection criteria and validate the fitting procedure. We assume that the signal component is negligible in the validation data. 6.4 Candidate Reconstruction We search for the reaction e+ e− → γΥD,ΥD → A0 A0 A0, where the dark photon subsequently decays into a lepton or pion pair. We conduct our analysis by first identifying events containing exactly six charged tracks, and we form all possible combinations of three dark photon candidates. The electrons, muons and pion pairs are selected from loose PID algorithms (a track may satisfy several PID 26 Category C0 C1 C2 Final state 3e+ e− 2e+ e− 1µ+ µ− 1e+ e− 2µ+ µ− 3µ+ µ− 2e+ e− 1π + π − 2µ+ µ− 1π + π − + 1e e− 1µ+ µ− 1π + π − 1e+ e− 2π + π − 1µ+ µ− 2π + π − Number of events generated 2040083 1103956 1103956 Table 6.1: Number of signal Monte Carlo events generated for each category of final states. Generated signal MC in (m D , mA ) Υ 3.0 ′ Dark photon mass / GeV 2.5 2.0 1.5 1.0 0.5 0.0 0 2 4 6 8 Dark Upsilon mass / GeV 10 Figure 6.5: Monte Carlo-generated ΥD and dark photon masses. algorithms and be used multiple times with different mass hypotheses). To limit the combinatoric background, we require the presence of at least two tracks loosely identified as either electrons or muons, and the maximum mass difference between the three dark photon candidates must be less than 0.5 GeV. Tighter PID and criteria are applied at a later stage. We then perform a kinematic fit with a beam constraint to each reconstructed ΥD candidate. No energy constraint is applied, to allow for the possibility of initial state radiation or Bremsstrahlung from electrons. The ISR photon kinematics can then be inferred, based on the dark photon and ΥD kinematics. For signal events, the recoil mass should be compatible with the photon hypothesis (∼ 0 GeV), while 27 background events produce a wide distribution. The ISR photon can be observed or not, and this analysis treats both cases simultaneously. If the extrapolated ISR photon is emitted in the EMC acceptance, a corresponding photon must be found. Otherwise, only a small amount of extra neutral energy should be seen. Same-sign A0 candidates are also reconstructed. For our analysis with 6 charged tracks grouped into 3 neutral dark photon pairs, the same-sign A0 candidates are obtained by swapping two particles of the same charge with opposite charge in different pairs. The mass difference between the same-sign candidates tends to be smaller for background events than signal events. The background is mostly combinatorics and the same-sign mass difference is distributed more broadly than signal events, i.e, the probability of having a small same-sign mass difference is larger for background than signal. 6.5 Event Selection To further improve the signal purity, we use a multivariate classifier based on the following variables. • χ2 : The probability of the constrained fit on the ΥD candidate, imposing a beam constraint to the six tracks. Signal events should have small value, as they should be compatible with the beam spot and their tracks can form a common vertex. Background events are spread continuously. • PIDpass: The particle type (e, µ, π) assigned to each track must be compatible with the PID requirements described in Table 6.2. We explored several PID requirements by changing the selectors and the required number of charged tracks, and choose the ones optimizing the signal significance. • massdiff : Maximum mass difference between the three dark photon candidates. Given the reconstructed masses of the three dark photons m1A0, m2A0, m3A0 , we compute massdiff = max( m1A0 − m2A0 , m2A0 − m3A0 , m1A0 − m3A0 ). Signal events should be peaked around 0, while background events are spread continuously. • ISR_costh: The cosine of the polar angle of the reconstructed ISR photon in the laboratory frame. 28 • bestISRPhoton: A categorical feature indicating whether the emitted ISR photon is found in the calorimeter. Its value can be {−2, −1, 0}. Its value is -2 when the reconstructed ISR photon polar angle is outside the detector acceptance. Its value is -1 when reconstructed ISR photon polar angle is within the detector acceptance, but not found in the calorimeter, i.e, there is a missing energy in this event. Its value is 0 when the reconstructed ISR photon polar angle is within the detector acceptance and detected by calorimeter. The ISR photon is considered to be detected if a neutral cluster is found in the calorimeter with an energy within 10% of that reconstructed from the ΥD measurement, and the cosine of the angle between the neutral cluster and reconstructed ISR photon directions smaller than 0.1. • ISR_mass2: Invariant mass of reconstructed ISR photon m2ISR = (pe+ e− − pΥD )2 . Signal candidates should be compatible with the photon hypothesis, peaking around 0, while background candidates span a large range of values. • E_extra: Sum of neutral clusters in the ECAL, excluding the ISR photon candidate and Bremsstrahlung photons. Signal events are expected to have small extra neutral energy, while the background events can have larger values. • A_helicity: The average helicity angle of three dark photons. The helicity angle is defined as the angle between the ΥD flight in CM frame and A0 flight in ΥD rest frame, as illustrated in Fig. 6.6. • A_angle: The average of angles between the dark photons in the ΥD rest frame, as illustrated in Fig. 6.6. • A_dihedral: The average dihedral angles between the dark photons. For two dark photons A01 → X1+ X1− and A02 → X2+ X2− , the dihedral angle between A01 and A02 is the angle between the planes defined by X1+ X1− and X2+ X2− , as illustrated in Fig. 6.6. • A_poslepton_helicity: The average helicity angle of the positive lepton for each dark photon decay, as illustrated in Fig. 6.6. • A_FlightLen_avg: The average flight length of the "zero-lifetime" dark photons. • same_sign_massdiff : For each candidate with three pairs of opposite charged particles, we construct its same sign candidates by swapping two same particles with opposite charge in different pairs. If there is only one type of 29 particle in final states, i.e, in the form of X1+ X1−, X2+ X2−, X3+ X3− where X indicates symbol of particles (e, µ, π) and integer indicates which dark photon it decays from, we have 3 × C32 = 9 possible same sign candidates: X1+ X2+, X1− X2−, X3+ X3− X1+ X2+, X1− X3−, X3+ X2− X1+ X2+, X3− X2−, X3+ X1− X1+ X3+, X2+ X1−, X2− X3− X1+ X3+, X2+ X2−, X1− X3− (6.8) X1+ X3+, X2+ X3−, X1− X2− X2+ X3+, X1+ X1−, X2− X3− X2+ X3+, X1+ X2−, X1− X3− X2+ X3+, X1+ X3−, X1− X2− . If there are two types of particles in final states, i.e, X1+ X1−, X2+ X2−,Y3+Y3− where X and Y are different particles from e, µ, π, we have only one possible candidate: X1+ X2+, X1− X2−,Y3+Y3− . If there are three types of particles in final states, namely e+ e−, µ+ µ−, π + π − , we will consider pions as muons and use the previous rule. The same sign mass candidate for tracks X1 and X2 is calculated as p A0 = p X1 + p X2 m A0 = E A2 0 − p2A0 . (6.9) We choose the minimum same-sign mass difference among all possible combinations for a candidate. The mass difference between the same-sign candidates tends to be smaller for background events than signal events. The background is mostly combinatorics and the same-sign mass difference is distributed more broadly than signal events, i.e, the probability of having a small same-sign mass difference is larger for background than signal. The distributions of these features for all three categories are shown in Fig. 6.7 and 6.8, in which each entry comes from one candidate. Correlations among these features are shown in Fig. 6.9. The background candidates are taken from Run3 data, 30 Final state 3e+ e− 2e+ e− 1µ+ µ− 1e+ e− 2µ+ µ− 3µ+ µ− 2e+ e− 1π + π − 1e+ e− 1µ+ µ− 1π + π − 2µ+ µ− 1π + π − 1e+ e− 2π + π − 1µ+ µ− 2π + π − PID Selection Criteria at least 4 tracks pass tight electron PID and 5 track passes loose electron PID at least 3 tracks pass loose electron PID and 1 track passes loose muon PID at least 1 track passes loose electron PID and 3 tracks pass loose muon PID at least 3 tracks pass tight muon PID and 5 tracks pass loose muon PID at least 3 tracks pass loose electron PID at least 1 track passes tight electron PID and 1 track passes loose muon PID at least 3 tracks pass loose muon PID and 2 tracks pass tight muon PID at least 1 track passes tight electron PID at least 1 track passes loose muon PID Table 6.2: PID selection rules of different final states. while the signal candidates are selected from a mixture of MC samples generated at all different masses. Machine learning classifier Using the features described above, our next step is to choose and train machine learning models for signal/background classification. The dataset we use is a mixture of Run3 data as negative (background) sample and signal MC as positive (signal) sample. We separate the dataset into two disjoint subsets acting as a training dataset and a validation dataset respectively. Models are first trained using the training dataset, while hyper-parameters are selected based on the models’ performance on validation dataset to reduce over-training. We tried both Random Forest (RF) and Gradient Boosting Decision Tree (GBDT) models, and we found that RF has better performance, so we use RF in this analysis. The classification probability of RF models is discrete by construction. For example, a RF classifier containing 100 trees will have score values separated by increments of 0.01, which is insufficient for our purpose. In order to obtain a continuous classification probability and determine the hypothesis testing rejection region, we perform a logistic regression on the RF output, as shown in Fig. 6.10. Specifically, we train the RF with our training dataset first. The outcome of each decision tree in the RF indicates which leaf node our event ends. This information is represented by a binary vector d = [0, 0, 0, ..., 1, ..., 0], 31 e+ γ e+ A' ΥD θh A' γ e− ΥD A_angle A' e− 𝘭 ΥD A' 𝘭 + A_poslepton_helicity − Figure 6.6: Illustration of angle-related physical variables used in the MVA. A_helicity (upper left), A_angle (upper right), A_dihedral (bottom left), and A_poslepton_helicity (bottom right). indicating that the event ended in the i th leaf node. Given a RF with n decision trees, the binary representation of the full RF is a stack of decision trees’ representations [d1, d2, ..., dn ]. To obtain a continuous score function, we fit a Logistic Regression function on top of this binary representation of training samples as our classification model. We denote this model as stacked RF. All of models are trained and validated using the scikit-learn package. We trained stacked RF under a different number of assembled decision trees (n_estimators). Other hyper-parameters such as the criterion (the function to measure the quality of a split) and the min_samples_split (minimum number of samples required to split an internal node), are kept as default values, as they are believed to have been well-tuned. Under the figure of merit defined below, we find that RF with 60 trees offers stable and close to optimal performance, and we select this 0.07 background signal 0.06 0.05 1.0 Frequency Density 32 background signal 0.8 0.6 0.04 0.03 0.4 0.02 0.2 0.01 20 40 60 80 background signal 50 True 100 χ2 Density Density 0.000 False PIDpass background signal 20 40 15 30 10 20 5 10 0.1 0.2 0.3 0 −1.00 −0.75 −0.50 −0.25 0.00 0.4 0.5 massdiff (GeV) 1.0 background signal 0.8 Density Frequency 0 0.0 1.4 0.25 0.50 0.75 1.00 ISR costh background signal 1.2 1.0 0.6 0.8 0.6 0.4 0.4 0.2 0.2 0 -1 3.5 0.0 bestISRPhoton background signal 3.0 Density Density -2 1.75 2 4 6 8 10 12 14 ISR mass square (GeV) background signal 1.50 2.5 1.25 2.0 1.00 1.5 0.75 1.0 0.50 0.5 0.25 0.00 0 2 4 6 8 10 Eextra (GeV) 0.00 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 0.8 1.0 Average helicity of dark photon Figure 6.7: Signal and background density distributions of input variables. 1.4 Density Density 33 background signal background signal 0.8 1.2 1.0 0.6 0.8 0.4 0.6 0.4 0.2 0.2 0.25 0.50 0.75 0.0 −1.00 −0.75 −0.50 −0.25 1.00 Costh of angle between A’ background signal 15 0.4 10 0.2 5 0.00 0.25 0.50 0.75 0.25 0.50 0.75 1.00 0 0.0 1.00 Costh of helicity angle of positive lepton background signal 20 0.6 0.0 −1.00 −0.75 −0.50 −0.25 0.00 Costh of dihedral angle between A’ Density 0.8 0.00 Density Density 0.0 −1.00 −0.75 −0.50 −0.25 0.1 0.2 0.3 0.4 0.5 Average dark photon flight length (mm) background signal 1.2 1.0 0.8 0.6 0.4 0.2 0.00 1 2 3 5 4 same sign massdiff (GeV) Figure 6.8: Signal and background density distributions of input variables. 34 Linear co elation mat ix among MVA va iables χ2 0.8 PIDpass massdiff ISR_costh bestISRPhoton ISR_mass2 E_extra A_helicity A_angle A_dihedral A_poslepton_helicity A_FlightLen_avg same_sign_massdiff 0.4 0.0 −0.4 PIDpass massdiff ISR_costh bestISRPhoton ISR_mass2 E_extra A_helicity A_angle A_dihedral A_poslepton_helicity A_FlightLen_avg same_sign_massdiff χ2 −0.8 Linear co elation mat ix among MVA va iables χ2 0.8 PIDpass massdiff ISR_costh bestISRPhoton ISR_mass2 E_extra A_helicity A_angle A_dihedral A_poslepton_helicity A_FlightLen_avg same_sign_massdiff 0.4 0.0 −0.4 PIDpass massdiff ISR_costh bestISRPhoton ISR_mass2 E_extra A_helicity A_angle A_dihedral A_poslepton_helicity A_FlightLen_avg same_sign_massdiff χ2 −0.8 Figure 6.9: Linear correlation matrix of variables used for MVA training for signal MC (upper) and Run3 data (bottom). 35 Figure 6.10: Illustration of the stacked Random Forest model. Figure taken from [20]. configuration. Selection criteria optimization The distribution of the classifier outputs is shown in Fig. 6.11 and Fig. 6.12 for both prompt and displaced dark photon decays. The optimal criteria on the MVA output (ml_score) is determined by optimizing the average cross-section upper limit over the whole mass range, assuming every observed event is a signal candidate. This approach is clearly conservative, but as it can be seen in Fig. 6.11 and 6.12, the background rises rapidly in the vicinity of the optimal cut, while the signal efficiency varies smoothly. Minimizing the background level (ideally completely suppressing it) therefore has a small effect on the signal sensitivity, but greatly reduces the associated systematic uncertainties and facilitates the determination of the signal significance as well. More specifically, we optimize the classifier by minimizing the figure of merit λ̂(n)  , where λ̂(n) is the 90% CL upper limit on the number of signal events given n observed events in the signal region (assumed to be signal), and  is the signal efficiency averaged over the phase space. We determined n and  as a function of the classifier cut by fitting the classifier output distributions, and integrated the density over the cut value. We apply a binned fit to the ml_score histogram of the data events in the range of 0 ≤ ml_score ≤ 8, by fitting 2 a function of the form f (x) = ea·x +b·x+c . The fit is performed for each category of C0 , C1 and C2 separately. The fitting results are shown in Fig. 6.13. The dashed 36 Entries vertical lines indicates the optimal cuts. The range of 0 ≤ ml_score ≤ 8 ensures that the fit is completely dominated by background events and still offers a good estimation of the data. Then, for a given cut c, we estimate the expected number of ∫ +∞ observed events by integrating the fitting density over the cuts: n = N · c f (x)dx, where N is the total number of data events in the corresponding category. C0 : Data C1 : Data C2 : Data 106 105 C0 : Signal C1 : Signal C2 : Signal 104 103 102 101 100 −20 −15 −10 −5 0 5 10 15 20 Score Figure 6.11: The distribution of the classifier scores for each event category for the data (markers) and signal Monte Carlo (solid lines) samples, for the prompt dark photon decays. The MC simulations are arbitrarily normalized. Best candidate selection The best candidate selection aims to assign each event with one single best matching candidate. If multiple candidates pass the MVA selection for one event, the best candidate is selected based on its final state, according to the following sequence of hypotheses: 6e, 4e2µ, 2e4µ, 6µ, 4e2π, 2e2µ2π, 4µ2π, 2e4π, 2µ4π. The sequence of hypotheses is ordered according to the purity of the final states to minimize mis-identification among channels. Entries 37 105 104 Data, cτ =0.1 mm Data, cτ =1 mm Data, cτ =10 mm Signal, cτ =0.1 mm Signal, cτ =1 mm Signal, cτ =10 mm 103 102 101 100 −30 −20 −10 0 10 20 30 40 Score Figure 6.12: The distribution of the classifier scores for data (markers) and signal MC simulations (solid lines) for dark photon lifetimes corresponding to (top) cτA0 = 0.1mm, (middle) cτA0 = 1mm, and (bottom) cτA0 = 10mm. The MC simulations are arbitrarily normalized. Event selection results For signal with prompt decays, a total of 69 events pass all the selection criteria. The corresponding (mΥD , m A0 ) distribution is shown in Fig. 6.14. Most events correspond to e+ e− → q q̄ process. The events near mΥD ∼ 0.1 GeV and m A0 ∼ 0.05 GeV arise from e+ e− → γγγ events in which all three photons convert to e+ e− pairs. For signal with displaced decays, a total of 56, 33, and 31 events are selected for the cτA0 = 0.1, 1, and 10 mm data sample, respectively. The resulting mass distributions are shown in Fig. 6.15. 6.6 Signal Modeling and Efficiency Mass resolution The A0 and ΥD mass resolutions are used for both determining the size of a given scan box, as well as calculating the scan step size. The mass resolutions are estimated by fitting the A0 (ΥD ) mass distributions with Crystal Ball functions (CBF) for each MC sample. One example of the fits is shown in Fig. 6.16, in which we fit the ΥD and dark photon mass spectrum at mΥD = 8 GeV, m A0 = 0.5 GeV. The CBF captures the energy loss in low tails and gives a better estimation of the (∆mΥ , ∆m A0 ) resolution 38 Data 105 104 103 102 10 1 − 10 −5 0 5 10 15 5 10 15 5 10 15 Data 106 105 104 103 102 10 1 − 10 −5 0 Data 106 105 104 103 102 10 1 − 10 −5 0 Figure 6.13: Fits to the classification score distributions of data with a function of 2 the form p(x) = ea·x +b·x+c , for C0 (upper), C1 (middle), and C2 (bottom) category of events. mA0 (GeV) 39 2.0 C0 C1 C2 1.5 1.0 0.5 0.0 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 mΥD (GeV) mA0 (GeV) Figure 6.14: The (mΥD , m A0 ) distribution for events passing all selection criteria for prompt dark photon decays. cτ = 0.1 mm cτ = 1 mm cτ = 10 mm 0.20 0.15 0.10 0.05 0.00 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 mΥD (GeV) Figure 6.15: The (mΥD , m A0 ) mass distribution of event candidates passing all selection criteria for the datasets optimized for each dark photon lifetime. 40 than a simple Gaussian fit. While the fits are not perfect, they are sufficient for the purpose of estimating the mass resolution. To evaluate the mass resolution as a function of (mΥD , m A0 ), we interpolate the results at known points using a 2-dimensional smooth spline. Smooth spline function aims to fit a polynomial function on data and guarantees the smoothness of the function at the same time. To minimize the fitting uncertainties, we fit each category of events individually. The results are shown in Fig. 6.17. The ΥD mass resolutions are at the level of 20 MeV. The dark photon mass resolutions are at the level of 3 MeV. Signal efficiency The signal efficiency is evaluated by counting the fraction of simulated events passing the selection criteria in the signal window at a given point, normalized to the total number of generated events. The size of the signal window is given by the mass resolutions at that point, namely 4∆mΥ and 4∆m A0 . To evaluate the signal efficiency as a function of (mΥD , m A0 ), we also use a 2dimensional smooth spline technique for interpolation, fitting each category separately. The degrees of the polynomial functions used for interpolation are increased until the model is stable. The results of the signal acceptance, selection efficiency, and signal efficiency are shown in Fig. 6.18. The signal acceptance is low when mΥD and m A0 are low, which stops us from improving the signal efficiency in the low mass region. The MVA prefers to select events in the range of 6 GeV ≤ mΥD ≤ 9 GeV, thus the selection efficiencies are higher in this region. For each signal window, the signal efficiencies from C0, C1 , and C2 are weighted by branching fractions to obtain the total signal efficiency of the signal window. The total signal efficiency is then used for signal extraction. The result is shown in Fig. 6.19. The two horizontal bands around m A0 = 0.78 GeV and 1.0 GeV correspond to the ω and φ resonances. The drop in the branching fraction A0 → X + X −, X = e, µ, π at m A0 = 0.78 GeV and 1.0 GeV is due to the A0 decaying predominantly into π + π − π 0 and K + K − , respectively. Lifetime dependency The signal efficiencies depend on the dark photon lifetime. When the dark photons have non-zero lifetime, the χ2 fit of assuming prompt decay becomes larger and the mass resolution worse. The larger the dark photon lifetime, the smaller the signal 41 ΥD mass C ystal Ball Fit Pe fo mance 10 C ystal Ball Fit M.C. Mass Histog am Density 8 6 4 2 0 6.0 120 Density 100 6.5 7.0 7.5 8.0 8.5 ΥD mass / GeV 9.0 9.5 10.0 Dark Photon mass Crystal Ball Fit Performance Crystal Ball Fit M.C. Mass Histogram 80 60 40 20 0 0.400 0.425 0.450 0.475 0.500 0.525 0.550 0.575 0.600 Dark Photon Mass / GeV Figure 6.16: Example of CBF fits to the ΥD and A0 mass spectrum (mΥD = 8.0 GeV and m A0 = 0.5 GeV). 42 C0 Dark Upsilon Mass Resolution (MeV) 0.5 0.0 2 4 6 ΥD mass / GeV 8 10 C1 Dark Upsilon Mass Resolution (MeV) Dark photon ma / GeV 2.00 1.75 1.50 1.25 1.00 0.75 0.50 0.25 0 4 6 ΥD ma / GeV 8 10 C2 Dark Upsilon Mass Resolution (MeV) 30 Dark photon ma / GeV 2.00 1.75 1.50 1.25 1.00 0.75 0.50 0.25 2 35 30 25 20 15 10 5 25 20 15 10 5 2 4 6 ΥD ma / GeV 8 10 C0 Dark Pho on Mass Resolu ion (MeV) 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 Dark pho on mass / GeV 1.0 2.0 1.5 1.0 0.5 0.0 2.00 1.75 1.50 1.25 1.00 0.75 0.50 0.25 0 2 4 6 ΥD mass / GeV 8 10 C1 Dark Photon Ma Re olution (MeV) Dark photon ma / GeV 1.5 35 30 25 20 15 10 5 2.00 1.75 1.50 1.25 1.00 0.75 0.50 0.25 2 4 6 ΥD ma / GeV 8 10 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 C2 Dark Photon Ma Re olution (MeV) Dark photon ma / GeV Dark pho on mass / GeV 2.0 2 4 6 ΥD ma / GeV 8 10 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 Figure 6.17: (Left) The ΥD and (right) the dark photon mass resolution as a function of mΥD and m A0 for each category of events. efficiency. We generate and reconstruct the signal MC with non-zero lifetime to estimate and interpolate the impact of lifetime on the signal efficiencies. We use the functional form log() = a · log(l) + b to fit the relation between the average signal efficiencies () and the dark photon lifetime (l). We checked that the results obtained for each mass hypothesis are compatible with those obtained by considering the average efficiency, but the later are more robust and were chosen to describe the global dependence on the lifetime. For a given mass hypothesis (mΥD , m A0 ) with prompt signal efficiency 0 , the relation between the signal efficiency and the dark photon lifetime under the mass hypothesis is reasonably described by: log() = a · (log(l) − log(l0 )) + log(0 ), (6.10) 43 Figure 6.18: The acceptance, selection efficiency, and signal efficiency as a function of mΥD and m A0 for each category of events. with l0 = 0.1 mm, and 0 the efficiency for prompt decays. For a dark photon flight lengths above 100 µm, this interpolated relation will be used to correct the signal efficiency and establish the kinetic mixing strength upper limit. When the dark photon flight length is below 100 µm, we use the efficiency determined for prompt decays, as the displaced decay vertices of the dark photon have negligible effect on the signal efficiencies. 6.7 Systematic Uncertainties The following source of systematic uncertainties are considered: • Efficiency MC statistics (σst at ): The uncertainties on the efficiency estimation due to the limited Monte Carlo simulation areÍtaken as a systematic uncertainty. N xi Our estimate of signal efficiency is ˆ = i=1 N , where xi ∼ Bernoulli() indicates whether the i th event is selected as signal, and it follows a Bernoulli distribution with a true signal efficiency  probability of taking 1. Therefore, 44 Figure 6.19: Signal efficiency (including the branching fractions) as a function of mΥD and m A0 . The two horizontal bands around m A0 = 0.78 GeV and 1.0 GeV correspond to the ω and φ resonances. q ˆ ) ˆ the statistical uncertainty of the signal efficiency estimate is σˆ = (1− N , where ˆ is the estimated signal efficiency and N is the number of simulated events. • Efficiency Interpolation (σint ): This is the uncertainty due to the interpolation of the signal efficiency, as a function of mΥD and m A0 . For each M.C. sample, we have both an estimated signal efficiency (est ) and a predicted signal efficiency (eval ). The former was obtained by counting the fraction of simulation events passing the classifiers and within our scanning windows, while the latter is the output of smooth spline functions obtained by removing this M.C. point and refitting the spline interpolation. The interpolation uncertainty at this point is defined as |est − eval |/est . To evaluate the interpolation uncertainty of any point, we fit a smooth spline function on the interpolation uncertainties for each category of events. The average interpolation uncertainty is within the scale of 8%, which indicates that the spline functions provide a good interpolation. 45 • PID: A systematic efficiency uncertainty of 2% per muon, 2% per pion, and 1% per electron is used. The PID uncertainty is added linearly for each final state, and we average the PID uncertainties within each category of events, assuming each final state within category occurs with equal probability. The final uncertainty is 9% for C0 , 10% for C1 and 11% for C2 . • Luminosity: The uncertainty on the luminosity for Run 1-6 is taken to be 0.5%, following the result from BAD 2186. A systematic of 1.2% is taken for the fraction of Run7 data missing B-counting information. The luminosity weighted uncertainty is 0.6%. • Tracking efficiency: The systematic uncertainty is taken to be 0.2% per track, determined from various control sample, added linearly for the 6 tracks. • Branching fraction: This fraction of uncertainties come from uncertainties of photon decay to X + X −, X = e, µ, π. The uncertainty on the measurement of the ratio R [21] is propagated to the product of branching fraction. This uncertainty is mostly at the level of 1% for all category of events. All systematic uncertainties listed above are summed in quadrature to obtain the total systematic uncertainty for each category. We also combine the systematic uncertainties from each of the categories together. The correlations arising from the uncertainties on the PID and branching fractions are taken into consideration in the combination process. The combined total systematic uncertainties are shown in Fig. 6.20. The average scale of the systematic uncertainties is ∼ 11%; the main contributor is the PID uncertainty. There is a horizontal band around 0.7 GeV ≤ m A0 ≤ 0.85 GeV that corresponds to the region where the A0 → π + π − branching fraction is higher. This band results from the higher C2 category of branching fractions in the region. Without considering the correlation among the three categories, the relation between combined relative rel ) and each category’s relative systematic uncertainty (σ rel , i = uncertainty (σcomb i C0, C1, C2 ) is: relative 2 (σcomb ) = Õ (σirelative )2 · wi (BFi · i )2 . wi = Í ( BFi · i )2 i 46 Figure 6.20: The total systematic uncertainties as a function of mΥD and m A0 for combined systematic uncertainties. 6.8 Signal Significance Estimation This section discusses the procedure to estimate the significance of a potential signal. There are three steps to accomplish this: 1. Sideband Dataset: Construct a sideband dataset to estimate the background characteristics, since no events were selected in the optimization sample. 2. Distribution Structure: Use the sideband dataset to study the background distribution as a function of (mΥD , m A0 ) and establish the background estimation procedure. 3. Signal significance: Determine the distribution of the number of observed signal events under the null hypothesis. Sideband dataset To obtain the distribution of background events in a signal region, we would need to study background events selected by our classifiers, which is infeasible for two rea- 47 sons. First, there are no generators to simulate the relevant backgrounds accurately. Second, no events are selected in the optimization sample. Instead, we construct a sideband dataset to solve this challenge. Given our classifier with an optimal selection criteria optimal_cut on the classifier score, our sideband region is defined as [optimal_cut − δ, optimal_cut], where δ is a parameter we set. If δ is too small, we won’t have enough statistics, while the sideband dataset may be too biased if δ is too large. Our method is to select δ such that each sideband contains at least 500 events for each category of events, while making sure that δ is no smaller than the prediction uncertainty of each classifier, which is the standard deviation of optimal_cut estimated via bootstrapping. Run3 data events falling into the sideband regions constitute the sideband datasets. We claim that sideband datasets are a good approximation of the background distribution in the signal region. Fig. 6.21 illustrates the sideband region for the C2 category of events. Figure 6.21: Distribution of the ml_score for the C2 category for data. The sideband region is shown as a red rectangle. The optimal cut on the classifier score is shown as a solid line. The width of signal region is the δ we set to contain at least 500 events. Background structure To understand the distribution of background events in the signal region, we plot the distributions of sideband events on mΥD − m A0 space, shown in Fig. 6.22. The dimuon and dipion thresholds are indicated by horizontal lines. For the C0 and C1 category of events, the distribution of sideband events are uniform. There is also a small band for photon conversions in the C0 sample. For the C2 category of events, 48 the distribution of sideband events is also uniform, but there is a one-dimensional band when the dark photon mass is in the range of 0.7-0.8 GeV, corresponding to the ρ meson. The distributions of sideband events provide a prescription on how to estimate the number of background events in the scanning windows. Under the assumption that the background distributions in the signal region are smooth, a nearest-neighbor method can be applied to estimate the number of background events. Specifically, the estimated number of background events B within a scanning window is obtained by calculating the average number of observed events of its nearest left and right scanning windows, as illustrated in Fig. 6.23. We do not take up and down windows due to the horizontal band structure in the background distribution. If the scanning window is at the left (right) boundary of mΥD − m A0 space, we only use the nearest right (left) scanning window for estimation. When the scanning window is above dimuon threshold, the uniformity of the background guarantees the unbiasedness of the estimated number of background events. Using the left and right windows also works when we scan over the one-dimensional band produced by ρ or ω meson decays. When the scanning window is under the dimuon threshold, it is still a good estimation as long as background distribution is smooth along the ΥD axis. Signal significance Generally, three factors are needed to estimate the signal significance: (1) Hypotheses: our null hypothesis (H0 ) is that there is no darkonium signal. (2) Statistic: a test statistic that has discrimination between the null hypothesis (no signal) and alternative hypothesis (have signal). The statistic we use is the maximum number max among all signal windows. For each signal window, of observed signal events Nsig the number of observed signal events equals the number of observed events minus the estimated background: Nsig = Nobs − B. (3) Distribution of statistic: we use max . The distribution a bootstrapping technique to estimate the distribution of Nsig max will be used to calculate the signal significance (p-value) of the observed of Nsig signal events in the final data. max . First, we estimate the expected There are two steps to obtain the distribution of Nsig number of background events based on the optimal expected number of observed max distribution with a bootstrap events on data. The second step estimates the Nsig procedure (i.e, Toy MC). For each category Ci , we sample a number of events Ni from a Poisson distribution. We then bootstrap Ni events from the sideband dataset, A′ Mass (GeV) 49 C0 final states background distribution 3.0 2.5 Di-Pion threshold Di-Muon threshold 2.0 1.5 1.0 0.5 A′ Mass (GeV) 0.0 0 2 6 8 10 8 10 ϒD mass (GeV) C1 final states background distribution 3.0 2.5 4 Di-Pion threshold Di-Muon threshold 2.0 1.5 1.0 0.5 A′ Mass (GeV) 0.0 0 3.0 2.5 2 4 6 ϒD mass (GeV) C2 final states background distribution Di-Pion threshold Di-Muon threshold 2.0 1.5 1.0 0.5 0.0 0 2 4 6 8 10 ϒD mass (GeV) Figure 6.22: The distributions of sideband events for the (upper) C0 , (middle) C1 , and (bottom) C2 categories of events. A′ Mass (GeV) 50 3.0 Illustration of the background estimation technique 2.5 2.0 1.5 1.0 0.5 0.0 0 2 4 6 8 10 ϒD mass (GeV) Figure 6.23: Illustration of the background estimate method for a given mass hypothesis in the mΥD − m A0 space. The dataset displayed here is sideband data for C0 category of events. and run our scanning procedure on the mΥD − m A0 space. Since all categories are summed together when establishing scanning windows, a single scan is necessary here. The maximum number of events among all signal windows gives one entry of max distribution. We run 10,000 bootstraps to obtain the distribution of N max the Nsig sig shown in Fig. 6.24. 51 Frequency Distribution of maximum number of events 100 43.260 % 56.650 % 10−1 10−2 0.090 % 10−3 10−4 0 1 2 Maximum number of events Figure 6.24: The distribution of the maximum number of signal events after scanning the parameter space 0 < mΥD < 10 GeV, 0 < m A0 < 3 GeV. 6.9 Signal Extraction and Upper Limit of Coupling In this section, we first describe the procedure to extract the potential signal, then the results of upper limit of the coupling constant we obtain from data. Signal extraction We first combine the three categories (C0, C1, C2 ) together before scanning for the signals: (1) The combined ΥD and dark photon mass resolutions are taken as the maximum mass resolutions of three categories: σmΥD = max{σmΥD (C0 ), σmΥD (C1 ), σmΥD (C2 )} σm A0 = max{σm A0 (C0 ), σm A0 (C1 ), σm A0 (C2 )}. (2) The signal efficiencies are combined as the average of a categories’ signal efficiencies weighted by their branching fractions ¯ = BF(C0 ) · C0 + BF(C1 ) · C1 + BF(C2 ) · C2 . (3) The systematic uncertainties from three categories are combined. The signal is extracted by scanning the ΥD mass versus the dark photon mass plane. As mentioned above, the signal region is defined as a rectangular window mΥD − 4σmΥD < mΥD < mΥD + 4σmΥD and m A0 − 4σm A0 < m A0 < m A0 + 4σm A0 for a 52 given ΥD and photon mass (mΥD , m A0 ). An estimate of the background is obtained from the nearest left and right signal windows. The signal significance is obtained from the distribution of the maximum number of signal events derived from the bootstrap procedure (Fig. 6.24). We do not see a clear signal; therefore upper limits on the signal cross section are extracted, as described below. Suppose for a certain scanning window, we observe ni events passing the selection criteria for each category Ci (i = 0, 1, 2). The combined signal efficiency is . ¯ The estimated number of background events is B̂. The combined systematic uncertainty Í is σtot . We assume the total number of observed events n = i=0,1,2 ni follows a Poisson distribution with parameter λ = Lσ + B, where L is luminosity, σ is cross section of dark Upsilon decaying to three dark photons,  is unknown true signal efficiency, and B is unknown true number of background. We also assume that the evaluated signal efficiency ¯ follows a Gaussian distribution whose mean value is the true signal efficiency  and the standard deviation is the estimated total 2 . The likelihood for this scanning window is: systematic uncertainty σtot 2 L(n|λ = Lσ + B) = Poisson(n|λ) · N (¯ |, σtot ) · Poisson( B̂|B). There are two nuisance parameters in the likelihood,  and B. We apply a profile likelihood method to obtain a 90% confidence level limit on the cross section. For a likelihood L with parameters (θ, ψ), where θ are the nuisance parameters and ψ are the parameters of interest, the profile likelihood method takes two steps to optimize the L: For each ψ, we have a curve Lψ (λ) over λ. We first evaluate the maximal L over λ, then we choose the ψ that is maximum over all these curves. To take the signal lifetime into consideration, we apply an iterative algorithm to find the converged signal efficiency and converged upper limits. We use Equation 6.7 to compute the theoretical dark photon lifetime and Equation 6.10 to correct the signal efficiencies. The detailed algorithm is described in Table 6.3. Results We do not observe a significant signal of dark matter bound state ΥD . The 90% CL results on cross sections for both prompt and displaced signal decays are shown in Fig. 6.25 and 6.26. We derive separate limits for αD values set to 0.1, 0.3, 0.5, 0.7, 0.9, and 1.1. Constraints for different values of αD , m A0 , and mΥD are shown in Fig. 6.27. The results compared with previous search are shown in Fig. 6.28 for different αD . Bounds on the mixing strength ε down to 5 × 10−5 − 10−3 are set for 53 Algorithm to derive the limit when the lifetime is non-zero Step 0: Compute the zero-lifetime kinetic mixing strength upper limit () using profile likelihood method. Step 1: Calculate the corresponding dark photon lifetime (l) using Equation 6.7. If the lifetime is smaller than 10 um, go to Step 5; otherwise go to Step 2. Step 2: Update the signal efficiency under the dark photon lifetime l using Equation 6.10. Step 3: Update the corresponding mixing strength upper limit  using the updated signal efficiency. Step 4: Repeat Step 1 until the procedure has converged. Step 5: Return the converged signal efficiency, the cross-section upper limit and the mixing strength upper limit. Table 6.3: Algorithm to find converged signal efficiency, cross-section upper limits, and kinetic mixing strength upper limits. a wide range of dark photon mass from MeV to few GeV. Compared with previous search, our analysis shows that the constraints on γ − A0 coupling strength can be improved by a large fraction of parameter space and can be more than one order of magnitude for 40 MeV ≤ m A0 ≤ 200 MeV with strong electromagnetic interaction in the dark sector, demonstrating the dark matter bound state as a sensitive probe of dark matter. Figure 6.25: The 90% CL upper limits on the e+ e− → γΥD cross section for prompt dark photon decays. 54 Figure 6.26: The 90% CL upper limits on the e+ e− → γΥD cross section for dark photon lifetimes corresponding to (top left) cτA0 = 0.1 mm, (top right) cτA0 = 1 mm, and (bottom) cτA0 = 10 mm. 6.10 Summary and Outlook In summary, we explore the sensitivity of collider probes on dark sector searches with the existence of dark sector structures. Dark matter bound states could exist under a simple dark photon model when the dark photon is light enough to generate an attractive force between dark fermions. We report the first search for a dark sector bound state decaying into three dark photons in the range 0.001 GeV < m A0 < 3.16 GeV and 0.05 GeV < mΥD < 9.5 GeV. We do not observe significant signals. The limits on the γ − A0 kinetic mixing  are derived at the level of 5 × 10−5 − 10−3 , depending on the values of the model parameters. These measurements improve upon existing constraints over a significant fraction of dark photon masses below 1 GeV for large values of the dark sector coupling constant. In the future, the search for ηD bound state can also be also included. With the current experiment data, the upper limits on the cross section (in the absence of a signal) could be improved by around a factor of 2, leading to an improvement on √ the constraints on the kinetic mixing strength by about a factor of 2. With 50 ab−1 luminosity collected by future B factory experiments (Belle II), the limits on the 55 Figure 6.27: The 90% CL upper limits on the kinetic mixing  2 as a function of the ΥD mass, mΥD , and dark photon mass, m A0 , assuming (top left) αD = 0.1, (top right) αD = 0.3, (middle left) αD = 0.5, (middle right) αD = 0.7, (bottom left) αD = 0.9, and (bottom right) αD = 1.1. kinematic mixing strength can potentially reduce by at least one order of magnitude. 56 ε 10−2 10−3 mΥD = 3 GeV mΥD = 6 GeV mΥD = 9 GeV 10−4 −3 10 10−2 10−1 100 101 102 mA0 (GeV) ε 10−2 10−3 αD = 0.3 αD = 0.5 αD = 0.7 αD = 0.9 10−4 10−5 −3 10 BaBar αD = 1.1 10−2 10−1 100 101 102 mA0 (GeV) Figure 6.28: The 90% CL upper limits on the kinetic mixing ε for (top) various ΥD masses assuming αD = 0.5 and (bottom) various αD values assuming mΥD = 9 GeV together with current constraints (gray area). Part IV Precise Measurement of Semileptonic B Meson Decays 57 58 Chapter 7 MEASUREMENT OF R(D) AND R(D∗ ) USING SEMILEPTONIC TAGGING AND LEPTONIC τ DECAYS 7.1 Introduction The excess of semitauonic decays B → D(∗) τν have been one of the most interesting puzzles in flavor physics in recent years [22]. The physical quantity to measure is the ratio between B → D(∗) τν and B → D(∗) lν: B(B → D(∗) τ − ντ ) R(D ) = (l = e, µ) B(B → D(∗) l − νl ) (∗) Figure 7.1: Feynman diagram of B → D(∗) τν decay. In the Standard Model, B → D(∗) τν decay involves the semileptonic quark transition b → cl ν̄ with l = e, µ, τ. The Feynman diagram is shown in Fig. 7.1. The b quark in the B meson decays to c quark via first-order electroweak interactions, and the mediator, a W boson, decays into pair of leptons l ν̄. The other quark q̄ in the B meson binds with the c quark in the final state to form a D(∗) meson. The Lagrangian of the system is GF +µ L = √ |Vcb |Jνl Jcb,µ + h.c., 2 µ which involves the quark currents Jcb = ψ̄c γ µ (1 − γ 5 )ψb and the leptonic currents µ Jνl = ψ̄ν γ µ (1 − γ 5 )ψl . The amplitudes of this scattering process are: Õ GF M λλlM (q2, x) = √ Vcb ηλW LλλWl HλλWM , 2 λW where λ is the particle helicities, M = D, D∗ indicates the final state meson. The leptonic amplitudes LλλWl can be analytically calculated with the standard framework of electroweak interactions [23]. Although the analytical calculation of hadronic 59 amplitudes is intractable because of complex interactions between quarks and the self-interactions of the gluons, they can be expressed in terms of six functions, or the form factors depending only on q2 : f± (q2 ) and fi (q2 ), i = 1, 2, 3, 4, where q = pB − pD(∗) is the four-momentum of the virtual W. One model to form factors is the Caprini-Lellouch-Neubert (CLN) model [24], which introduces dispersive constraints from heavy quark symmetry to provide relations between the form factors near zero recoil. The SM predictions for R(D(∗) ) are: R(D)SM = 0.299 ± 0.003, R(D∗ )SM = 0.254 ± 0.005. The uncertainties on R(D(∗) )SM come from the uncertainties of the form factor parameters. If new particles couple proportionally with their masses, B → D(∗) τν are sensitive probes for non-SM contributions due to the heavy τ mass. The individual decay modes will suffer from large hadronic uncertainties related to the form factors and the Cabibbo-Kobayashi-Maskawa (CKM) elements. Normalizing the B → D(∗) τν to the corresponding decay with light leptons provides better sensitivity to new physics (NP). Since the energy scale of NP should be far above the scale of B meson, we can integrate out higher degrees of freedom and form the following effective Hamiltonian: cb cb Heff = CSM OSM + CRcbO Rcb + CLcbO Lcb, where cb OSM = q̄γ µ PL bτ̄γ µ PL ντ O Rcb = q̄PR bτ̄PL ντ O Lcb = q̄PL bτ̄PL ντ . The first term is the SM charged current. The second and third terms are the fourfermion operators allowed by Lorentz invariance [25]. The SM Wilson coefficient cb = 4G√ F Vcb CSM . The corresponding Wilson coefficients CRcb and CLcb affect the 2 observables in the following way: " © R(D) = RSM (D) ­1 + 1.5R « CRcb + CLcb cb CSM # C cb + C cb ª + 1.0 R cb L ® CSM ¬ 2 60 C cb − C cb C cb − C cb ª © R(D ) = RSM (D ) ­1 + 0.12R R cb L − 0.05 R cb L ® . CSM CSM « ¬ " ∗ # 2 ∗ Therefore, any deviation between the SM prediction and observed value of R(D(∗) ) may imply existence of NP, including • Supersymmetry: In the Minimal Supersymmetric Standard Model (MSSM), a charged Higgs boson which couples proportionally to the masses of the fermions is introduced. The Higgs only couples significantly to τ, so that it contributes to the deviation between SM prediction and observed value of measured quantity. • Lepton flavor universality (LU) violation: The b → clν transition is a FlavorChanging Charged Current (FCCC) test. If LU is violated, the Wilson coefficients would be different for the operators with the same structure but different lepton flavors, contributing additional term to R(D(∗) ). The first measurement of R(D(∗) ) were at BABAR in 2008 [26], and the first observation of an excess was on 2012 [27]. The measurements use hadronic tags to fully reconstruct B mesons and leptonic τ decays. The measured value has a 3.4σ deviation from SM predictions. Since then, many other experiments have contributed their measurements to the quantity of interest [28, 29, 30]. As of 2019, the R(D) and R(D∗ ) exceeds the SM prediction by 1.4σ and 2.5σ respectively. Considering the R(D) − R(D∗ )) correlation of -0.38, The difference with the SM predictions reported above, corresponds to about 3.08σ. Fig. 7.2 and Table 7.1 show a summary of previous measurements of this quantity. The goal of this BABAR analysis is to provide another measurement of this quantity using, for the first time, a semileptonic tagging method. 61 Figure 7.2: Visualization of previous measurements of R(D) and R(D∗ ). The plot is from HFLAV [22]. BABAR 2013 [23] Belle 2015 [28] LHCb 2015 [29] Belle 2016 [31] Belle 2017 [32] LHCb 2018 [33] Belle 2019 [30] Average (HFLAV Summer 2018) Standard Model [22] R(D) 0.440 ± 0.058 ± 0.042 0.375 ± 0.064 ± 0.026 0.307 ± 0.037 ± 0.016 0.340 ± 0.027 ± 0.013 0.299 ± 0.003 R(D∗ ) 0.332 ± 0.024 ± 0.018 0.293 ± 0.038 ± 0.015 0.336 ± 0.027 ± 0.03 0.302 ± 0.030 ± 0.011 0.270 ± 0.035+0.028 −0.025 0.291 ± 0.019 ± 0.029 0.283 ± 0.018 ± 0.014 0.295 ± 0.011 ± 0.008 0.258 ± 0.005 Table 7.1: Previous results. 7.2 Analysis Strategy Overview Event types definition For the B B̄ system, the signal events are defined as Btag → D(∗) lν, Bsig → D(∗) τν. On the tag B side, the B meson decays to D(∗) lν, where l = e, µ. On the signal B side, the B meson decays to D(∗) τν, and the τ subsequently decays to leptons. We have two types of signal events, based on whether there is a D or D∗ meson. The normalization events are defined as events with both the tag side B meson and 62 signal side B meson decaying to D(∗) lν. One of the most important category of backgrounds arises from a B → D∗∗ lν events where D∗∗ denotes the excited charm meson states heavier than D∗ because of the similar decay topology to signal events. Two other sources of background are combinatorial B B̄ events and continuum events. The definition of these events is listed in Table 7.2. Event type Signal event Normalization event signal D signal D∗ norm D norm D∗ D∗∗ event Combinatorial B B̄ event Continuum event Description One B decays to D(∗) lν, the other B decays to Dτν, τ → leptons One B decays to D(∗) lν, the other B decays to D∗ τν, τ → leptons One B decays to D(∗) lν, the other B decays to Dlν Both B decay to D∗ lν At least one B decays to D∗∗ (l/τ)ν, where D∗∗ denotes any excited charmed meson states that are not in the ground state 1S doublet. In this analysis, it includes 1P states D0∗, D1, D10 , D2∗ , and 2S states. We also allow non-resonant final states consisting of a D(∗) and one pion. Any B B̄ events that are not signal and not normalization and not D∗∗ . Non-B B̄ events produced in the detector Table 7.2: Definition of event types in the B B̄ system. Reconstruction The first step of this analysis is to identify primitive particles (e.g, e, µ, π ±, K, γ), and reconstruct composite particles (e.g, π 0 , KS , D(∗) , B) by applying a series of selection criteria. When we reconstruct signal events, we first reconstruct a Btag , then search for D(∗) l in the remaining tracks and calorimeter clusters. The signal events can be categorized into four disjoint subsets based on the reconstructed type of D meson: D+ l, D0 l, D∗+ l, D∗0 l. We denote these as channel_labels for convenience. The details of the reconstruction and selection criteria are described in Section 7.4. Estimate of R(D(∗) ) To measure R(D(∗) ), let us denote P := B(B → Dτν) P∗ := B(B → D∗ τν) Q := B(B → Dlν)(average for l = e or µ) (7.1) Q∗ := B(B → D∗ lν)(average for l = e or µ). Therefore, R(D) = QP and R(D∗ ) = QP ∗ . If we denote the number of generated B B̄ event as N, the relationship between the B decay branching fractions and the ∗ 63 expected number of signal events generated in the detector would be N(signal D) = 2N · (2Q + 2Q∗ ) · P · B(τ → leptons) N(signal D∗ ) = 2N · (2Q + 2Q∗ ) · P∗ · B(τ → leptons) N(norm D) = 4N · (Q2 + 2QQ∗ ) (7.2) N(norm D∗ ) = 4N · Q∗ 2 . Given the estimated number of generated signal events N̂(signal D) and N̂(signal D∗ ) and normalization events N̂(norm D) and N̂(norm D∗ ), the estimated P(∗) can be derived from Equation (7.2): N̂(signal D) P̂ = √ p 2 N · B(τ → leptons) · N̂(norm D) + N̂(norm D∗ ) N̂(signal D∗ ) Pˆ∗ = √ p 2 N · B(τ → leptons) · N̂(norm D) + N̂(norm D∗ ) p p N̂(norm D) + N̂(norm D∗ ) − N̂(norm D∗ ) Q̂ = √ 4N s N̂(norm D∗ ) Q̂∗ = . 4N (7.3) Therefore, the estimated R(D(∗) ), as functions of the estimated number of generated events N̂(signal D) and N̂(signal D∗ ), would be R(D) = 1 · 2B(τ → leptons) p R(D∗ ) = N̂(signal D) p  p ∗ ∗ ∗ N̂(norm D) + N̂(norm D ) · N̂(norm D) + N̂(norm D ) − N̂(norm D ) 1 · 2B(τ → leptons) N̂(signal D∗ ) . p p N̂(norm D) + N̂(norm D∗ ) · N̂(norm D∗ ) (7.4) The branching fraction of τ decays to leptons is taken to be [34] B(τ → leptons) = B(τ → e ν̄e ντ ) + B(τ → µν̄ µ ντ ) = (35.24 ± 0.05)% (7.5) N̂(signal D), N̂(signal D∗ ), N̂(norm D), and N̂(norm D∗ ) can be estimated as the extracted yields divided by the corresponding efficiencies. 64 7.3 Simulation Samples The simulated data includes generic MC, signal MC, as well as normalization MC events. The generic MC aims to represent the data collected in the detector by faithfully generating the possible results of the e+ e− collisions. The generic MC includes two types of events: B B̄ events and continuum events. B B̄ events emulate both charged and neutral B B̄ meson pairs produced from Υ(4S). Continuum sample simulates e+ e− → q q̄. Other constributions (QED, two-photon processes) are negligible. The weighted collection of these types of events gives our best reproduction of what an on-peak collider run produces. To weight the generic MC, we need to consider the cross section of each type of event by the following formula: σi , Ni where L is the integrated detector data luminosity, σi is the corresponding cross section for event component i, and Ni is the number of event components i generated. The cross sections for each event type are listed in Table 7.3. The components of the generic MC dataset are listed in Table 7.4. wi = L SP Mode 1235 1237 1005 998 Mode type B+B− B0B0 cc uds Cross section (pb) 525.0 525.0 1300.0 2090.0 Table 7.3: Cross sections used to convert the sizes generic simulated data to the equivalent on-peak dataset. The signal (normalization) MC is generated by forcing every B meson decay in the signal (normalization) modes. The datasets we use are listed in Table 7.5 and Table 7.6, in which both charged and neutral B B̄ meson pairs are simulated. The signal (normalization) MC datasets are used in the analysis described later. 65 Simulated Dataset Name Mode Type SP-1235-AllEventsSkim-Run1-R24a1 B+B− SP-1235-AllEventsSkim-Run2-R24a1 B+B− SP-1235-AllEventsSkim-Run3-R24a1 B+B− SP-1235-AllEventsSkim-Run4-R24a1 B+B− SP-1235-AllEventsSkim-Run5-R24a1 B+B− SP-1235-AllEventsSkim-Run6-R24a1 B+B− SP-1237-AllEventsSkim-Run1-R24a1 B0B0 SP-1237-AllEventsSkim-Run2-R24a1 B0B0 SP-1237-AllEventsSkim-Run3-R24a1 B0B0 SP-1237-AllEventsSkim-Run4-R24a1 B0B0 SP-1237-AllEventsSkim-Run5-R24a1 B0B0 SP-1237-AllEventsSkim-Run6-R24a1 B0B0 SP-1005-AllEventsSkim-Run1-R24a1 cc SP-1005-AllEventsSkim-Run2-R24a1 cc SP-1005-AllEventsSkim-Run3-R24a1 cc SP-1005-AllEventsSkim-Run4-R24a1 cc SP-1005-AllEventsSkim-Run5-R24a1 cc SP-1005-AllEventsSkim-Run6-R24a1 cc SP-998-AllEventsSkim-Run1-R24a1 uds SP-998-AllEventsSkim-Run2-R24a1 uds SP-998-AllEventsSkim-Run3-R24a1 uds SP-998-AllEventsSkim-Run4-R24a1 uds SP-998-AllEventsSkim-Run5-R24a1 uds SP-998-AllEventsSkim-Run6-R24a1 uds Collisions Generated 34878000 105561000 56035000 166784000 215168000 130336000 34941000 104188000 57888000 169801000 215953000 135224000 55254000 164722000 88321000 267308000 344275000 208664000 176404000 525504000 276381000 845899000 1110944000 655152000 Multiplier 0.306 0.305 0.303 0.314 0.323 0.316 0.306 0.308 0.292 0.307 0.321 0.304 0.479 0.483 0.475 0.484 0.499 0.488 0.241 0.243 0.244 0.246 0.249 0.250 Table 7.4: Generic simulated data. Multiplier is the factor by which the size of the corresponding on-peak dataset exceeds that of the given simulated dataset. More simulated datasets are usually generated to better study the decay characteristics. 7.4 Event Reconstruction The event reconstruction procedure includes three steps, all of which are implemented using the BABAR offline analysis framework (version-52). 1. Pre-screening: We use a collection of loose criteria to remove obvious background events. This step helps reduce the amount of data to analyze. 2. Candidate reconstruction: We then attempt to reconstruct the Btag and Bsig of signal events from a set of tracks of final state particles. The event reconstruction is a hierarchical process: we first identify primitive particles (e.g, 66 Simulated Dataset Name Mode Type SP-11440-Run1-R24 SP-11440-Run2-R24 SP-11440-Run3-R24 SP-11440-Run4-R24 SP-11440-Run5-R24 SP-11440-Run6-R24 SP-11441-Run1-R24 SP-11441-Run2-R24 SP-11441-Run3-R24 SP-11441-Run4-R24 SP-11441-Run5-R24 SP-11441-Run6-R24 SP-11442-Run1-R24 SP-11442-Run2-R24 SP-11442-Run3-R24 SP-11442-Run4-R24 SP-11442-Run5-R24 SP-11442-Run6-R24 SP-11443-Run1-R24 SP-11443-Run2-R24 SP-11443-Run3-R24 SP-11443-Run4-R24 SP-11443-Run5-R24 SP-11443-Run6-R24 B0 → D(∗) `ν, B0 → Dτ(e, µ)ν B0 → D(∗) `ν, B0 → Dτ(e, µ)ν B0 → D(∗) `ν, B0 → Dτ(e, µ)ν B0 → D(∗) `ν, B0 → Dτ(e, µ)ν B0 → D(∗) `ν, B0 → Dτ(e, µ)ν B0 → D(∗) `ν, B0 → Dτ(e, µ)ν B0 → D(∗) `ν, B0 → D∗ τ(e, µ)ν B0 → D(∗) `ν, B0 → D∗ τ(e, µ)ν B0 → D(∗) `ν, B0 → D∗ τ(e, µ)ν B0 → D(∗) `ν, B0 → D∗ τ(e, µ)ν B0 → D(∗) `ν, B0 → D∗ τ(e, µ)ν B0 → D(∗) `ν, B0 → D∗ τ(e, µ)ν B+ → D(∗) `ν, B− → Dτ(e, µ)ν B+ → D(∗) `ν, B− → Dτ(e, µ)ν B+ → D(∗) `ν, B− → Dτ(e, µ)ν B+ → D(∗) `ν, B− → Dτ(e, µ)ν B+ → D(∗) `ν, B− → Dτ(e, µ)ν B+ → D(∗) `ν, B− → Dτ(e, µ)ν B+ → D(∗) `ν, B− → D∗ τ(e, µ)ν B+ → D(∗) `ν, B− → D∗ τ(e, µ)ν B+ → D(∗) `ν, B− → D∗ τ(e, µ)ν B+ → D(∗) `ν, B− → D∗ τ(e, µ)ν B+ → D(∗) `ν, B− → D∗ τ(e, µ)ν B+ → D(∗) `ν, B− → D∗ τ(e, µ)ν Nevent Collected 2621000 2756000 2789000 2356000 Table 7.5: Simulated signal data. e, µ, γ), then reconstruct D(∗) and B mesons, and then combine B and B̄ pairs to reconstruct Υ(4S) candidates. 3. Candidate selection: When multiple candidates are involved in an event after reconstruction, we select a single best candidate, as defined below. Pre-screening The following set of broad selection criteria is first applied to loosely pre-select signal events. These selection criteria, studied from MC samples, are set to exclude obvious backgrounds but keep most signal and normalization events. • Size of ChargedTracks ≤ 14, 67 Simulated Dataset Name SP-11438-Run1-R24 SP-11438-Run2-R24 SP-11438-Run3-R24 SP-11438-Run4-R24 SP-11438-Run5-R24 SP-11438-Run6-R24 SP-11439-Run1-R24 SP-11439-Run2-R24 SP-11439-Run3-R24 SP-11439-Run4-R24 SP-11439-Run5-R24 SP-11439-Run6-R24 Mode Type Nevent Collected B+ → D(∗) `ν, B− → D(∗) `ν B+ → D(∗) `ν, B− → D(∗) `ν B+ → D(∗) `ν, B− → D(∗) `ν B+ → D(∗) `ν, B− → D(∗) `ν B+ → D(∗) `ν, B− → D(∗) `ν B+ → D(∗) `ν, B− → D(∗) `ν B0 → D(∗) `ν, B0 → D(∗) `ν B0 → D(∗) `ν, B0 → D(∗) `ν B0 → D(∗) `ν, B0 → D(∗) `ν B0 → D(∗) `ν, B0 → D(∗) `ν B0 → D(∗) `ν, B0 → D(∗) `ν B0 → D(∗) `ν, B0 → D(∗) `ν 2443000 2869000 Table 7.6: Simulated normalization data. • Size of GoodPhotonsLoose ≤ 10, • −2 ≤ Qtotal ≤ 2, • Apply tag filter BGFMultiHadron, • Apply tag filter TagL3. Candidate reconstruction Events passing the pre-screening are then used to reconstruct Btag and Bsig . Initially, each event is a collection of tracks of final-state particles. To reconstruct the signal candidates, we first identify final-state particles using built-in Particle Identification algorithms. The loosest selectors are applied to identify e, µ, π ±, K, and γ. The π 0 mesons are then reconstructed from π 0 → γγ by combining a pair of photons, and KS is reconstructed from the KS → π + π − mode. The following are applied to select light mesons: • π 0 : pi0AllDefault and pi0SoftDefaultMass. Using the photon list above, reconstruct π 0 → γγ with the mass of the pion candidate constrained to be between [0.115, 0.15] GeV. • KS : KsDefault. Reconstruct KS → π + π − using TreeFitter. 68 The D mesons are reconstructed based on the identified final state particles and light mesons. The D decay modes listed in Table 7.7 are considered in reconstruction, which sum up to 19.2% and 30.1% of D+ and D0 branching fractions. Other decay modes are not considered in this analysis due to higher backgrounds. Decay Modes D+ → K − π + π + D + → KS π + D + → KS π + π 0 D+ → K − K + π + D0 → K − π + D0 → K − π + π 0 D0 → K − π + π − π + D 0 → KS π + π − D 0 → KS π 0 D0 → K − K + Branching Fraction (%) 9.46 ± 0.24 1.53 ± 0.06 7.24 ± 0.17 0.99 ± 0.026 3.88 ± 0.05 14.3 ± 0.80 8.06 ± 0.23 2.85 ± 0.20 1.20 ± 0.04 0.40 ± 0.007 Number in Generic MC 142202 22691 430750 62859 158718 1235214 646768 641954 94151 49909 Table 7.7: D meson decay modes used in the analysis. The mass of D mesons (m(D)) is used to select correctly reconstructed candidates. Since the mass resolutions vary for different types of D meson decay modes, our criteria on m(D) are different depending on the D decay modes. The distribution of m(D) and the criteria is shown in Fig. 7.3. – D decay without a π 0 in the final state: The mass resolution (σmD ) for this type of D decay modes is around 5 MeV. Therefore, we require m(D) within [-15, 15] MeV of the nominal mass, which corresponds to ±3σmD . – D+ decay with a π 0 in the final state: The mass resolutions for this type of D decay modes is around 12 MeV. We require m(D± ) within [-36, 24] MeV of the nominal charged D mass, corresponding to [-3σmD , +2σmD ]. – D0 decay with a π 0 in the final state: The mass resolutions for this type of D decay modes is around 15 MeV. We require m(D± ) within [-45, 30] MeV of the nominal neutral D mass, corresponding to [-3σmD , +2σmD ]. Density 50 B̄ → D(∗)l −ν̄l Continuum B̄ → D(∗)τ −ν̄τ BB̄ bkgs B̄ → D∗∗(τ /l)−ν̄ Criteria 50 40 30 30 20 20 10 10 Density 20.0 1.83 1.84 1.86 1.87 1.88 BB̄ bkgs B̄ → D∗∗(τ /l)−ν̄ 0 1.82 1.89 msig D (GeV) B̄ → D(∗)l − ν̄l Continuum 17.5 1.85 20.0 B̄ → D(∗)τ −ν̄τ Criteria 1.83 Criteria 1.84 1.85 1.86 1.87 1.88 1.89 msig D (GeV) B̄ → D(∗)l − ν̄l Continuum BB̄ bkgs 17.5 15.0 B̄ → D(∗)τ −ν̄τ B̄ → D∗∗(τ /l)−ν̄ 40 0 1.82 B̄ → D(∗)l −ν̄l Continuum BB̄ bkgs Density Density 69 B̄ → D∗∗(τ /l)−ν̄ B̄ → D(∗)τ −ν̄τ Criteria 15.0 12.5 12.5 10.0 10.0 7.5 7.5 5.0 5.0 2.5 2.5 0.0 1.80 1.82 1.84 1.86 1.88 1.90 msig D (GeV) 0.0 1.80 1.82 1.84 1.86 1.88 1.90 msig D (GeV) Figure 7.3: The distributions of reconstructed D meson mass for (upper left) D+ decays without π 0 , (upper right) D0 decays without π 0 , (bottom left) D+ decays with π 0 , and (bottom right) D0 decays with π 0 . The D∗ mesons are reconstructed in three decay modes. Other decay modes are not considered, due to the high backgrounds. The mass difference between reconstructed D∗ meson and its daughter D meson (∆M) is used to select well-reconstructed D∗ mesons. The resolution on ∆M is around 2 MeV. Tighter criteria are applied when there is a π 0 in the final state, as large backgrounds arise from mis-reconstructed neutral pions. The distributions of mass differences for each type of decays are shown in Fig. 7.4. • D∗+ → D0 π + : We require ∆M to be within 2.5 MeV of the nominal mass difference. • D∗+ → D+ π 0 : We require ∆M to be within 2.0 MeV of the nominal mass difference. • D∗0 → D0 π 0 : We require ∆M to be within 2.0 MeV of the nominal mass difference. BB̄ bkgs B̄ → D∗∗(τ /l)−ν̄ B̄ → D(∗)τ −ν̄τ Criteria 400 Continuum BB̄ bkgs 200 B̄ → D∗∗(τ /l)−ν̄ B̄ → D(∗)l −ν̄l B̄ → D(∗)τ −ν̄τ Criteria Density B̄ → D(∗)l −ν̄l Continuum 500 Density Density 70 160 B̄ → D(∗)l −ν̄l Continuum BB̄ bkgs B̄ → D∗∗(τ /l)−ν̄ 140 B̄ → D(∗)τ −ν̄τ Criteria 120 150 100 300 80 100 60 200 40 50 100 20 0 0.142 0.143 0.144 0.145 0.146 0.147 0.148 0.149 ∆msig (GeV) 0 0.136 0.138 0.140 0.142 0.144 0.146 0.148 ∆msig (GeV) 0 0.136 0.138 0.140 0.142 0.144 0.146 0.148 0.150 ∆msig (GeV) Figure 7.4: The distribution of mass difference between reconstructed D∗ meson and its daughter D meson for (left) D∗+ → D0 π + , (middle) D∗+ → D+ π 0 , and (right) D∗0 → D0 π 0 decays. The Btag mesons are reconstructed by combining a D(∗) and an electron or muon candidate. The D(∗) l invariant mass is required to be at most 5.2791 GeV, and the p-value for the χ2 test of the kinematic fit must be at least 0.001. The variable used to select a correctly reconstructed Btag is the angle between the 3-momentum of tag the Btag and the 3-momentum sum of its D(∗) and lepton daughters (cos θ B−D(∗) l ), defined as tag cos θ B−D(∗) l = tag 2Ebeam ED(∗) l − m2B − m2D(∗) l 2|pB | · |pD(∗) l | . The value of cos θ B−D(∗) l should be in the range of [-1, 1], as it has only one missing neutrino. Other events (signal, combinatorial, continuum events) tend to have more tag tag negative cos θ B−D(∗) l values. We select events requiring cos θ B−D(∗) l ∈ [−2, 1] to take the detector resolutions into consideration. • B + → D 0 e+ . • B+ → D0 µ+ . • B 0 → D − e+ . • B0 → D− µ+ . • B+ → D∗0 e+ . • B+ → D∗0 µ+ . • B0 → D∗− e+ . • B0 → D∗− µ+ . 71 In each event with a selected Btag candidate, we search for D(∗) l among the remaining tracks and calorimeter clusters to reconstruct Bsig , using the B decay modes described above. The invariant mass of D(∗) l candidates is required to be at most 5.2791 GeV, and the p-value of the χ2 test must be at least 0.001. The Υ(4S) candidates are reconstructed by combining Btag and Bsig , requiring that the two B mesons must conserve charge. To further suppress backgrounds, we require: • No extra charged tracks, KS0 or π 0 particles must be present. • The extra neutral energy in the calorimeter Eextra ≤ 1.2 GeV. • The second Fox-Wolfram moment R2 ≤ 0.4. • The 3-momentum magnitude of the Bsig ’s lepton daughter in the CM frame sig |pl | ≤ 2 GeV. Candidate selection We use the following algorithm to assign a truth-matched candidate for each signal MC event. The criteria for the best candidate is clear: if the reconstruction graph of a candidate matches exactly that of the truth, that is our best candidate. The truth matching algorithm is described in detail in [35]. About 91% of the reconstructed events have only a single candidate. If multiple candidates are observed in an event, we choose the candidate with the minimal Eextra as the best candidate to represent the event. Following this method, 95.8% of the reconstructed signal candidates are the truth-matched candidate. Reconstruction results Table 7.8 shows the composition of events selected by the reconstruction strategy evaluated on generic MC. 1.56% (3.18%) of our reconstructed generic MC sample are B → Dτν (B → D∗ τν) events. The normalization events (B → D(∗) lν) sum up to around 38%. Combinatorial B B̄ background dominates and is roughly 28%, and we also have around 6.4% of continuum events. 72 Component B → Dτν B → D∗ τν B → Dlν B → D∗ lν B → D∗∗ lν Combinatorial B B̄ Continuum Proportion (%) 1.56 3.18 19.80 18.11 22.74 28.17 6.44 Table 7.8: Proportions of each component after reconstruction, evaluated using the generic MC. 7.5 Signal Detection We apply machine learning models to classify each type of event after reconstruction. Two binary classifiers C1, C2 are used, the corresponding scores denoted as z1 and z2 . C1 is used to identify signal events and normalization events from background D∗∗ , B B̄ combinatorial and continuum events. C2 is used to separate signal and normalization events. The combination of z1 and z2 enables us to measure signal events based on their density shape difference from other event types. We do not apply any selection criteria on z1 or z2 , however, the (z1, z2 ) spectrum of all event types is used to extract the signal yields. C1 classifier The sample used to train and validate the C1 classifier includes 350K events from the generic MC. The sample is divided into two parts using random sampling without replacement: training sample (240K) and validation sample (110K). We first use the training sample to train different classifiers and then use the validation sample to choose the best classifier. Generally, a binary classifier aims to classify two categories of events (often labeled as positive/negative), based on a set of input variables. C1 is a binary classifier to distinguish signal events and normalization events, from all types of background events. When training the C1 classifier, we label signal and normalization events as positive, and label D∗∗ , combinatorial B B̄, and continuum events as negative. The following variables are used as inputs to train the C1 classifier. The histograms of these variables for positive label (signal and normalization events) and each type of background event are shown in Fig. 7.5, 7.6, and 7.7. • Ntracks : Number of charged tracks. 73 • R2 All: Second Fox-Wolfram moment. 2 : Square of the missing 4-vector of the event in the CM frame. • Mmiss • Eextra : Extra neutral energy in the calorimeter. • cos θT : Cosine of the angle between the thrust and the beam momentum. tag • |pl |: 3-momentum magnitude of the Btag lepton in the CM frame. tag • cos θ B−D(∗) l : Cosine of the angle between the 3-momentum of the Btag and the 3-momentum sum of its D and lepton daughters in the CM frame. tag • cos θ D−l : Cosine of the angle between the 3-momentum of the D meson and the lepton daughter in the tag side. tag • mD : Mass of the Btag D meson daughter. • ∆mtag : Mass difference between D∗ and D meson in the tag side, if exists. tag • cos θ D so f t : Cosine of the angle between the D∗ mesons daughters in the tag side in the CM frame. tag • |pso f t |: 3-momentum magnitude of the D∗ ’s soft daughter in the tag side in the CM frame. sig • |pl |: 3-momentum magnitude of the Bsig lepton daughter in the CM frame. sig • cos θ D−l : Cosine of the angle between the 3-momentum of the D meson and the lepton daughter in the sig side in the CM frame. • χ2 : χ2 of the Bsig vertex fit. sig • mD : Mass of the Bsig D meson daughter. • ∆m sig : Mass difference between D∗ and D meson in the sig side, if it exists. sig • cos θ D so f t : Cosine of the angle between the D∗ mesons’ daughters in the sig side in the CM frame. sig • |pso f t |: 3-momentum magnitude of the D∗ ’s soft daughter in the sig side in the CM frame. • cos θ D(∗) l−D(∗) l : Cosine of the angle between the two Dl Dl systems in the CM frame. 74 • tag l electron PID: Btag lepton daughter’s electron PID level. • tag l muon PID: Btag lepton daughter’s muon PID level. • sig l electron PID: Bsig lepton daughter’s electron PID level. • sig l muon PID: Bsig lepton daughter’s muon PID level. Figure 7.5: Histograms of variables used for the C1 classifier. 75 Figure 7.6: Histograms of variables used for the C1 classifier. 76 Figure 7.7: Histograms of variables used for the C1 classifier. We use a Gradient Boosting Decision Tree (BDT) [36] for the C1 classifier. The number of decision trees ranges from 20 to 600. The model is implemented using the scikit-learn package. The metric used to evaluate the classification performance is the area under the ROC curve. The higher the score, the higher the classification power. Fig. 7.8 shows the relationship between the area under ROC curve and the number of trees. The classification performance becomes stable when the number 77 of decision trees is above 500. We use a BDT with 600 trees as the final model. Fig. sig 7.9 shows the importance of variables for classification. Eextra and |pl | are the most powerful variables to identify B → D(∗) (τ/l)ν from all types of backgrounds. For Eextra , the signal and normalization events usually have near-zero extra neutral energy, while background events have a wide distribution. The normalization decay sig has an energetic lepton produced by the D decay, leading to a higher |pl | value than signal events, as well as all types of backgrounds. The output of the BDT score p is transformed using a logit function: z1 = logit(p) = log p . 1−p The z1 distribution for all types of events is shown in Fig. 7.10. Signal and normalization events tend to have higher z1 values than backgrounds. Figure 7.8: Area under the ROC curve for BDT classifiers with different numbers of trees. 78 Figure 7.9: Importance of each variable for learning the C1 classifier. Figure 7.10: z1 distribution for signal, normalization, D∗∗ lν, B B̄ combinatorial, and continuum events. 79 C2 classifier The C2 Classifier aims to classify signal events and normalization events. Similar to the C1 classifier, we divide the sample into training and validation samples. The training sample is first used to train the classifier, the validation sample is then applied to evaluate the performance of the classifiers. Signal events are labeled positive and normalization events are labeled negative before training. We use the same variables used for the C1 classifier, with the addition of the following quantity: sig • cos θ B−D(∗) l : Cosine of the angle between the 3-momentum of the Bsig and the 3-momentum sum of its D and lepton daughters. The histogram of these variables for signal and normalization events are shown in Fig. 7.11. Figure 7.11: Histograms of variables used for the C2 classifier. Similar to the C1 classifier, we use a BDT for the C2 classifier. The number of decision trees ranges from 100 to 600. The classification performance is stable when the number of decision trees is above 100, as shown in Fig. 7.12. We use a BDT with 600 trees as the final model. Fig. 7.13 shows the importance of the sig sig variables for classification: cos θ B−D(∗) l and |pl | are the most powerful variables to sig distinguish signal from normalization events. For cos θ B−D(∗) l , normalization events have only a single neutrino, and the value should be from -1 to 1. However, signal 80 events tend to have more negative values, due to the presence of three neutrinos in the final state. The e/µ produced from the τ decays in signal events have a softer sig |pl | spectrum than the leptons produced from D decays for normalization events. The output of the BDT score is transformed to z2 using a logit function. The z2 distribution for all types of events is shown in Fig. 7.14. Signal events tend to have a higher z2 score, while normalization events tend to have a lower z2 score. Backgrounds have z2 scores in between. Figure 7.12: The area under the ROC curve for BDT classifiers with different number of trees. 81 Figure 7.13: Importance of each variable for learning the C2 classifier. Figure 7.14: z2 distribution for signal, normalization, D∗∗ lν, B B̄ combinatorial, and continuum events. 82 7.6 Signal Extraction In this section we describe how the signal yields are estimated. Our procedure includes two steps, and is applied to each subset individually to extract the corresponding signal yields: 1. Density Estimation: estimate the 2-dimensional density f (z1, z2 ) for each type of event. 2. Maximum Likelihood Estimation: estimate signal yields based on estimated densities, by solving a convex optimization problem. Density estimation Suppose a random variable X follows an unknown distribution f . If we have N observations of X, the objective of the density estimation procedure is to estimate f based on the N observations. Many approaches have been developed, both parametric and non-parametric, to solve this problem. One popular non-parametric approach is histogramming. Histogramming divides the range of X into m equallysized bins, and assigns each observation to one of these bins. The normalized counts of each bin gives us the estimated density f . Kernel density estimation (KDE) [37] is a non-parametric approach, which estimates N of X, the KDE estimate of f f via data smoothing. Given N observations {xi }i=1 is: N 1 Õ x − xi ), K( fˆ(x) = N h i=1 h where K is a kernel function, h is the bandwidth parameter to control the level of 1 2 smoothing. General choices of kernel function can be Gaussian (K(ν) = √1 e− 2 ν ) 2π or Epanechnikov (K(ν) = 34 (1 − ν 2 )), though the kernel function K can be any nonnegative bounded function satisfying 1. ∫1 −1 K(ν)dν = 1 2. K(ν) = K(−ν) ∫1 3. −1 ν 2 K(ν)dν < ∞. 83 The choice of bandwidth parameter h involves minimizing the mean squared error (MISE) between the estimated fˆ and f : M ISE = ∫ ( fˆ − f ) = ∫ 2 f −2 ∫ ˆ2 fˆ · f + ∫ f 2, which can be evaluated by summing over the observed data points of X. Practically, we use cross validation to choose the optimal bandwidth, which gives the same asymptotic accuracy: CV(h) = where fˆ−i (xi ) = N h bution from point i. ∫ fˆ2 − 2N −1 Õ fˆ−i (xi ), i xi −x j i, j K( h ) is used to evaluate MISE by removing the contri- 1 Í In this analysis, we use an adaptive KDE to evaluate the 2-dimensional probability density function (PDF) of each event type. Adaptive KDE changes the bandwidth parameters depending on the densities. In lower data density regions, adaptive KDE chooses wider bandwidths to reduce the effect of outliers on the overall KDE, and provides better overall performance. The adaptive KDE of f is: N Õ 1 x − xi ˆf (x) = 1 K( ) N i=1 hi hi hi = h × λi G α λi = ( ) , f˜(xi ) (7.6) where f˜(xi ) is a pilot estimate of f (xi ), λi is local bandwidth factor to control the local smoothing, G is the geometric mean of f (xi ), α is a sensitivity parameter between 0 and 1. We normally set its values to be 0.5. We use an Epanechnikov function as the kernel function, and use grid search to find the optimal bandwidth h. The KDE using an Epanechnikov kernel is claimed to have better mean integrated squared error under the same amount of data, and is one of the popular kernels used in practice. We use the bbrcit_kde library to reduce the time complexity of performing KDE from O(n2 ) to O(n). It applies kd-tree [38] algorithm and applies GPU computation to speedup the KDE computation. The learned densities for all event components for each of four subsets are shown in Fig. 7.15, 7.16, 7.17, and 7.18. B̄ → Dτ −ν̄τ 0.06 0.05 4 2 0 0. 01 B̄ → D∗τ −ν̄τ 0.06 0.05 2 0.05 0 −2 0.04 0.01 0.03 2 0.0 0.03 6 4 0.04 0.04 0. 0 005 6 8 0.03 6 0.04 8 z1 score z1 score 84 −2 0.03 0.02 0.02 0.02 −4 −4 0.01 −6 −4 0 −2 2 4 6 −8 −8 8 z2 score 8 6 B̄ → Dl −ν̄l 0.16 0.14 4 0 6 0.10 0.08 0.08 0.10 0.06 0 0.04 −2 −4 0.04 −4 −6 0.02 −6 −4 0 −2 2 4 6 −8 −8 8 z2 score 8 0.035 B̄ → D∗∗(τ /l)−ν̄ 0.030 4 0.025 2 z1 score z1 score 4 6 8 z2 score 0.12 B̄ → D∗l −ν̄l 2 0.06 6 2 0.10 0.12 0.05 −6 0 −2 8 −2 −8 −8 −4 4 0.1 0 0.12 0 .08 0.03 2 −6 0.08 0.0 2 z1 score −8 −8 0.01 −6 z1 score −6 0.06 0.04 0.02 −6 −4 0 −2 2 4 6 8 z2 score 8 6 0.055 Other bkg. 0.050 0.045 4 0.040 0.035 2 0.020 −2 0.03 0.010 0.0 2 −4 0.005 −6 −8 −8 −6 −4 0.025 0.015 0.0 1 −2 0 2 4 6 8 z2 score 0.020 −4 0.015 4 0.0 0.02 −6 −8 −8 0.05 0.03 −2 0.030 0 0.01 0.02 0.0 2 0 0.010 0.01 0.005 −6 −4 −2 0 2 4 6 8 z2 score Figure 7.15: KDE learned densities for event components in the D+ l subset. 8 B̄ → Dτ −ν̄τ 6 0.06 4 z1 score z1 score 85 8 0.06 B̄ → D∗τ −ν̄τ 6 0.05 4 0.05 0.02 2 0 0.01 00.0.04 5 0 2 0.0 0.03 −2 3 0.0 0.03 0.04 2 0.04 0.0 0. 0.056 4 0.0 1 −2 0.03 0.02 0.02 0.02 −4 −4 0.01 z1 score −8 −8 −6 −4 −2 0 2 4 6 −8 −8 8 z2 score 8 B̄ → Dl −ν̄l 6 0.18 0.16 0 0.15..1029 4 3 0.0 0.06 0.12 0.10 0 0.06 −4 0.04 8 z2 score 8 6 B̄ → D∗∗(τ /l)−ν̄ 0.028 0.024 4 0.0 1 2 3 0.0 0.0 2 0.01 0 0.12 0.10 0.05 0.08 0.02 z1 score z1 score 6 −6 −4 6 2 0.016 0 −6 0.004 −6 −2 0 2 4 6 8 z2 score 6 8 z2 score 0.045 0.040 0.035 0.030 0.025 0.020 −8 −8 0.015 0.04 0.02 −4 −4 4 Other bkg. −2 0.008 −6 2 0.02 −4 −8 −8 0 4 0.020 0.012 0.0 1 −2 8 0.03 −2 0.16 0.04 −8 −8 4 8 z2 score −4 0.00 2 6 0.06 0.02 0 4 −2 −8 −8 −2 2 0.14 −6 −4 0 00.1 .100. 2 08 −6 −6 −2 B̄ → D∗l −ν̄l 6 2 0.08 −2 −4 8 4 0.14 0 −6 0.03 2 0.01 −6 z1 score −6 0.03 0.010 0.01 0.005 −6 −4 −2 0 2 4 6 8 0.000 z2 score Figure 7.16: KDE learned densities for event components in the D0 l subset. 6 B̄ → D∗τ −ν̄τ 0.054 0.22 8 6 B̄ → Dl −ν̄l 0.20 00..12 6 0.16 0.18 0.048 4 4 0.0 3 0.10 −2 0.08 0.018 −4 −6 −4 −2 0 2 6 4 8 0.18 B̄ → D∗l −ν̄l 4 0.16 0.14 0.01.09 2 −6 −4 0 2 4 6 8 z2 score 8 6 0.036 B̄ → D∗∗(τ /l)−ν̄ 0.032 0.028 0.10 2 0.01 0 3 0.0 0.03 0.06 −2 0.01 0.08 −2 −2 4 0.12 0 0.02 −8 −8 z2 score 8 0.04 −6 0.006 −8 −8 0.06 −4 0.012 −6 2 0.12 0 0.024 −2 6 0.14 0.04 0.030 0.0 2 0.0 1 0 2 0.036 0.0 0.05 4 0.0 8 0.042 2 z1 score z1 score 8 z1 score z1 score 86 0.06 0.024 0.020 0.03 0.02 0.016 0.02 0.012 −4 0.04 −4 0.008 −6 0.02 −6 0.004 −6 −4 −2 0 2 z1 score −8 −8 6 4 −8 −8 8 z2 score −6 −4 −2 0 2 4 6 8 0.000 z2 score 8 6 Other bkg. 0.036 0.032 4 0.028 2 0 0.024 0.0 1 0.020 0.02 0.016 −2 4 0.0 01 0. −4 0.0 3 0.02 0.012 02 0. 0.008 −6 −8 −8 0.004 −6 −4 −2 0 2 4 6 8 0.000 z2 score Figure 7.17: KDE learned densities for event components in the D∗+ l subset. B̄ → D∗τ −ν̄τ 6 0.06 0.05 4 0.01 2 0 0.03 −2 0.22 0.20 0.18 00.12 .16 2 0.16 0.14 0.04 0.12 0 0.03 0.02 B̄ → Dl −ν̄l 6 4 0.04 0.0 0.0 0.5064 8 0.08 8 z1 score z1 score 87 0.10 −2 0.08 0.02 −4 0.01 −6 −6 −4 −2 0 2 6 4 z2 score 0.18 B̄ → D∗l −ν̄l 4 0.16 0.14 09 0. 0.12 0.1 5 2 0.02 −8 −8 8 8 6 0.04 −6 z1 score −8 −8 z1 score 0.06 −4 −6 −4 −2 0 2 4 6 8 z2 score 8 B̄ → D∗∗(τ /l)−ν̄ 6 0.028 0.024 4 0.12 0.10 0 0.020 1 0.0 3 0.0 0.0 3 2 0.06 0.0 1 0 0.016 0.08 0.06 −2 −4 0.04 −4 0.02 −6 −6 −6 −4 −2 0 2 z1 score −8 −8 6 4 0.008 0.004 −8 −8 8 z2 score 0.012 1 0.0 −2 0.02 −6 −4 −2 0 2 4 6 8 z2 score 8 6 0.054 Other bkg. 0.048 4 0.042 2 0.036 0.030 0 −2 0.02 −4 0.01 0.024 0.04 0.05 0.018 0.03 0.012 −6 −8 −8 0.006 −6 −4 −2 0 2 4 6 8 z2 score Figure 7.18: KDE learned densities for event components in the D∗0 l subset. Maximum likelihood estimation We extract the signal and normalization yields from each subset (channel_label) individually, and then combine them to get the overall yields. Given a subset with C components and the corresponding learned densities f j (z1, z2 ) ( j = 1, 2, ..., C), the z = (z1, z2 ) distribution of the subset can be written as Prob(z) = C Õ p j f j (z), j=1 where p j is the proportion of j’s component. If we observe N events in the subset N , the likelihood is {z}i=1 88 C N Õ Ö © ª p j · f j (zi )® . L= Prob(zi ) = ­ i=1 i=1 « j=1 ¬ The MLE estimate of component proportions p j ( j = 1, 2, ..., C) can be obtained by solving the following convex optimization problem: N Ö C N Õ Ö © ª p j · f j (zi )® ­ p1,p2,...,pC i=1 « j=1 ¬ Õ s.t. p j = 1. max (7.7) j B̄ → Dτ −ν̄τ B̄ → D∗τ −ν̄τ B̄ → Dl −ν̄l B̄ → D∗l −ν̄l 0.175 0.150 Density Density We use the CVXOPT [39] python package to implement the above estimation. The procedure is applied to the generic MC to evaluate its performance. To estimate the statistical sensitivity, we bootstrap the generic MC 900 times and solve the MLE for each bootstrap sample. The standard deviations of the 900 estimated proportions are used as absolute statistical uncertainties for the MLE estimation. The extracted yields and the absolute statistical uncertainties for all the four subsets are listed in Table 7.9. The projected fitting results on z1 and z2 scores are shown in Fig. 7.19 and Fig. 7.20. B̄ → D∗∗(τ /l)−ν̄ Other bkg. Data 0.150 0.125 0.125 0.100 0.100 0.075 0.075 0.050 0.050 0.025 0.025 0.000 0.25 B̄ → Dτ −ν̄τ B̄ → D∗τ −ν̄τ B̄ → Dl −ν̄l B̄ → D∗l −ν̄l B̄ → D∗∗(τ /l)−ν̄ Other bkg. Data 0.000 −6 B̄ → D∗τ −ν̄τ B̄ → Dl −ν̄l B̄ → D∗l −ν̄l −4 −2 0 2 4 6 −8 z1 score Density −8 Density 0.175 B̄ → D∗∗(τ /l)−ν̄ Other bkg. Data 0.20 −6 B̄ → D∗τ −ν̄τ B̄ → Dl −ν̄l B̄ → D∗l −ν̄l −4 −2 0 2 4 0 2 4 6 z1 score B̄ → D∗∗(τ /l)−ν̄ Other bkg. Data 0.15 0.20 0.15 0.10 0.10 0.05 0.05 0.00 −8 0.00 −6 −4 −2 0 2 4 6 z1 score −8 −6 −4 −2 6 z1 score Figure 7.19: Comparison of z1 score distributions of the data with the projections of fit results for (upper left) D+ l, (upper right) D0 l, (bottom left) D∗+ l, and (bottom right) D∗0 l subsets. 89 Subset Component B → Dτν B → D∗ τν B → Dlν B → D∗ lν B → D∗∗ lν Other Bkgs B → Dτν B → D∗ τν B → Dlν B → D∗ lν B → D∗∗ lν Other Bkgs B → D∗ τν B → Dlν B → D∗ lν B → D∗∗ lν Other Bkgs B → D∗ τν B → Dlν B → D∗ lν B → D∗∗ lν Other Bkgs D+ l D0 l D∗+ l D∗0 l Extracted Yield 117 ± 44 191 ± 57 1965 ± 110 1261 ± 146 3248 ± 233 7258 ± 182 726 ± 273 2452 ± 285 12012 ± 555 22808 ± 648 15594 ± 445 22065 ± 307 110 ± 17 614 ± 127 1246 ± 139 292 ± 54 665 ± 42 175 ± 25 479 ± 126 2154 ± 149 1387 ± 101 1415 ± 70 0.14 B̄ → Dτ −ν̄τ B̄ → D∗τ −ν̄τ B̄ → Dl −ν̄l B̄ → D∗l −ν̄l 0.12 Density Density Table 7.9: Fit results for the yields of all the components in four subsets. B̄ → D∗∗(τ /l)−ν̄ Other bkg. Data B̄ → Dτ −ν̄τ B̄ → D∗τ −ν̄τ B̄ → Dl −ν̄l B̄ → D∗l −ν̄l 0.20 B̄ → D∗∗(τ /l)−ν̄ Other bkg. Data 0.10 0.15 0.08 0.10 0.06 0.04 0.05 0.02 0.00 0.00 −6 −4 −2 0 2 4 6 −8 z2 score −6 −4 −2 0 2 4 6 z2 score 0.30 0.40 B̄ → D∗τ −ν̄τ B̄ → Dl −ν̄l B̄ → D∗l −ν̄l 0.35 Density Density −8 B̄ → D∗∗(τ /l)−ν̄ Other bkg. Data B̄ → D∗τ −ν̄τ B̄ → Dl −ν̄l B̄ → D∗l −ν̄l 0.25 B̄ → D∗∗(τ /l)−ν̄ Other bkg. Data 0.30 0.20 0.25 0.15 0.20 0.15 0.10 0.10 0.05 0.05 0.00 0.00 −8 −6 −4 −2 0 2 4 6 z2 score −8 −6 −4 −2 0 2 4 6 z2 score Figure 7.20: Comparison of z2 score distributions of the data with the projections of fit results for (upper left) D+ l, (upper right) D0 l, (bottom left) D∗+ l, and (bottom right) D∗0 l subsets. 90 Cross-check on PDFs modeling Cross-checks are performed to study the modeling of PDF shapes. We first check the PDF modeling of normalization events by comparing the PDFs of MC with data on the normalization enriched region. The normalization enriched region is defined in (z1, z2 ) space as z1 > 2, z2 < −4, in which 94% events are normalization decays. The MC and data comparison shown in Fig. 7.21 indicates the normalization PDF shapes are well modeled. B̄ → D(∗)τ −ν̄τ 0.7 0.5 B̄ → D(∗)τ −ν̄τ 0.4 B̄ → D∗∗(τ /l)−ν̄ BB̄ combinatorial Continuum B̄ → D(∗)l −ν̄l B̄ → D∗∗(τ /l)−ν̄ Density 0.5 Density BB̄ combinatorial Continuum B̄ → D(∗)l −ν̄l 0.6 0.4 0.3 0.3 0.2 0.2 0.1 0.1 2.5 3.0 3.5 4.0 4.5 1 0 −1 2.0 2.5 3.0 3.5 4.0 z1 score 0.0 −8.0 5.0 4.5 5.0 (Data - MC) / Data (Data - MC) / Data 0.0 2.0 −7.5 −7.0 −6.5 −7.5 −7.0 −6.5 −6.0 −5.5 −5.0 −4.5 −4.0 −6.0 −5.5 −5.0 −4.5 −4.0 0.5 0.0 −0.5 −8.0 z2 score Figure 7.21: MC/data comparison of z1 and z2 score distributions in normalization enriched region. We then check the PDF modeling of backgrounds (after on-peak background calibration) by comparing the PDFs of MC with data on background enriched region. The background enriched region is defined as −4 ≤ z1 ≤ −2, in which 97% events are either B → D∗∗ lν, B B̄ combinatorial, or continuum events. The MC and data comparison shown in Fig. 7.22 indicates the PDF shapes of backgrounds are also well-modeled. 1.0 B̄ → D(∗)τ −ν̄τ B̄ → D∗∗(τ /l)−ν̄ B̄ → D∗∗(τ /l)−ν̄ 0.125 0.6 0.4 BB̄ combinatorial Continuum B̄ → D(∗)l −ν̄l 0.150 Density Density B̄ → D(∗)τ −ν̄τ 0.175 BB̄ combinatorial Continuum B̄ → D(∗)l −ν̄l 0.8 0.100 0.075 0.050 0.2 0.025 −3.75 −3.50 −3.25 −3.00 −2.75 −2.50 −2.25 −2.00 0.25 0.00 −0.25 −4.00 −3.75 −3.50 −3.25 −3.00 z1 score −2.75 −2.50 −2.25 −2.00 0.000 −8 (Data - MC) / Data (Data - MC) / Data 0.0 −4.00 −6 −4 −2 −6 −4 −2 0 2 4 6 8 0 2 4 6 8 1 0 −1 −8 z2 score Figure 7.22: MC/data comparison of z1 and z2 score distributions in backgrounds enriched region. 91 7.7 Systematic Uncertainties Procedure Many quantities are required to precisely measure R(D(∗) ), for instance, the PDF shape of each event type, the physical parameter setup for the Monte Carlo simulations, and the detector efficiencies. We use our best knowledge of this information during the measurement, however, this information has associated uncertainties. Consequently, our measurement might be biased due to these sources of uncertainty. This section describes the systematic uncertainties that are important in this analysis, and how we evaluate them. Different procedures are developed to evaluate each source of uncertainty. Uncertainties from form factors, branching fractions of B̄ → D∗∗ (τ/l)− ν̄, b → c c̄ decays, and D meson decays are evaluated based on the delta method [40]. We fluctuate every uncertain parameter up and down according to their 1σ uncertainties. Each time the parameter is changed, we re-evaluate the kernel densities and repeat the fit with the new PDFs. The changes in R(D(∗) ) are quoted as the corresponding systematic uncertainty. Uncertainties due to the limited size of the MC sample and background calibration are evaluated using a bootstrap algorithm [41]. B → D(∗) lν form factors One source of uncertainty comes from the model behind the Monte Carlo simulation of B meson decays, in which form factors are required. However, the form factors used for the MC simulation have associated uncertainties, which contribute to the differential decay rate of B → Dτν through angular distributions, and therefore the shapes of the PDFs. To transform the uncertainties on the form factors to R(D(∗) ), we first transform the default model i, which is used to generate simulation events, to an updated parameterization model j, and use model j to evaluate the corresponding systematics. The default form factor model used to simulate B → Dlνl , l = e, µ, τ events is the 92 CLN [24] model. It uses dispersion relations to expand the form factors [42] V1 (ω) = 1 − 8ρ2 z + (51ρ2 − 10)z 2 − (262ρ2 − 84)z3 + O(z4 ) V1 (1) S1 (ω) = (1 + ∆(ω)) V1 (ω) ρ2 = 1.186 ± 0.054 (7.8) V1 = 1.0816 ∆ = 1.0. The expansions of the CLN model used to simulate B → D∗ lν, l = e, µ, τ and baseline parameter settings are:: h A1 (ω) = 1 − 8ρ2 z + (53ρ2 − 15)z2 − (231ρ2 − 91)z3 + O(z4 ) F1 R0 (ω) = R0 (1) − 0.11(ω − 1) + 0.01(ω − 1)2 + O(ω3 ) R1 (ω) = R1 (1) − 0.12(ω − 1) + 0.05(ω − 1)2 + O(ω3 ) R2 (ω) = R2 (1) + 0.11(ω − 1) − 0.06(ω − 1)2 + O(ω3 ) F1 = 0.921 (7.9) ρ2 = 1.207 ± 0.026 R0 = 1.14 R1 = 1.401 ± 0.033 R2 = 0.854 ± 0.02. To evaluate the systematic uncertainty due to the B → D(∗) (l/τ)νl form factors, we measure the central values of R(D(∗) ) using the world average values of parameters in form factor model. Then we change the form factor parameters with the corresponding uncertainties by ±σ. The difference between the values obtained before and after fluctuation are listed as systematic. The results are shown in Table 7.10. Source B → D(l/τ)ν form factor B → D∗ (l/τ)ν form factor ∆R(D) (%) ∆R(D∗ ) (%) 0.42 0.13 0.92 0.31 Table 7.10: Evaluated relative systematic uncertainties from the B → Dlν and B → D∗ lν form factors. 93 B → D∗∗ (l/τ)νl form factors The systematic uncertainties from semileptonic B meson decays involving D∗∗ are also evaluated, where D∗∗ denotes the L = 1 excitations of the ground state D meson: D0∗, D1, D10 , and D2∗ . In this case, we use the LLSW B1 model [43] as baseline, fluctuate it to the LLSW B2 [43] model and take the difference between the two models as systematic uncertainties. The resulting systematic uncertainty is shown in Table 7.11. Source B → D∗∗ (l/τ)ν form factor ∆R(D) (%) ∆R(D∗ ) (%) 0.48 0.18 Table 7.11: Evaluated relative systematic uncertainties from the B → D∗∗ lν form factors. B → D(∗) lν branching fractions Branching fractions from semileptonic B meson decays are another source of systematic uncertainties, as they affect the relative abundance of the decays and thus the density distributions of event components, as well as the estimate of signal and normalization efficiencies. Some of the branching fractions used for the MC generation are different than the current world average. Therefore, before evaluating the systematic uncertainties associated with the branching fractions, we assign correction factors to re-weight the MC events: ω(x) = ωW .A. (x) , ωDECAY.DEC (x) where ω(x) is the assigned weight for event x, ωW .A. is the world average value of event x’s branching fraction, and ωDECAY.DEC (x) is the value used in the MC simulation. The B → D(∗) lν decays are one of the dominant B meson decays in this analysis; we corrected and fluctuated their uncertainties to evaluate the systematic uncertainties according to the values listed in Table 7.12. The resulting systematic uncertainties are listed in Table 7.13. B → D∗∗ (l/τ)ν branching fractions Semileptonic decays of B → D∗∗ lνl and B → D∗∗ τντ are an important background as they have topologies to the signal events. In general, D∗∗ is defined as any excited 94 Decay Mode MC Simulation Value World Average Value Uncertainty B+ → D∗0 µ+ ν µ B+ → D∗0 e+ νe B+ → D0 µ+ ν µ B+ → D0 e+ νe B0 → D∗− µ+ ν µ B0 → D∗− e+ νe B0 → D− µ+ ν µ B0 → D− e+ νe 0.0617 0.0617 0.0224 0.0224 0.057 0.057 0.0207 0.0207 0.0566 0.0566 0.0235 0.0235 0.0506 0.0506 0.0231 0.0231 0.0022 0.0022 0.0009 0.0009 0.0012 0.0012 0.0010 0.0010 Table 7.12: Decay modes fluctuated to evaluate the B → D(∗) lν branching fractions on generic MC. The world average values are from [34]. Source B → D(∗) lν branching fraction ∆R(D) (%) ∆R(D∗ ) (%) 0.47 0.38 Table 7.13: Evaluated relative systematic uncertainties from B → D(∗) lν branching fractions on MC. charmed meson states that is not in the ground-state 1S doublet [44]. In this analysis, we consider the following three D∗∗ decay types: • Resonant D∗∗ (1P) state, which includes the four lightest excited charmed meson states D0∗ (0+ ), D1 (1+ ), D10 (1+ ) and D2∗ (2+ ). • Resonant D∗∗ (2S) state, the radially-excited modes of D(∗) . • Non-resonant D∗∗ states. We allow unmeasured D∗∗ states whose final states consist of a D(∗) and one pion. We estimate the systematic uncertainties from both B → D∗∗ lνl and B → D∗∗ τντ for these three types of D∗∗ decays. We do not have MC samples for some B → D∗∗ τντ decays, so we use a 3-body phase space model for estimation. Table 7.14 shows all the decay modes taken into consideration. We do not consider other charmed states heavier than D(∗) (2S) as their smaller phase space suppresses the branching fractions. 95 B → D∗∗ lνl B → D∗∗ τντ Resonant D∗∗ (1P) MC MC Resonant D∗∗ (2S) MC Phase Space Non-resonant B → D(∗) πlνl MC Phase Space Table 7.14: B → D∗∗ (D(∗) π)(l/τ)ν decays and method used to estimate the systematic uncertainties. "MC" means the uncertainties can be directly estimated using MC sample. "Phase Space" means we do not have the corresponding MC samples, therefore a phase-space-based estimation is applied. To evaluate the resonant D∗∗ (1P) for B → D∗∗ (l/τ)ν decays, the following decays in the generic MC sample are fluctuated as listed in Table 7.15. The corresponding systematic uncertainties are obtained in Table 7.16. 96 Decay Mode 0 B+ → D1 e+ νe ∗0 B+ → D0 e+ νe ∗0 B+ → D2 e+ νe 00 B+ → D1 e+ νe 0 B+ → D1 µ+ ν µ ∗0 B+ → D0 µ+ ν µ ∗0 B+ → D2 µ+ ν µ 00 B+ → D1 µ+ ν µ 0 B+ → D1 τ + ντ ∗0 B+ → D0 τ + ντ 00 B+ → D1 τ + ντ ∗0 B+ → D2 τ + ντ B0 → D2∗− e+ νe 0 B0 → D1− e+ νe B0 → D0∗− e+ νe B0 → D1− e+ νe B0 → D1− µ+ ν µ B0 → D0∗− µ+ ν µ 0 B0 → D1− µ+ ν µ B0 → D2∗− µ+ ν µ B0 → D1− τ + ντ B0 → D0∗− τ + ντ 0 B0 → D1− τ + ντ B0 → D2∗− τ + ντ MC Simulation Value Setup Value Uncertainty 0.0056 0.0096 0.001 0.0049 0.0044 0.0008 0.003 0.003 0.0004 0.009 0.002 0.0005 0.0056 0.0096 0.001 0.0049 0.0044 0.0008 0.003 0.003 0.0004 0.009 0.002 0.0005 0.0013 0.001 0.00014 0.0013 0.0004 0.00015 0.002 0.00012 0.00005 0.002 0.0023 0.0083 0.0045 0.0052 0.0052 0.0045 0.0083 0.0023 0.0013 0.0013 0.002 0.002 0.00021 0.0028 0.0019 0.00408 0.0089 0.0089 0.00408 0.0019 0.0028 0.0009 0.0003 0.00017 0.00013 0.00004 0.0004 0.00046 0.00074 0.000911 0.000911 0.00074 0.00046 0.0004 0.00013 0.00014 0.00005 0.00004 Table 7.15: Decay modes fluctuated to evaluate the resonant B → D∗∗ (1P)(l/τ)ν branching fractions on generic MC. The setup values are from [34, 45]. Source ∗∗ Resonant B → D (1P)(l/τ)ν branching fraction ∆R(D) (%) ∆R(D∗ ) (%) 1.87 0.35 Table 7.16: Evaluated relative systematic uncertainties from resonant B → D∗∗ (1P)(l/τ)ν branching fractions. The samples simulated B → D∗∗ lν decays are used to evaluate the systematic uncertainties from resonant B → D∗∗ (2S)lν decays, as listed in Table 7.17. The masses for D∗∗ (2S) and D∗∗ (2S)∗ are around 2.47 GeV and 2.7 GeV. These events are 97 then reconstructed, selected, and assigned (z1, z2 ) scores using the same procedure as for the generic MC. It is well known that there is a discrepancy of 1.5% between the inclusive branching fraction of semileptonic B decays and the sum of exclusive branching fractions. Conservatively, we use this sample to make up the 1.5% discrepancy. The densities of each event type are then re-evaluated. The difference between the measured R(D(∗) ) on whether the B → D∗∗ (2S)lν decays are included are taken as a systematic uncertainty. The resulting systematic uncertainty is shown in Table 7.18. Decay Type + B → D∗∗ (2S)(Dππ)`ν B0 → D∗∗ (2S)(Dππ)`ν B+ → D∗∗ (2S)(D∗ ππ)`ν B0 → D∗∗ (2S)(D∗ ππ)`ν B+ → D∗∗ (2S)∗ (Dππ)`ν B0 → D∗∗ (2S)∗ (Dππ)`ν B+ → D∗∗ (2S)∗ (D∗ ππ)`ν B0 → D∗∗ (2S)∗ (D∗ ππ)`ν # event [106 ] 6.776 6.826 6.530 6.769 6.369 6.552 6.425 6.616 BABAR Dataset Name SP-11461-R24 SP-11467-R24 SP-11462-R24 SP-11468-R24 SP-11463-R24 SP-11469-R24 SP-11464-R24 SP-11470-R24 Table 7.17: MC samples used for assessing the gap between inclusive and the sum of exclusive B → Xc `ν branching fractions. Source B → D∗∗ (2S)lν branching fraction ∆R(D) (%) ∆R(D∗ ) (%) 0.58 0.56 Table 7.18: Evaluated relative systematic uncertainties from the resonant B → D∗∗ (2S)lν branching fractions. For non-resonant B → D∗∗ lν decays, we fluctuated both B → D(∗) π + lν and B → D(∗) π 0 lν decay modes simultaneously and the corresponding systematic uncertainties are listed in Table 7.19. Source Non-resonant B → D(∗) πlν branching fraction ∆R(D) (%) ∆R(D∗ ) (%) 2.05 1.46 Table 7.19: Evaluated relative systematic uncertainties from non-resonant B → D(∗) πlν branching fractions. 98 We do not have MC samples for resonant B → D∗∗ (2S)τν and non-resonant B → D(∗) πτν decays. However, since B → D∗∗ τν and B → D∗∗ lν have similar topology in this analysis, we will use the branching ratio between B(B → D∗∗ τν) and B(B → D∗∗ lν) to estimate the systematic uncertainties. We estimate the relative branching ratio R(D∗∗ ) with the available phase space Φ R(D∗∗ ) = B(B → D∗∗ τν) Φ(B → D∗∗ τν) ≈ . B(B → D∗∗ lν) Φ(B → D∗∗ lν) The phase space of the three body decay B → Mlν is given by the integral d 3 p M d 3 pl d 3 p ν 4 δ (pB − p M − pl − pν ) 2E M 2El 2Eν s s ∫ (mB −mM )2 2 4 2m m m2 − m2M + q2 )2 − q2 . ∝ dq2 1 − 2l + 4l ( B 2mB q q ml2 Φ(B → Mlν) = ∫ (7.10) We calculate R(D∗∗ ) as well as the Φ(B → D∗∗ lν) integrals numerically for a wide range of m M , as shown in Fig. 7.23. The mass of the discovered excited D mesons are roughly in the range of [2300, 3300] MeV, so we use R(D∗∗ ) = 0.13 at mD∗∗ = 2300 MeV to conservatively estimate the systematic uncertainty from resonant B → D∗∗ (2S)τν decays. The results are shown in Table 7.20. For heavier excited D mesons, the phase space decreases and thus has a smaller effect on the measurement. HFLAV 2+ D ∗ (3000)0 2 1− Natural 3− D ∗ (2760)0 1 D ∗ (2760)± 3 D ∗ (2760)0 3 1− 2− 3− 1+ D(2740)0 D(2750)0 1− D ∗ (2650)0 D ∗ (2680)0 0− Natural 2+ D ∗ (2640)± 1+ D(2600)0 D(2600)± 0+ Unnatural 2021 D(2580)0 3700 3600 3500 3400 3300 3200 3100 3000 2900 2800 2700 2600 2500 2400 2300 D ∗ (2400)0 0 D ∗ (2400)± 0 D1 (2430)0 D (2420)0 1 D (2420)± 1 D ∗ (2460)0 2 D ∗ (2460)± 2 D(2550)0 Mass [MeV/c 2 ] 99 (b) 1e13 0.200 2.5 0.175 0.150 Φ(B→D * * lν) 2.0 R(D**) 0.125 0.100 1.5 0.075 1.0 0.050 0.025 0.000 1800 2000 2200 2400 2600 2800 mD * * / MeV 3000 3200 0.5 1800 2000 2200 2400 2600 2800 3000 3200 mD * * / MeV Figure 7.23: (upper) The current mass ranges for excited D meson states. (bottom left) The R(D∗∗ ) for different excited D meson masses. (bottom right) The Φ(B → D∗∗ lν) for different excited D meson masses. Source Resonant B → D∗∗ (2S)τν Non-resonant B → D(∗) πτν ∆R(D) (%) ∆R(D∗ ) (%) 0.08 0.07 0.27 0.19 Table 7.20: Evaluated relative systematic uncertainties from resonant B → D∗∗ (2S)τν and non-resonant B → D(∗) πτν branching fractions. D meson decay branching fractions The uncertainties from D meson decay branching fractions are also considered. Since performing fluctuations on all D meson decay channels is impractical, we 100 simply fluctuate the most common decay D → Kππ. Specifically, for each event involving a D → Kππ decay, we assign a factor of ω = ωW .A./ωDECAY.DEC = 0.992. Table 7.21 shows the resulting systematic uncertainties. Source D → Kππ branching fraction ∆R(D) (%) ∆R(D∗ ) (%) 0.84 0.70 Table 7.21: Evaluated relative systematic uncertainties from the D → Kππ branching fractions. Υ(4S)) → B B̄ branching fractions The branching fractions of B(Υ(4S)) decaying to charged or neutral B mesons are B(Υ(4S) → B+ B− ) = 51.4% ± 0.6%, and B(Υ(4S) → B0 B̄0 ) = 48.6% ± 0.6%. To evaluate the effect of the uncertainty of the neutral/charged B meson ratio, we fluctuate the relative number of events between B0 B̄0 and B+ B− in the amount of 1.1% and propagate it to the uncertainty of R(D(∗) ). The results are in Table 7.22. Source B(Υ(4S)) ∆R(D) (%) ∆R(D∗ ) (%) 0.48 0.43 Table 7.22: Evaluated relative systematic uncertainties from B(Υ(4S)) decays to charged or neutral B mesons. Lepton efficiency We assign a 1% relative uncertainty per electron or muon in an event, and assign a 1.1% relative uncertainty per kaon in an event. There are two light leptons and two kaons for a signal or normalization event, one from the Btag and the other from Bsig . However, most sources of PID efficiency uncertainties will cancel out, since we measure the ratio between signal and normalization event. The only source left is the signal side lepton PID efficiency, because the signal side lepton momentum varies for signal and normalization events. norm ∼ N (µ , σ 2 ), then If we assume PID ∼ N (µ x , σx2 ), PID y y signal    2 σ  σ  σx/y 2  σx  2 σy y x = + −2· · . x/y x y x y Based on [46, 47], the corresponding R(D(∗) ) uncertainties are 0.29% and 0.40%. 101 Neutral pion efficiency The soft π 0 efficiency uncertainty does not cancel in the ratio. Studies have been performed by the BABAR collaboration to measure the difference of the π 0 efficiency in data and simulation. Following BAD [48], an average π 0 correction ωπ0 = 0.958 ± 0.009 is applied, and the corresponding uncertainties are propagated to the results of R(D∗ ). Although the soft pion efficiency is momentum dependent, the momentum for nearly all soft pions is in the range of [0.1, 1] GeV in this analysis. Given the momentum-dependent correction factor as shown in Fig. 7.24, this is a second-order effect and can be neglected. Figure 7.24: The π 0 efficiency correction values in dependence of the π 0 momentum. The figure is from [48]. PDF shape of on-peak backgrounds Among the three types of backgrounds, the continuum background is characterized using off-peak data. However, the shape of B → D∗∗ lν and combinatorial B B̄ background are characterized using generic MC, which may be different from what happened in the detector. We denote the latter two backgrounds as B B̄ backgrounds. To compare the shape difference of B B̄ backgrounds between MC and data, we construct a sideband in which the data consists only of background events. The criteria to define the sideband is: S1: Eextra > 1.0 GeV 102 sig S2: |pl | < 0.4 GeV We assume that the signal and normalization events are negligible in the sideband region (estimated proportions are 0.4% for signal events and 3.8% for normalization events). A sample of 1% of on-peak data is used to demonstrate and calibrate the difference between MC and data, as shown in Fig. 7.25. Since the continuum background is well modeled using off-peak data, the discrepancy in the figure is a result of mis-modeling of B B̄ background. 0.6 BB background Data 0.5 1.0 BB background Data 0.8 0.4 0.6 0.3 0.4 0.2 0.2 0.1 0.0 4 2 0 Z1 2 4 6 0.0 3 2 Z2 1 0 1 Figure 7.25: z1 and z2 distribution of sidebands for generic MC and sideband data. The difference between data and histogram indicates the discrepancy between MC and data. The discrepancy of the (z1, z2 ) distribution between MC and data may introduce a bias on this measurement, therefore, we apply a correction factor based on the 2-dimensional shape difference as follows: 1. Let g(z1, z2 ) be the (z1, z2 ) distribution of on-peak data in the sideband region. Since the signal and normalization events are negligible, its density can be decomposed as the sum of B B̄ background and continuum components: g(z1, z2 ) = pB B̄ gB B̄ (z1, z2 ) + pcont gcont (z1, z2 ) pB B̄ + pcont = 1. 2. Let fB B̄ (z1, z2 ) be the density function of the sideband MC B B̄ background. 3. The correction factors are: ω(z1, z2 ) = gB B̄ (z1, z2 ) g(z1, z2 ) − pcont gcont (z1, z2 ) = . fB B̄ (z1, z2 ) pB B̄ fB B̄ (z1, z2 ) (7.11) 103 The systematic uncertainty from these correction factors is estimated using a bootstrap technique. Since the variance of ω(z1, z2 ) is due to the limited sample, to model the PDF shapes, we bootstrap the sideband sample of B B̄ background, on-peak data, and off-peak data to capture its effect on the measurement. 1. We bootstrap the sideband sample of B B̄ background, on-peak data, and off-peak data 100 times. 2. For each bootstrap sample, we calculate Equation 7.11 and apply the sideband calibration for each B B̄ background events. 3. For each bootstrap sample, we use updated generic MC events to extract the signal and normalization yields and measure R(D(∗) ). 4. We take the sample standard deviation of the measured R(D(∗) ) as an estimate of the systematic uncertainty. The results are shown in Table 7.23. Source B B̄ background shape ∆R(D) (%) ∆R(D∗ ) (%) 0.73 0.37 Table 7.23: Evaluated relative systematic uncertainties from the calibration factor on the B B̄ background shape difference between generic MC and data. Sideband selection efficiency Our definition of sideband does not guarantee the same selection efficiency for B B̄ events and D∗∗ lν events. Therefore, the corresponding systematic is estimated. We apply the correction only for B B̄ events and D∗∗ lν events, and take the half of the discrepancy as an estimate of the systematic uncertainty. The results are shown in Table 7.24. Source Sideband Region ∆R(D) (%) ∆R(D∗ ) (%) 2.34 0.23 Table 7.24: Evaluated relative systematic uncertainties from the calibration factor on the B B̄ background shape difference between generic MC and data. 104 Correlation between the uncertainties on R(D) and R(D∗ ) The correlations between R(D) and R(D∗ ) are used to combine all the uncertainties together and compare with the theoretical predictions. We use the MC samples to estimate the correlations for additive systematic uncertainties, while the correlations for multiplicative uncertainties can be clearly derived. For the additive systematic uncertainties, a standard approach to estimate correlations is bootstrapping: we should bootstrap several samples and measure R(D(∗) ) for each sample, and then estimate the correlations using the sample results. However, the above approach is impractical in this analysis, due to the lengthy computation required to evaluate the kernel densities for each bootstrapping sample. Therefore, we use importance sampling as an approximate method. For a given systematic uncertainty source X with standard deviation σ, we assume it follows a Gaussian distribution X ∼ N (µ, σ 2 ), and affects the R(D(∗) ) measurements in the unknown form of r (∗) (X). We fluctuate X up and down on an amount of ±σ/2, ±σ, and then we measure R(D(∗) ) for each fluctuated sample. The correlation between R(D) and R(D∗ ) under the systematic uncertainty source X can be obtained by computing the weighted Pearson correlations between the following two samples: [r(X − σ), r(X − σ/2), r(X), r(X + σ/2), r(X + σ)] [r ∗ (X − σ), r ∗ (X − σ/2), r ∗ (X), r ∗ (X + σ/2), r ∗ (X + σ)], with the weighting [Prob[X = µ − σ], Prob[X = µ − σ/2], Prob[X = µ], Prob[X = µ + σ/2], andProb[X = µ + σ]], which are standard Gaussian densities at the corresponding ±σ/2, ±σ away from mean value. This approach is applied to estimate correlations for all the additive systematic uncertainties as well as the B(Υ(4S)) uncertainty. Among the multiplicative uncertainties, the PID efficiency, soft π 0 efficiency and tracking efficiency affect R(D) and R(D∗ ) equally, so it has a 100% correlation. The efficiency uncertainty due to limited statistics (MC Efficiency) is clearly uncorrelated. The uncertainty on B(τ → l − ν̄l ντ ) affects all channels equally, so it also has 100% correlation. The correlation on the statistical uncertainty is evaluated by fitting the MLE multiple times, and is found to be negatively correlated. Summary Table 7.25 summarizes all the uncertainties taken into consideration, as well as the correlations between the uncertainties on R(D) and R(D∗ ). The effective correlation 105 coefficient ρtot is calculated by adding the covariance matrices as follows: ! ! ∗ 2 Õ ρtot σtot σtot σi2 ρi σi σi∗ σtot , = ∗2 ∗ σtot ρtot σtot σtot ρi σi σi∗ σi∗2 i (7.12) where σi(∗) refers to uncertainties on R(D(∗) ), and i runs over each source of uncertainty. Source B → Dlν form factor B → D∗ lν form factor B → D∗∗ lν form factor B(B → D(∗) lν) B(b → c c̄) B(B → D∗∗ lν) B(D) PDF shapes MC statistics On-peak background calibration B(Υ(4S)) PID efficiency Soft π 0 efficiency B(τ → l − ν̄l ντ ) Systematic Total Statistical Uncertainty Total ∆R(D) (%) ∆R(D∗ ) (%) 0.42 0.13 0.92 0.31 0.48 0.18 0.47 0.38 0.34 0.13 2.83 1.60 0.84 0.70 4.12 4.37 2.45 0.44 0.48 0.43 0.29 0.40 0.84 1.25 0.16 0.16 5.86 4.96 19.8 9.9 20.65 11.07 Correlation -0.18 -0.19 -0.90 0.97 1 -0.97 -0.40 -0.15 -0.05 1 1 1 1 -0.20 -0.92 -0.82 Table 7.25: Summary of evaluated uncertainties (preliminary). 7.8 Results Our preliminary results are R(D) = 0.316 ± 0.062 ± 0.019 R(D∗ ) = 0.226 ± 0.022 ± 0.012, (7.13) where the first uncertainties are statistical and the second are systematic. The correlation between R(D) and R(D∗ ) is -0.82. The fitted yields for each event type in all of four subsets are listed in Table 7.9. The signal and normalization efficiencies are listed in Table 7.26. The comparison of our result with previous measurments is shown in Fig. 7.26 106 Figure 7.26: The comparison of preliminary results with previous measurments for R(D(∗) ). Event type B → Dτν B → D∗ τν B → Dlν B → D∗ lν Efficiency (%) 0.41 0.47 0.46 0.43 Table 7.26: Efficiencies for signal and normalization events. 7.9 Summary In summary, we measure R(D(∗) ) using a semileptonic tagging method and leptonic τ decays. Our preliminary results are R(D) = 0.316 ± 0.062 ± 0.019 and R(D∗ ) = 0.226 ± 0.022 ± 0.012, where the first uncertainties are statistical and the second are systematic. The measured R(D), R(D∗ ), and their combinations agree with the Standard Model predictions by 0.26σ, 1.10σ, and 1.51σ. This analysis is BABAR’s first measurement of R(D(∗) ) using a semileptonic tagging method. The measurement strategy imposes as minimal assumptions from Monte Carlo samples as possible. For instance, instead of fixing background yields, they are also treated as free parameters when fitting the data. By doing this, our strategy 107 can overcome potential simulation bias. It is also feasible because of the customer package developed to best model PDFs of decays with millions of MC samples. Due to this, our measurement strategy is unique and reliable. Our combined results on R(D(∗) ) agree with the Belle’s semileptonic tagging measurement by 2.23σ. If we average the two semileptonic tagging measurements together and compare with the SM prediction, the combined difference is about 0.46σ. Similar to previous measurements, this measurement is also statistically dominant. However, with future Belle II experiment target luminosity of 50 ab−1 , the statistical uncertainty on R(D(∗) ) measurements can be reduced to only few percent, making it comparable with the systematic uncertainties and better to probe the new physics. Part V Conclusion 108 109 Chapter 8 CONCLUSION This thesis focuses on collider searches of beyond Standard Model physics at the intensity frontier. With an integrated luminosity of 531 fb−1 collected at the BABAR experiment, both direct and indirect methods to probe for New Physics are performed with this unique experimental environment. We first briefly review the Standard Model, which is the current best theory to describe the interactions between fundamental particles, and has been extensively tested by a series of experiments. However, the Standard Model is not a complete theory. Besides the fact that it does not include gravity, we also summarize its limitation in both experimental and theoretical aspects. The main goal of high energy physics is to search and understand the beyond Standard Model physics. We then describe the details of the PEP-II accelerator and BABAR detector. It is an e+ e− collider operating at the center-of-mass energy around 10.58 GeV. The asymmetric design enables the B meson decay length to be measurable in the lab frame. We also briefly summarize the main components of the BABAR detector. We present our direct search for beyond Standard Model physics by looking for dark sectors. Dark sectors are new particle(s) interacting only feebly with ordinary matter mediated via portals, and have become an intriguing framework to explain the presence of dark matter in the Universe. While previous collider searches focus on identifying new mediators, we investigate the possibility of dark matter bound states as a probe for dark sector. This is also the first search for darkonium. In an absence of signal, we show that this search improves the existing constraints on the γ − A0 mixing strength over a significant fraction of dark photon masses below 1 GeV for large values of the dark sector coupling constant. We then present our indirect search for beyond Standard Model physics by precision measurement of semileptonic B meson decays. We measure R(D) and R(D∗ ), which are sensitive probes for BSM physics, using a semileptonic tagging method. The results are dominated by statistical uncertainties, and agree with the Standard Model prediction within 1.51σ. This measurement applies a data-driven approach and is therefore robust to potential bias from simulation. The unique strategy is enabled by our custom-developed software for fast kernel density estimation powered by GPU 110 technology, which enables us to best characterize signatures of decays with large samples. In the appendix, we demonstrate the idea of applying deep transfer learning algorithms to reduce systematic uncertainties. With the adversarial neural network architecture, the multivariate classifiers can be trained to be insensitive to the small variance of event distributions. This framework can be applied to most high energy physics analysis to reduce the systematics from mis-modeling of signal or background samples. 111 BIBLIOGRAPHY [1] “https://www.quantumdiaries.org/2014/03/14/the-standard-model-a-beautifulbut-flawed-theory”. [2] MS Zisman. “PEP-II: An Asymmetric B Factory: Conceptual Design Report”. In: SLAC Report 418 (1993), pp. 2–93. [3] D Boutigny et al. “The BABAR physics book: Physics at an asymmetric B factory”. In: (1999). [4] Gerard Jungman, Marc Kamionkowski, and Kim Griest. “Supersymmetric dark matter”. In: Physics Reports 267.5-6 (1996), pp. 195–373. [5] SA Abel et al. “Kinetic mixing of the photon with hidden U (1) s in string phenomenology”. In: Journal of High Energy Physics 2008.07 (2008), p. 124. [6] Mark Goodsell et al. “Naturally light hidden photons in LARGE volume string compactifications”. In: Journal of High Energy Physics 2009.11 (2009), p. 027. [7] Mark Goodsell, Saúl Ramos-Sánchez, and Andreas Ringwald. “Kinetic mixing of U (1) s in heterotic orbifolds”. In: Journal of High Energy Physics 2012.1 (2012), p. 21. [8] Nima Arkani-Hamed and Neal Weiner. “LHC signals for a superunified theory of dark matter”. In: Journal of High Energy Physics 2008.12 (2008), p. 104. [9] Pierre Fayet. “U-boson production in e+ e- annihilations, ψ and Y decays, and Light Dark Matter”. In: Physical Review D 75.11 (2007), p. 115017. [10] C Cheung. “C. Cheung, JT Ruderman, L.-T. Wang, and I. Yavin, Phys. Rev. D 80, 035008 (2009).” In: Phys. Rev. D 80 (2009), p. 035008. [11] David E Morrissey, David Poland, and Kathryn M Zurek. “Abelian hidden sectors at a GeV”. In: Journal of High Energy Physics 2009.07 (2009), p. 050. [12] Haipeng An et al. “Direct detection constraints on dark photon dark matter”. In: Physics Letters B 747 (2015), pp. 331–338. [13] Haipeng An et al. “Probing the dark sector with dark matter bound states”. In: Physical review letters 116.15 (2016), p. 151801. [14] FJ Rogers, HC Graboske Jr, and DJ Harwood. “Bound eigenstates of the static screened Coulomb potential”. In: Physical Review A 1.6 (1970), p. 1577. [15] R Agnese et al. “Search for low-mass weakly interacting massive particles using voltage-assisted calorimetric ionization detection in the SuperCDMS experiment”. In: Physical review letters 112.4 (2014), p. 041302. 112 [16] Daniel S Akerib et al. “First results from the LUX dark matter experiment at the Sanford Underground Research Facility”. In: Physical review letters 112.9 (2014), p. 091303. [17] Geoffrey T Bodwin, Eric Braaten, and G Peter Lepage. “P hys. Rev., D51, 1125 (1995),[Erratum”. In: Phys. Rev. D 55 (1997), p. 5853. [18] Andrea Petrelli et al. “NLO production and decay of quarkonium”. In: Nuclear Physics B 514.1-2 (1998), pp. 245–309. [19] Sea Agostinelli et al. “GEANT4—a simulation toolkit”. In: Nuclear instruments and methods in physics research section A: Accelerators, Spectrometers, Detectors and Associated Equipment 506.3 (2003), pp. 250–303. [20] Xinran He et al. “Practical lessons from predicting clicks on ads at facebook”. In: Proceedings of the Eighth International Workshop on Data Mining for Online Advertising. ACM. 2014, pp. 1–9. [21] Kenzo Nakamura, Particle Data Group, et al. “Review of particle physics”. In: Journal of Physics G: Nuclear and Particle Physics 37.7A (2010), p. 075021. [22] Yasmine Sara Amhis et al. “Averages of b-hadron, c-hadron, and τ-lepton properties as of 2018”. In: (2019). updated results and plots available at https://hflav.web.cern.ch/. arXiv: 1909.12524 [hep-ex]. [23] BaBar Collaboration et al. “Measurement of an excess of B → D(∗) τντ decays and implications for charged Higgs bosons”. In: Phys. Rev. D 88.072012 (2013), pp. 1303–0571. [24] Irinel Caprini, M Neubert, and LP Lellouch. “Dispersive Bounds on the Shape of B → D(∗) lν Form Factors”. In: Nucl. Phys. B 530.hep-ph/9712417 (1998), pp. 153–181. [25] Simone Bifani et al. “Review of Lepton Universality tests in B decays”. In: Journal of Physics G: Nuclear and Particle Physics 46.2 (2018), p. 023001. [26] Bernard Aubert et al. “Observation of the Semileptonic Decays B → D∗ τντ and Evidence for B → Dτντ ”. In: Physical review letters 100.2 (2008), p. 021801. [27] JP Lees et al. “Evidence for an excess of B̄ → D(∗) τ − ν̄τ decays”. In: Physical review letters 109.10 (2012), p. 101802. [28] Matthias Huschle et al. “Measurement of the branching ratio of B̄ → D(∗) τ − ν̄τ relative to B̄ → D(∗) l − ν̄l decays with hadronic tagging at Belle”. In: Physical Review D 92.7 (2015), p. 072014. [29] Roel Aaij et al. “Measurement of the ratio of branching fractions B( B̄0 → D∗+ τ − ν̄τ )/B( B̄0 → D∗+ µ− ν̄ µ )”. In: Phys. Rev. Lett. 115.LHCB-PAPER2015-025 (2015), p. 111803. [30] Giacomo Caria et al. “Measurement of R(D) and R(D∗ ) with a Semileptonic Tagging Method”. In: Physical review letters 124.16 (2020), p. 161803. 113 [31] A. Abdesselam et al. In: (2016). eprint: arXiv:1603.06711[hep-ex]. [32] S. Hirose et al. “Measurement of the τ Lepton Polarization and R(D∗ ) in the Decay B → D∗ τ − ν τ ”. In: Phys. Rev. Lett. 118 (21 2017), p. 211801. doi: 10.1103/PhysRevLett.118.211801. url: https://link.aps.org/ doi/10.1103/PhysRevLett.118.211801. [33] R. Aaij et al. “Measurement of the Ratio of the B0 → D∗− τ + ντ and B0 → D∗− µ+ ν µ Branching Fractions Using Three-Prong τ-Lepton Decays”. In: Phys. Rev. Lett. 120 (17 2018), p. 171802. doi: 10.1103/PhysRevLett. 120.171802. url: https://link.aps.org/doi/10.1103/PhysRevLett. 120.171802. [34] PARTICLE DATA GROUP Collab. K. Olive et al. In: Chin. Phys. C 38 (2014), p. 090001. [35] Daniel Shuteh Chao. “Measuring R(D(∗) ) for B → D(∗) τντ using Semileptonic Tags and Tau Decays to Hadrons”. PhD thesis. California Institute of Technology, 2018. [36] Jerome H Friedman. “Greedy function approximation: a gradient boosting machine”. In: Annals of statistics (2001), pp. 1189–1232. [37] Bernard W Silverman. Density estimation for statistics and data analysis. Vol. 26. CRC press, 1986. [38] Jon Louis Bentley. “Multidimensional binary search trees used for associative searching”. In: Communications of the ACM 18.9 (1975), pp. 509–517. [39] Martin S Andersen, Joachim Dahl, Lieven Vandenberghe, et al. “CVXOPT: A Python package for convex optimization”. In: abel. ee. ucla. edu/cvxopt 88 (2013). [40] Joseph L Doob. “The limiting distributions of certain statistics”. In: The Annals of Mathematical Statistics 6.3 (1935), pp. 160–169. [41] Bradley Efron. “Bootstrap methods: another look at the jackknife”. In: Breakthroughs in statistics. Springer, 1992, pp. 569–593. [42] Y Amhis et al. “Averages of b-hadron, c-hadron, and τ-lepton properties as of summer 2014”. In: arXiv preprint arXiv:1412.7515 (2014). [43] Adam K Leibovich et al. “Semileptonic B decays to excited charmed mesons”. In: Physical Review D 57.1 (1998), p. 308. [44] Florian Bernlochner. “Impact of D∗∗ on B → D∗ lν and R(D(∗) )”. In: Amplitude Analyses for Flavour Anomalies (Bristol 2019). url: https://indico. cern . ch / event / 810429 / contributions / 3376226 / attachments / 1875340/3087759/Amplitudes.pdf. [45] Florian U Bernlochner and Zoltan Ligeti. “Semileptonic B(s) decays to excited charmed mesons with e, µ, τ and searching for new physics with R(D∗∗ )”. In: Physical Review D 95.1 (2017), p. 014022. 114 [46] Vuosalo C. O., Telnov A. V., and K. T. Flood1. “Muon Identification Using Decision Trees”. In: BaBar Analysis Document (2010). url: https : / / babar.heprc.uvic.ca/BFROOT/www/Physics/BAD/vol18/01853. 003.pdf. [47] Alessandro Gaz et al. “Particle Identification Using Error Correcting Output Code Multiclass Classifier”. In: BaBar Analysis Document (2011). url: https : / / www . slac . stanford . edu / babar - internal / BAD / doc / detail.html?docNum=2199. [48] BaBar Collaboration et al. “Branching fractions measurement of τ − → K − nπ 0 ντ with n = 0, 1, 2, 3 and τ − → π − nπ 0 ντ with n = 3, 4”. In: BaBar Analysis Document 2318 (2012). [49] Yaroslav Ganin and Victor Lempitsky. “Unsupervised domain adaptation by backpropagation”. In: International conference on machine learning. PMLR. 2015, pp. 1180–1189. [50] Eric Tzeng et al. “Adversarial discriminative domain adaptation”. In: Proceedings of the IEEE conference on computer vision and pattern recognition. 2017, pp. 7167–7176. [51] Gilles Louppe, Michael Kagan, and Kyle Cranmer. “Learning to pivot with adversarial networks”. In: Advances in neural information processing systems 30 (2017). 115 APPENDIX A: DEEP ADVERSARIAL NETWORK FOR SYSTEMATIC UNCERTAINTY REDUCTION Monte Carlo simulation samples are almost always used in high energy physics to study the characteristics of particle decays and establish the measurement strategies. They are used to simulate the response of detectors in the event reconstruction stage, and to optimize criteria in the event selection stage. However, simulations are not always perfect, especially for background with high-multiplicity decay modes. The imperfection from simulations is usually considered as a source of systematic uncertainties in the measurements. For a simple case in which experimental data is categorized into only two labels, signal and background, we usually generate and simulate signal MC sample S mc and background MC sample Bmc . A typical and widely used strategy to establish event selection strategy is to optimize criteria using the combined samples of S mc and Bmc . To estimate corresponding systematics, people compare Bmc distribution with experimental data on signal suppressed regions, and evaluate the difference as systematics. This strategy is reliable but not optimal, as it does not take into account the information of MC/data discrepancy when optimizing selection criteria. The systematics due to MC/data discrepancy can be potentially reduced if the adversarial framework can be introduced. One architecture we played with is unsupervised domain adaptation [49], as shown in Fig. 8.1, while other adversarial architectures about unsupervised domain adaptation [50, 51] are also applicable. In this framework, the input variables x (normally physical variables of an event) are first transformed to f using G f (x; θ f ), and then used to train two classifiers together. The first classifier Gt (x; θ t ), similar to traditional multivariate analysis algorithms in event selection, aims to distinguish between signal and backgrounds. The second classifier G d (x; θ d ) aims to distinguish whether the input event is from MC simulation or experimental data. Each event has two labels, yt to label signal/background, and yd to label it is from MC simulation or experimental data. Therefore, a more robust classifier can be obtained by minimizing the combined loss function L = Lt − λL d Lt = Lt (Gt (G f (x)), yt ) (8.1) L d = L d (G d (G f (x)), yd ), where Lt and L d are the loss functions for signal/background classification and 116 MC/data classification. λ is a hyper-parameter to control the tradeoff between accuracy and robustness. Figure Figure 8.1: Architecture of unsupervised domain adaptation. The figure A.1: Diagram of the DANN algorithm. [31]is taken from [49]. for target domain points. The task classifier Gt is also a supervised learning algorithm, but the algorithm does not give significant improvement in our analysis, as uses (zs , ys ) as Although its training data instead. our analysis is dominated by statistical uncertainties, we believe this framework is We denote the domain loss and the prediction loss as sufficient general to reduce systematics when the mis-modeling of particle decays is not negligible. Ld (✓f , ✓d ) = Ld (Gd (Gf (x)), yd ), Lt (✓f , ✓t ) = Lt (Gt (Gf (xs )), ys ), (A.1) (A.2) where L can be any loss function such as the cross-entropy loss. Then, the overall loss function of the algorithm is L = Lt Ld , (A.3) where is a hyperparameter controlling the amount of adaptation. If = 0, it is as if we do not perform any feature transformation; the optimizer solely focuses on the task classifier performance. On the other hand, as becomes large, the domain adaptation will be the main focus, disregarding the task classifier performance. The optimal parameters are the saddle point solution to the following optimization problem: min max L(xs , xt , ys ; ✓f , ✓d , ✓t ). (A.4) ✓f ,✓t ✓d Figure A.1 shows the diagram of the DANN algorithm. A.3 Reverse validation Given that the test set has no labels, how can we quantify the performance of the task classifier? This is a problem that the vanilla supervised learning algorithms also face given