Open Collections

UBC Undergraduate Research

Applications of the Wavelet Transform to B Mixing Analysis Cadien, Adam Samuel 2008-06-30

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
52966-Cadien_Adam_UBC_2008_Physics_449_Honours_Thesis.pdf [ 1.48MB ]
Metadata
JSON: 52966-1.0085956.json
JSON-LD: 52966-1.0085956-ld.json
RDF/XML (Pretty): 52966-1.0085956-rdf.xml
RDF/JSON: 52966-1.0085956-rdf.json
Turtle: 52966-1.0085956-turtle.txt
N-Triples: 52966-1.0085956-rdf-ntriples.txt
Original Record: 52966-1.0085956-source.json
Full Text
52966-1.0085956-fulltext.txt
Citation
52966-1.0085956.ris

Full Text

Applications of the Wavelet Transformto B Mixing AnalysisbyAdam Samuel CadienA THESIS SUBMITTED IN PARTIAL FULFILMENT OFTHE REQUIREMENTS FOR THE DEGREE OFBachelor of Science (Hons)inThe Faculty of Science(Physics)The University Of British Columbia(Vancouver, Canada)April, 2008c Adam Samuel Cadien 2008AbstractThe neutral B mesons B0 and B0s can under go flavor changing oscillations due to interactionsby the weak force. Experiments which measure the frequency of these state transitions produceextremely noisy results that are difficult to analyse. A method for extracting the frequency of Bmesons oscillations using the continuous wavelet transform is developed here.In this paper the physics of B meson mixing is related, leading to the derivation of a func-tion describing the expected amount of mixing present in B0 and B0s meson decays. This resultis then used to develop a new method for analysing the underlying frequency of oscillation inB mixing. An introduction to wavelet theory is provided in addition to details on interpretingdaughter wavelet coefficient diagrams. Finally, the effectiveness of the analysis technique pro-duced, referred to as the Template Fitting Method, is investigated through an application todata generated using Monte Carlo methods.Adam Cadienadamcadien@gmail.comiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 B- ? Oscillations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.1 CP Violation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.1.1 K Mixing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2 B Mixing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2.1 Quantum Treatment of the B ? System . . . . . . . . . . . . . . . . . . . 42.2.2 The Probabilities of Mixed and Un-mixed Decay . . . . . . . . . . . . . . . 52.2.3 Sources of Error in Detection . . . . . . . . . . . . . . . . . . . . . . . . . . 72.3 The CKM Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.3.1 Finding CKM matrix elements Vtd and Vts . . . . . . . . . . . . . . . . . . 83 Applications of the Wavelet Transform in Data Analysis . . . . . . . . . . . . . 103.1 Mother Wavelets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103.2 The Continuous Wavelet Transform . . . . . . . . . . . . . . . . . . . . . . . . . . 103.2.1 Wavelet Coefficient Diagrams . . . . . . . . . . . . . . . . . . . . . . . . . 114 Finding deltam with the Continuous Wavelet Transform . . . . . . . . . . . . . . . 154.1 Data Generation and Error Simulation . . . . . . . . . . . . . . . . . . . . . . . . 154.1.1 Selecting the Mother Wavelet . . . . . . . . . . . . . . . . . . . . . . . . . 154.2 The Template Fitting Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184.2.1 The Bounding Box . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184.2.2 Calculating chi for deltam Extraction . . . . . . . . . . . . . . . . . . . . . . . 194.3 Results & Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26iiiTable of ContentsAppendicesA Programs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27A.1 Program to Generate Pu,m(t) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27A.2 Program that Applies Simulated Detection Errors . . . . . . . . . . . . . . . . . . 29A.3 Program to Extract deltam using the Template Fitting Method . . . . . . . . . . . . 32ivList of Figures2.1 Lowest Order B Mixing Diagrams . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 The Unitarity Triangle with deltam constraints. . . . . . . . . . . . . . . . . . . . . . 93.1 The Morlet and Meyer mother wavelets. . . . . . . . . . . . . . . . . . . . . . . . . 113.2 Applying the Continuous Wavelet Transform. . . . . . . . . . . . . . . . . . . . . . 133.3 Coping with noise in the continuous wavelet transform. . . . . . . . . . . . . . . . 144.1 Simulated B decay data using the function Pu,m(t). . . . . . . . . . . . . . . . . . . 164.2 Simulated B decay data with detection errors. . . . . . . . . . . . . . . . . . . . . . 164.3 Coefficient diagrams using the Morlet and Meyer mother wavelets. . . . . . . . . . 174.4 The decomposing the Pu,m(t) coefficient diagrams. . . . . . . . . . . . . . . . . . . 194.5 An example of how to set the bounding box on a coefficient diagram. . . . . . . . . 204.6 Setting the Bounding box on test data. . . . . . . . . . . . . . . . . . . . . . . . . . 224.7 Results: First Iteration of the Template Fitting Method . . . . . . . . . . . . . . . 234.8 Results: Second Iteration of the Template Fitting Method. . . . . . . . . . . . . . 24vAcknowledgementsColin Gay, under whose supervision I have learned what it means to be a researcher. In additionto his support of my undertaking this project, Colin has been an essential source of advice and aguide to all of my work.Alex Muir for always showing interest in my work, and doing his best to answer my constantbarrage of questions regarding everything from the inner workings of B detectors to the propermethods of statistical analysis.Cheng-Wei Loh and Bill Mills who proof-read all of the derivations presented in Chapter 2 andwho both presented me with great advice on how to write this paper.Bryce Baril for proof reading the entire paper.Everyone at UBC?s ATLAS TRT lab, who provided me with a friendly and productive work-ing environment.viChapter 1IntroductionParticle physics appears to characterize our world through symmetries. When evidence for some-thing as elegant as Charge Parity (CP) symmetry is found it is typically embraced by the scientificcommunity at large. So in 1964 when Val Fitch and James Cronin discovered proof of CP vi-olation in the decay of the K mesons it had profound effects. The effects of CP violation aredescribed experimentally by approximately 30 quantities. Accurately measuring such values fu-els the creation of many particle accelerators, and is an active area of research for thousands ofphysicists today.The Cabbibo-Kobayashi-Maskawa (CKM) matrix is one way of representing the amount ofCP violation that exists in quark to quark interactions. The neutral B meson1 exhibits CPviolation when it oscillates between its particle-antiparticle states which is used to bound three ofthe CKM matrix values, the magnitude of each element is proportional to the oscillation frequencyof the B which is determined by the mass difference (deltam) between the two states. The massdifference between the B and ? is the largest uncertainty currently present in the CKM Matrix.Specifically the measurement of deltams is very difficult to obtain since the Bs oscillates sev-eral times during its lifetime. An experiment involving B mixing measurements may measureupwards of a half-million events [1] before seeing B particle/anti-particle transitions. The initialmeasurement of deltams took nearly 20 years to gather sufficient statistics to be able to extract thefrequency of oscillation. Precision on the current measurement of deltams = 17.32+0.35-0.18 ? 0.07ps-1[2] needs to be greatly increased in order to verify the unitarity of the CKM matrix. This is thecurrent method used by experimentalists for verifying the Standard Model, which may also leadto the discovery of new particles or interactions.The extraction of deltam from measurements of the B?s lifetime is complicated because of itshigh decay rate and unusually high frequency of oscillation. Upwards of 30% of the measurementstaken are due to errors in the detection process or background produced in a collision. Thetypical methods of frequency extraction, Fourier and Discrete Wavelet analysis, are ineffectiveat identifying the oscillation frequency to a high degree of certainty. Fourier methods are mosteffective when transforming a signal which has a large number of periods. The dataset analyzedcontains only a few visible periods, meaning that when transformed errors are introduced. Similarproblems arrise when applying the discrete wavelet transform which also has difficulties workingwith finite length signals.The Template Fitting Method is introduced here, a new data analysis technique for extractingthe frequency of curves with low statistics and a low signal to noise ratio which is directlyapplicable to the extraction of deltam. The continuous wavelet transform is used to decomposea signal into its frequency content and a fitting process which uses predefined templates then1The symbol B will be used when results pertaining to both the B0 and B0s are discussed, the same applies to?B. Similarly deltam refers to both deltams and deltamd.1Chapter 1. Introductionextracts deltam. Using Monte Carlo data it is demonstrated that the Template Fitting Method canextract the mixing frequency to the same degree of accuracy using half of the statistics used incurrent analysis methods.This document is split primarily into 3 parts. Chapter 2 provides the background and theoryof B mixing. While in Chapter 3 the process behind and effects of the wavelet transform arediscussed. The Template Fitting Method along with results demonstrating its effectiveness arepresented in Chapter 4. Finally a conclusion summarizing the information presented is given inChapter 5.2Chapter 2B- ? OscillationsThe instability of some neutral mesons due to the weak force causes them to undergo a changeof state into their antiparticle. This rare event is the only source of CP violation that has beenexperimentally observed. The B meson?s rate of oscillation is of interest because it is a measurablequantity of the amount of CP violation that takes place in some quark to quark interactions.This chapter begins with a background on CP Violation, whose effect on K mesons is givenin section 2.1.1. The K experiments are used as a basis for providing details of B mixing. Ananalysis of the B- ? system is provided in section 2.2.1, which leads into the derivation of theprobability of mixing presented in section 2.2.2. Details on the errors introduced in B detectionare given in section 2.2.3. Finally how these results are used to describe the likeliness of CPViolation in quark interactions is in section 2.3.2.1 CP ViolationCP is the combined operation of charge and parity conjugation. C is the charge conjugationoperator, which represents the exchange of a particle with its anti-particle. P, parity conjugation,is a simultaneous reflection about all axes.Before 1955 it was thought that all particles obey symmetry through parity and charge indi-vidually. When Gell-Mann and Pais pointed out that interactions between the K0 and ?0 occurthrough a weak interaction[3] an investigation into the properties of K mesons ensued. Shortlyafter in 1956 it was demonstrated that parity symmetry is violated by the K+[4]. So the theoryof the time changed so that symmetry should truly be conserved under the combination of chargeand parity (CP symmetry). However even this to was found to be broken in 1964 when CPviolation was found to occur in of decaying neutral K mesons, which due to weak interactionsinitially discovered by Gell-Mann and Pais, can mix in a unique way [5].2.1.1 K MixingUnder CP symmetry the K meson has the following normalized CP eigenstates|K1angbracketright = 12(|K0angbracketright -| ?0angbracketright) and |K2angbracketright = 12(|K0angbracketright + | ?0angbracketright) (2.1)The K1, with a CP of +1, decays into two pions while the K2, with a CP of -1, decays intothree pions, which conserves CP symmetry. Accordingly the K1 when decaying purely in thetwo-pion configuration has a much shorter lifetime than the K2 decaying purely in its long livingthree-pion configuration. The K meson doesn?t always decay into its K1 or K2 states. It hasbeen observed that a small fraction (represented by epsilon1) of the long living KL decay as K1 and3Chapter 2. B- ? Oscillationsconversely some of the short living KS decay as K2. When the K decays, it decays not as aneigenstate of the CP but as a combination of these states:|KLangbracketright = 11 + |epsilon1|2 (epsilon1|K1angbracketright + |K2angbracketright) and |KSangbracketright = 11 + |epsilon1|2 (|K1angbracketright + epsilon1|K2angbracketright) (2.2)Similarly CP violation is seen in the particle-antiparticle mixing which also occurs2 in B0 andB0s.The lifetimes of the K1 and K2 states differ by approximately[6] deltatauK2-1 = 5.2 ? 104ps.Consequently experimentally detecting whether a K has decayed as a K1 or K2 is simplified.A beam of K1?s will decay over a few centimeters and K2?s over several meters thus making itvery easy to detect an event with CP violation. In contrast the Bs meson has a mean lifetimeof tauB = 1.5ps and a lifetime difference on the order of deltatauB proportional 1e-2ps resulting in a difference ofdecay distances on the order of ?m?s.2.2 B Mixing2.2.1 Quantum Treatment of the B ? SystemMuch like the K0- ?0 system, the Hamiltonian of the B- ? system under the strong force is simplyH0 =parenleftBiggm0 00 m0parenrightBigg. (2.3)B mesons are made up of a b quark and d or s quark with flavor eigenstates (where g = (d,s) forthe B0d and B0s respectively):|Bangbracketright = |? angbracketright and | ?angbracketright = |b?angbracketright (2.4)Now taking into account weak interactions results in a Hamiltonian which is the sum of H0 withterms due to an energy shift of the particle-antiparticle system and a set of new states accessiblethrough the coupling of |Bangbracketright and | ?angbracketright (defined by the element W12). It has the form3:H =parenleftBiggM M12Masteriskmath12 MparenrightBigg- i2parenleftBiggG G12Gasteriskmath12 GparenrightBigg(2.5)where M = m0 +deltaE is the mass eigenstate of the strong force summed with the mass shift due tothe weak force. The off diagonal elements M12 = W12 + deltaE12 and G12 is the mass and decay rateof the coupled |Bangbracketright and | ?angbracketright states. Its these off diagonal elements which explain the energy shiftduring transition between B and ?. This is shown in figure 2.1 as the exchange in W bosons asexpected from the Standard Model.The B?s mass eigenstates are not states with definite flavor. These mixed eigenstates areconventionally named light and heavy in parallel to the K?s short and long states. The eigenstatesare|BH,Langbracketright = p|Bangbracketright minusplusq| ?angbracketright, (2.6)2It is expected that the D0 meson will also exhibit this phenomena in measurements to be taken at the LHC.3A detailed derivation of this Hamiltonian is available in [7].4Chapter 2. B- ? OscillationsFigure 2.1: The Feynman diagrams describing the coupling of the B0 and ?0 particles. The B0sand ?0s interact in a similar manner.withqp = -radicaltpradicalvertexradicalvertexradicalbtM12 - i2G12Masteriskmath12 - i2Gasteriskmath12 , (2.7)and with eigenvalueslambdaH,L = M - i2G ?radicalBigg(M12 - i2G12)(Masteriskmath12 - i2Gasteriskmath12). (2.8)Knowing that H = M - i2G and then diagonalizing H revealsHdiag =parenleftBigglambdaH 00 lambdaLparenrightBigg= Mdiag - i2Gdiag, (2.9)which means thatMdiag = Re(Hdiag ) and Gdiag = -2Im(Hdiag ). (2.10)We can now definedeltam = Re(lambdaH) -Re(lambdaL) and deltaG = -2Im(lambdaH) + 2Im(lambdaL). (2.11)These values are needed to describe the wave function of particles, so that we can see how theyevolve with time through Schrodinger?s Equation.2.2.2 The Probabilities of Mixed and Un-mixed DecayWe are now prepared to derive the probabilities of a B/ ? decaying as its anti-particle, referred toas the probability of mixing. The probability of mixed/un-mixed decay is used for data generationpurposes, presented in section 4.1, since it matches experimental results for the probability ofmixed and un-mixed B decays.We begin with the time independent states of B and ?:|Bangbracketright = 12p|BHangbracketright + 12p|BLangbracketright and | ?angbracketright = - 12q|BHangbracketright + 12q|BLangbracketright (2.12)5Chapter 2. B- ? Oscillationswhich were found by inverting equation 2.6 and substituting in |q|2 + |p|2 = 1. The time depen-dence for BH,L is thenphiH,L(t) = e-ilambdaH,Lt = e-iMteG2 te?deltaG4 te?ideltam2 t, (2.13)which is tacked onto the stationary state to obtain the time dependent wave equations|B(t)angbracketright = 12e-iMteG2 t[parenleftBigedeltaG4 teideltam2 t + e-deltaG4 te-ideltam2 tparenrightBig|Bangbracketright+qpparenleftBigedeltaG4 teideltam2 t -e-deltaG4 te-ideltam2 tparenrightBig| ?angbracketright]| ?(t)angbracketright = 12e-iMteG2 t[pqparenleftBigedeltaG4 teideltam2 t -e-deltaG4 te-ideltam2 tparenrightBig|Bangbracketright+parenleftBigedeltaG4 teideltam2 t + e-deltaG4 te-ideltam2 tparenrightBig| ?angbracketright].(2.14)The normalization of these states is theneta2 =integraldisplay infinity0angbracketleftB(t)|B(t)angbracketright dt = G2bracketleftBigg1 + |q/p|2G2 -deltaG2/4 +1 -|q/p|2G2 + deltam2bracketrightBigg?2 =integraldisplay infinity0angbracketleft ?(t)| ?(t)angbracketright dt = G2vextendsinglevextendsinglevextendsinglevextendsingle qpvextendsinglevextendsinglevextendsinglevextendsingle2bracketleftBigg1 + |q/p|2G2 -deltaG2/4 +1 -|q/p|2G2 + deltam2bracketrightBigg. (2.15)The likeliness of a particle-antiparticle transition is described byPm(B arrowdblboth ?)(t) = 1?eta2 |angbracketleft ?|B(t)angbracketright|2 = 12?eta2 e-Gt(coshdeltaGt2 -cosdeltamt) andPm( ? arrowdblbothB)(t) = 1eta2 |angbracketleftB| ?(t)angbracketright|2 = 12eta2 e-Gt(coshdeltaGt2 -cosdeltamt), (2.16)since the cross-terms cancel out in both cases these are greatly simplified. The correspondingprobability that a transition does not occur and the particle is then un-mixed isPu(B arrowdblbothB)(t) = 1eta2 |angbracketleftB|B(t)angbracketright|2 = 12eta2 e-Gt(coshdeltaGt2 + cosdeltamt) andPu( ? arrowdblboth ?)(t) = 1?eta2 |angbracketleft ?| ?(t)angbracketright|2 = 12?etae-Gt(coshdeltaGt2 + cosdeltamt). (2.17)The current measured values[2] of the ratio of decay difference (between the particle-antiparticlestates) to mean decay rate isdeltaGB0sGB0s approxequal 0.31 ,deltaGB0GB0 approxequal 0.08 (2.18)andx0s = deltamB0sGB0s> 19.9 , x0 = deltamB0GB0> 0.78. (2.19)It follows that deltaG lessmuch deltam for both B mesons [8]. Therefore G12 lessmuch M12 and from equation 2.7|q| = |p| thus neglecting CP violation in B mixing, this turns out to be a good approximation.In this limit deltaG = 0, simplifying the eta term such that eta = ? and the mixed and un-mixed decay6Chapter 2. B- ? Oscillationsprobabilities become equal for the B and ? and equations 2.16, 2.17 becomePu,m(t) = 12Ge-Gt(1 ?cosdeltamt). (2.20)The lifetime, G, of a particular B or ? is typically measured to be on the order of G = 1.5ps.The value of interest is however the oscillation frequency, deltam. It is the mass difference betweenthe heavy and light states of the B which is the frequency of oscillation between the matterarrowdblbothanti-matter states of the B meson. As discussed in section 2.3.1 the value deltam is relevant to determiningthe CKM elements Vtd and Vts.2.2.3 Sources of Error in DetectionWhen a collision occurs inside of a detector, given that the proper particles4 are used and if thereis enough energy between them, there is a small possibility that B or ? mesons may be produced.For every collision that occurs within the detector?s range the position and energy of anyemitted particles from the collision are tracked. Should a B or ? be produced, reconstruction ofthe event using emission data produces the presumed points of production and decay of the meson.Similar reconstruction methods are used to find the momentum of the particle and its decayproducts. Knowing the production and decay points gives the total distance, L, traveled by themeson before decaying. It is then possible to calculate the lifetime of the particle: t = LMB/pBwhere MB is the mass of the meson.Faster moving particles produce tracks in the detector which require a position resolutiongreater than what most detectors are able to provide. This decreases the accuracy in the positionmeasurement of the particle, an effect known as momentum smearing.Errors are also introduced through false detection of B mesons. This background error iscaused by misidentification of the decay products of the B. Incorrect tags can occur from im-properly detecting the lepton present in the decay resulting in detection of the opposite chargethat is actually present. Depending on the method of detection and particle of interest, fake rateis typically responsible for 5% of B detections.Determining whether the particle detected is a B or ? is typically done by flavor tagging theb quark at production and decay points. Of course the B will be detected using whatever methodis designed by the experiment. The SLD group uses the angle between the b quark and incomingelectron to tag the flavor of the quark, with an accuracy of 66%[7]. At BaBar and Belle thecombined tagging efficiency has been estimated to be 28% [10], this is the percentage of swappedevents used in generating simulation data.2.3 The CKM MatrixThe CKM matrix holds experimentally found values which describe the amplitude of flavour-changing weak interactions. Each element is defined as Vij, where i is one of u,c, or t and j isd,s, or b, and the value |Vij|2 is the probability for a j-quark to transition to a i-quark. Applying4Experimentally B mesons are most commonly produced from an electron positron collision, although CDF usesa proton anti-proton collision7Chapter 2. B- ? Oscillationsthe unitary constraint to the matrix, meaningVCKMV ?CKM = 1, (2.21)requires that the probabilities of quark transitions be normalizable. This stands to reason sinceeach value squared is a probability, the sum of any row or column should be 1.There are numerous ways of representing the CKM matrix, most common is Wolfenstein?sparameterisation [9] shown here to O(lambda3):VCKM =parenlefttpparenleftexparenleftbt Vud Vus VubVcd Vcs VcbVtd Vts Vtbparenrighttpparenrightexparenrightbt =parenlefttpparenleftexparenleftbt 1 -12lambda2 lambda Alambda3(rho -ieta)-lambda 1 - 12lambda2 Alambda2Alambda3(1 -rho -ieta) -Alambda2 1parenrighttpparenrightexparenrightbt (2.22)where lambda is a function of the Cabbibo angle (lambda = sin(thetac)), and A,rho and eta are real parameters.Specifically of interest here is eta which defines the CP violating phase present in the flavor changingtransitions. Imposing unitarity of VCKM, the CKM matrix is constrained by 6 conditions. Oneof these normalization constraints,VudV asteriskmathub = -VcdV asteriskmathcb -VtdV asteriskmathtb, (2.23)is directly related to observables found from B meson decay. Each of the 6 equations arising from2.21 are of the form of equation 2.23 and each makes a Unitarity Triangle. The magnitude ofeach term from equation 2.23 is the length of one side of a Unitarity Triangle which exists in thecomplex plane.The area of each Unitarity Triangle is the same. Existing in the complex plane means thatonly terms with a complex phase contribute to the area of a triangle. Thus only the elements Vtdand Vub, which have a complex phase due to CP Violation, contribute to the area.The Unitarity Triangle?s area is a measure of the magnitude of CP Violation which occursin a quark-quark transition. It also allows for a visual representation of the constraints placedon the CKM matrix elements from experimental results. The experimental constraints on thetriangle (described by equation 2.23) are shown in figure 2.2. The uncertainty of the measuredvalue of deltams and deltamd directly correlate to errors in the lengths and angles of the triangle. Bydecreasing the uncertainty of an experimental result it is possible to put tighter constraints onthe Unitarity Triangle which may show that the CKM matrix does not actually meet the unitarycondition (equation 2.21, in this way the standard model is verified. Should the CKM elementsnot meet the unitary condition this would violate the current structure of the standard model.2.3.1 Finding CKM matrix elements Vtd and VtsThe B oscillation frequency has been shown to be the mass difference between the B and itsanti-particle ?.Omitting lengthy calculations beyond the scope of this paper the end result isdeltams = KsmBs|V asteriskmathtsVtb|2, (2.24)deltamd = KdmBd|V asteriskmathtdVtb|2, (2.25)8Chapter 2. B- ? OscillationsFigure 2.2: The Unitarity Triangle plotted in the complex plane. This is one of the two UnitarityTriangle whose sides all have lengths that are of the same order and are sensible to representgraphically. The dashed line represents the bounds placed on the error of the current placementof the vertex with angle alpha. The experimentally measured values of deltams and deltamd are responsiblefor setting the bounds on the CKM matrix elements Vtd and Vts.where the values mBs and mBd are the mass of the Bs and Bd mesons respectively. The valuesKs and Kd are functions of the bag parameter and the B decay constant [7], theoretically derivedvalues which will be treated as constants here. Additionally the ratio of |Vts/Vtd|2 can be relatedwhen approximating the diagonal element |Vtb| approxequal 1. So from equations 2.24 and 2.25deltamsdeltamd =KsmBsKdmBdvextendsinglevextendsinglevextendsinglevextendsingle VtsVtdvextendsinglevextendsinglevextendsinglevextendsingle2 . (2.26)This is the relationship between the B mixing oscillation frequency deltam and the CKM matrixelements Vtd.9Chapter 3Applications of the WaveletTransform in Data AnalysisAlthough wavelet theory has been a topic of study since the 1930?s, their application to dataanalysis began only 20 years ago. When Grossman and Morlet proposed the wavelet transform[11] it was used as an improvement to windowed Fourier analysis in searching for oil using sonarmethods. Because it provides a sturdy platform to carry out frequency analysis of previously un-usable datasets (due to a high noise ratio, discontinuities, or low statistics) the wavelet transformhas become a popular technique for data analysis.3.1 Mother WaveletsThe wavelet transform involves the use of a master wavelet function, referred to as the motherwavelet, and some function or dataset of interest. Mother wavelets must exist in the HilbertSpace for square integrable functions, that isintegraldisplay infinity-infinity|psi(t)|2dt < infinity. (3.1)They must also be differentiable, this ensures the wavelet will not have any discontinuities whenevaluated at any scale or position (a,b,t). Any function that satisfies these constraints can beused for the continuous wavelet transform. Several wavelet functions arise that satisfy theseconstraints, the mother wavelets shown in figure 3.1 are the primary candidates for use in thedata analysis section and are discussed more there.3.2 The Continuous Wavelet TransformThe only constraint on a the function to be transformed f(x) is that it must occupy the samefunctional space as the wavelet psi. The wavelet transform of f(x) is defined asWpsi = angbracketleftpsi,f(x)angbracketright, (3.2)where the inner product isangbracketleftf,gangbracketright =integraldisplay infinity-infinityf(x)asteriskmathg(x)dx. (3.3)The transformed function Wpsi is called a daughter wavelet. The daughter wavelet is meaninglesshowever, without a specified scale and position of the mother wavelet. These scale and translationparameters are labeled a and b respectively, and the daughter wavelet is a function of these10Chapter 3. Applications of the Wavelet Transform in Data Analysis(a) (b)Figure 3.1: (a): The Morlet mother wavelet. (b): The Meyer mother wavelet, the accompanyingscaling function is not shown.Wpsi(a,b). A wavelet coefficient CW (a,b) is then definedCW (a,b) =integraldisplay infinity-infinityWpsi(a,b)dx. (3.4)The continuous wavelet transform results in a set of wavelet coefficients corresponding to eachdaughter wavelet within the bounds x0 <=< b <=< xlscript and 0 <=< a <=< hlscript. The wavelet coefficient is ameasure of how well the mother wavelet fits the function at a certain scale and position.3.2.1 Wavelet Coefficient DiagramsEvaluating a 1-D function using the wavelet transform results in a 2-D description of the functionnow as a function of scale and position. Larger scale values correspond to a wider analysingfunction, that is psia(x) = psi(x/a), and lower frequencies are detected. These dilated waveletshave no dependence on the high frequency content of the function being analysed. The oppositeis true for small scale values, any low frequency content is lost when the transform is taken. Sothe daughter wavelet?s coefficient is maximized when it?s central frequency is near some frequencycontent of the dataset and the wavelet is in phase with these oscillations. These are the primaryconsiderations taken into account when evaluating the 2 dimensional wavelet coefficient diagram.An example of finding the wavelet coefficients for a continuous wavelet transform is givenin figure 3.2. Shown in figure 3.2(a), are the initial wavelet (dashed) with its parameters with[a,b] = [1,0] and the sinusoid (solid). The wavelet is then shifted through all of the combinationsof translation and scale finding the coefficient CW (a,b). In figure 3.2(b) the mother wavelethas parameters [a,b] = [5,0]. This appears to be the central frequency of the sinusoid, whichthen becomes clear in figure 3.2(c) when the wavelet is drawn with parameters [a,b] = [5,pi/2].When the wavelet is transitioned to be in phase with the function and matches its frequency,the largest values of CW (a,b) are found. With the wavelet completely out of phase, as in figure11Chapter 3. Applications of the Wavelet Transform in Data Analysis3.2(d) where [a,b] = [1,pi], the smallest coefficients are produced. As the scale of the waveletis increased, coefficients will never be as greater than when it was perfectly fit at a = 5. Thewavelet coefficient diagram is then shown in figure 3.2(f), which is made by plotting all of thecoefficients as a function of their scale and position. The light areas represent large coefficients,black areas are smaller values of the daughter wavelet coefficient and they represent a poor fit.Peaks exist when the wavelet has a scale near a = 5 and when it is in phase with the sinusoid.The coefficient diagram alone makes it possible to interpret much about a dataset. Periodiccontent which may have been hidden in noise or background becomes apparent when frequencymatching occurs. This property is especially useful when there are multiple periodic functions in asingle dataset, as shown in figure 3.3. Clearly the continuous wavelet transform is still applicablefor datasets with low statistics and high noise content.12Chapter 3. Applications of the Wavelet Transform in Data Analysis(a) (b)(c) (d)(e) (f)Figure 3.2: (a)-(e) The mother (dashed line) and a sinusoid function (solid line). (f) The contin-uous wavelet transform coefficient diagram. Light areas represent large coefficients, black areasare small coefficients CW (a,b). 13Chapter 3. Applications of the Wavelet Transform in Data Analysis(a) (b)(c) (d)Figure 3.3: (a) The function f(t) = cos(t/2)sin(5t)sin(10t). (b) The continuous wavelet trans-form of f(t), each component of the function is visible in the diagram. The low frequency cosineterm is present in the broad scale, its inflection point is visible at t = pi. The higher frequencyterm sin(5t) exists over 90 < a < 240,the slope inwards due to interactions with the cosine term.The highest frequency term sin(10t) exists in the small peaks at 25 < a < 65. (c) The functionf(t) = cos(t/2)sin(5t)sin(10t) with each point perturbed by a random Gaussian number withwidth sigma = 0.1. (d) The coefficient diagram of the data from (c). Although not as smooth asdiagram (b) all frequency content in the transform of the noisy function is still clearly visible.The function f(t) could easily be extracted from this.14Chapter 4Finding deltam with the ContinuousWavelet TransformThe data analysis method for extracting the oscillation frequencies of mixed and un-mixed Bmesons is described in this chapter. A brief overview of data generation methods is given first,followed by a discussion on how a mother wavelet was selected in section 4.1.1. In section 4.2 thedata analysis method developed using the wavelet coefficient diagrams is presented 5.4.1 Data Generation and Error SimulationExperimental datasets for B meson oscillations are small and too noisy to work with while devel-oping a frequency extraction method. However an extremely accurate model for the probabilityof decay as a mixed or un-mixed B, equation 2.20 presents the opportunity for an alternativemethod for algorithm development. For a given frequency deltam a large master set of data is gen-erated using Monte Carlo methods and the master function Pu,m(t), shown in figure 4.1. This isalso performed by the Matlab script given in appendix sections A.1.To ensure that an algorithm developed on test data will work equally as well on experimentaldata it is necessary to mimic the errors involved in detection. Background detection errors areintroduced by sampling multiple curves which dilute the dataset with false B detections. Thelargest source of errors, the swapping of B and ? events, is achieved in the test data by simplyswitching a number of mixed points for un-mixed. The finite resolution of a particular detectoris mimicked by disturbing every sample point by a Gaussian random number centered at themeasured time of decay. Additionally, to account for the smearing affect caused by the widerange of time dilation experienced by the particles, the width of the Gaussian used to shift eachsample is increased as a function of the decay time. An implementation of error simulation andinjection is provided in the appendix section A.2. The data presented in figure 4.1 convolutedwith the simulation errors described above is shown in figure 4.2.4.1.1 Selecting the Mother WaveletChoosing which mother wavelet is appropriate for an analysis is dependent upon the dataset tobe analyzed and the type of analysis to be performed. A wavelet which is able to extract theunderlying function in a dataset is the primary goal in choosing the mother wavelet. For thepurpose of frequency extraction from a dataset with a large background ratio, selecting a waveletwhich closely matches the underlying function will increase the likeliness of a good decomposition.A mother wavelet which is close in shape to the underlying curve will highlight its broad structure5The code required for obtaining these results is given in the Appendices.15Chapter 4. Finding deltam with the Continuous Wavelet Transform(a) (b)Figure 4.1: (a) Histogram of the simulated probability of decay for un-mixed B mesons, Pu(t).(b) The corresponding histogram of the probability of decay as a mixed particle Pm(t). Bothplots contain 50,000 points divided into 1,000 evenly spaced bins. These samples were generatedusing Monte Carlo methods and have no simulated errors injected.(a) (b)Figure 4.2: (a) The Pu(t) data with detection errors. (b) The corresponding Pm(t) curve withthe similar errors. Both datasets contain 30% swapped events, which is introduced by switchingPm(t) samples with Pu(t) samples. Additionally 5% of the points are replaced by samples fromrandomly generated functions, this simulates background errors. Finally all of the points areshifted by a gaussian random number with a varying width sigma = 0.1t, this is a function of thetime of their decay and mimics the effect of smearing described in section 2.2.3. Both of theseplots were made using the program in appendix section A.2.16Chapter 4. Finding deltam with the Continuous Wavelet TransformFigure 4.3: (a) Coefficient diagram made using the Morlet wavelet. (b) Coefficient diagram madeusing the Meyer wavelet. Both are transforms of the noise injected Pmt data given in figure4.2(b).resulting in clear sections of periodic behaviour in the wavelet coefficient diagram. The oscillationspresent in the coefficient diagram are caused by the mother wavelet coming in and out of phasewith the function being transformed.There are some subtle qualities of mother wavelets which may not become clear for a certaindataset until an inspection of the coefficient diagram is done. In figure 4.3 the wavelet transformfor the sample dataset from figure 4.2, using both the Morlet and Meyer mother wavelets, isgiven.The Morlet wavelet proves to be effective in providing a stable series of oscillations in thefrequency regime of interest. The Morlet wavelet is defined by the functionpsi(x) = Ce-x2/2cos(5x) (4.1)where C is the normalization constant. This clearly bears a resemblance to the mixed and un-mixed decay probability (function 2.20) as shown in figure 3.1(a).A second mother wavelet candidate, the Meyer wavelet, is ineffective at high swapping rates,see figure 4.3(b). When more then similar30% of the points become swapped the mixing amplitude isnot visible to the eye. For the transform to be effective with such a large amount of background, awavelet which very closely matches the function to be extracted is needed. Additional propertiessuch as exisiting for several periods at a significant amplitude also makes the Morlet wavelet moresensitive local periodic behaviour.17Chapter 4. Finding deltam with the Continuous Wavelet Transform4.2 The Template Fitting MethodTo begin with, a template is defined here as the array of wavelet coefficients for a particulardaughter wavelet. The quality of interest in these templates is the frequency of the master curveused in data generation. The extraction of deltam is performed following the procedure referred toas the Template Fitting Method:1. Take the Continuous Wavelet Transform of the dataset to be analysed.2. Define a bounding box which isolates the frequency peaks of the deltam oscillations.3. Generate a set of templates over a desired frequency range.4. Calculate the chi value for each template against the experimental data template, within thebounding box.5. Iterate over steps 3 and 4 reducing the step size of the frequency and the bounds as necessary,concentrating on the frequencies near the minimum chi2 value.6. The oscillation frequency is the frequency of the template with the minimum chi value.The underlying concept is that signals with like frequencies will result in similar peak distri-butions in the same scale region of the coefficient diagrams. We can generate datasets which veryclosely mimic the experimentally sampled one, this is taken advantage of by generating templatesat varying frequencies. By taking a point by point difference between the wavelet coefficientsinside of a predetermined bounding box, exactly how well the generated data matches the exper-imental data is found. The template with the closest frequency to the experimental data shouldthen have the smallest chi value.Data generation for the templates can take a significant amount of time6 , sometimes up toan hour for each template. When a high degree of accuracy is desired it has been shown thatperforming Template Fitting Method in a recursive manner reduces the number of templatesneeded. With each iteration the bounds for the frequencies of the next iteration should be set tobe similar?10% of the current frequency with the minimum chi value. As the minimum and maximumfrequency approach a common value, the stepsize deltaftmp of the template frequency is reduced toprovide a sufficient number of templates to compare with, typically on the order of 10 templatesper iteration. All results shown in this paper were generated using the Matlab scripts given inthe appendix. The program in appendix section A.3 provides an example implementation of theTemplate Fitting Method.4.2.1 The Bounding BoxThe most crucial portion of the Template Fitting Method is setting the bounding box on theexperimental data properly. The goal is to completely isolate any part of the template whichcorresponds to the oscillation frequency. This is the only parameter that the fit should be de-pendent on. Setting the bounds too large will result in the fit being dependent on random noiseor a divergent exponential, thus skewing the measured frequency. Bounding boxes which are toosmall or off-center may not provide enough data for an accurate fit.6This is dependent on sample size, area of the bounding box and speed of the computer used.18Chapter 4. Finding deltam with the Continuous Wavelet Transform(a) (b)Figure 4.4: (a) The function Pu,m(t) is decomposed into an exponential term in addition to(b) a sinusoidal term. The peaks in the coefficient diagram which are purely dependent on thecos(deltamt) term are what need to be isolated by setting the bounding box.The template created by the transform of Pu,m(t) has 3 components. High frequency noise canbe described at smaller scales. At larger scales the exponential portion of Pu,m(t) which divergesas the scale increases (see figure 4.4(a)) can be characterized. Between these is the sinusoid thatis deltam dependent (see figure 4.4(b)). The bounding box should encompass the peaks in this thirdfrequency range. Since the data contains a minimum of %25 swapped samples only the first fewpeaks will clearly be visible in the coefficient diagram. Although it is possible to successfullyextract the frequency using a fit on just the first peak, it will be much more sensitive to any biasin detection at short lifetimes, tau. Figure 4.5 shows an appropriate bounding box which containsthe minimum number of periods, 3, at the scale with the largest coefficient values.4.2.2 Calculating chi for deltam ExtractionBy truth matching the templates with the experimental data template a measure of the goodnessof the fit is calculated. In this case chi is usedchitmp =radicaltpradicalvertexradicalvertexradicalbtana=1bnsummationdisplayb=0(Wa,b,tmp -Xa,b)2sigma2b , (4.2)where tmp specifies the template, a and b are the wavelet scale and position coefficients with anand bn being their largest values. X is the wavelet coefficient of the experimental data and Wis the coefficient for the template. The variance sigma2b is constant over all scale values and is setat sigmab = radicalnb, where nb is the number of samples that decayed at the time corresponding to thecenter of the translated wavelet (b).Because the experimental data doesn?t always have the same number of samples as generated19Chapter 4. Finding deltam with the Continuous Wavelet TransformFigure 4.5: The combination of the cosine and exponential terms result in the coefficient diagramshown here. The bounding box, represented by the dashed line, isolates out the cosine term fromthe high frequency noise and diverging exponential term.data it is also necessary to first minimize the functionchitmp =radicaltpradicalvertexradicalvertexradicalbtana=1bnsummationdisplayb=0(Wa,b,tmp -alphatmpXa,b)2sigma2b (4.3)for the value of alphatmp. This fixes any problems that may arise from an inherent skew in the valuesof Wa,b,tmp and Xa,b.4.3 Results & DiscussionThe results presented here provide a proof of concept for the Template Fitting Method anddemonstrate its degree of accuracy along with any limitations found in the fitting process. Allexperimental data presented here was generated using the technique described in section 4.1.Plotting chitmp as a function of the underlying oscillation frequency for the template, as shownin figure 4.8, will typically result in a parabolic curve centered about the frequency of interest.The frequency deltam is taken to be the minimum of this curve with error bounds set at the 95%confidence level (CL) bounds of a cubic polynomial fit.The experimental data set whose frequency is to be extracted is generated with 20,000 sam-ple points and a randomly selected oscillation frequency of deltam = 6.418. Detection errors aresimulated with 28% of the mixed and un-mixed points being swapped, as well as 5% of the sam-ples being due to injected background errors. Once generated, the samples are binned and the20Chapter 4. Finding deltam with the Continuous Wavelet Transformcontinuous wavelet transform is taken over a broad sweep of scales to assist in deciding upon theplacement of the bounding box. The transform of both the Pu and Pm datasets is shown in figure4.6 with the bounding box highlighted.Since the Template Fitting Method is used in a recursive manner, two iterations are necessaryto isolate deltam. The first iteration is performed using nine templates evenly spaced over a frequencyrange of 3 < deltam < 12 oscillations per lifetime. The goal of the first iteration is to find a generalarea over which a more accurate second iteration can be performed. The templates are generatedusing the same number of sampling points as the experimental data, but with 30% swappedevents and 2% backround.The calculated value of chitmp is found to be minimized within a frequency range of 5 < deltam < 7oscillations per lifetime, as seen in figure 4.7. This is the range of frequencies used for the seconditeration, now using 100 evenly spaced templates to obtain a higher degree of accuracy whenextrapolating the minimum. The results from the second iteration are summarized in figure 4.8.The extracted frequencies are 6.52 ?0.14 for the Pm data and 6.44 ?0.14 for the Pu data, witha mean frequency of 6.48 ?0.14 oscillations per lifetime. The extracted frequency shows a goodagreement with the experimental dataset, demonstrating that the Template Fitting Method isapplicable to this kind of dataset.The results produced here are at the same (scaled) level of precision shown by the experimentalresults produced by CDF from Spring 2006 which place the B0s mass difference at deltamB0s =17.33ps-1 with a statistical error of +0.42-0.21 [2]. However it is necessary that the Template FittingMethod be compared directly to such methods on the same platform to ensure that it producesresults which are comparable.Test analyses produced thus far do not provide enough evidence that the Template FittingMethod is applicable to data sets with a high frequency. Results using the method demonstratedhere show that the precision begins to degrade at deltam > 12, at which point this method failsto extract a relevant frequency. This issue is caused by the overlay of noise with the mixingoscillations in the continuous wavelet transform?s coefficient diagrams. This makes it impossibleto set a bounding box which isolates the cos(deltamt) term (from the equation 2.20 for Pu,m) fromthe background errors.There is some hope that these problems may be overcome by averaging the amplitude ofthe coefficient diagrams made at the same frequency but with varying levels of noise. In thisway it may be possible to isolate out those characteristics in the coefficient diagrams which areunique to a specific frequency no matter what level noise is injected. It is also possible thattemplates generated with no noise, or generated from the cosine term alone may provide a betterfit. The adaptability of the Template Fitting Method is a promising quality, it also means thata significant amount of time is required to perfect the method used for analysing a dataset ascomplex as the flavor oscillations of B mesons.21Chapter 4. Finding deltam with the Continuous Wavelet Transform(a)(b)Figure 4.6: The bounding box for both the mixed, figure (a), and un-mixed, figure (b) testdatasets. The scales selected are 26 < am < 40 and 23 < au < 33 for the mixed and un-mixeddatasets respectively. In order to provide enough points to get a good measure of chitmp a scalespacing of deltaa = 0.5 is used to double the number of points generated.22Chapter 4. Finding deltam with the Continuous Wavelet Transform(a)(b)Figure 4.7: (a) The first iteration results of the Template Fitting Method, with chitmp values as afunction of the frequency for Pm(t). (b) The corresponding results for the first iteration of Pu(t).The Template Fitting Method shows a good fit in the vicinity of deltam = 6. The first iterationhas successfully identified the closest frequency to the one present in the simulated experimentaldata. 23Chapter 4. Finding deltam with the Continuous Wavelet Transform(a)(b)Figure 4.8: (a) The second iteration results for Pm(t). (b) The corresponding second iterationresults for Pu(t). Each point corresponds to a chitmp value, the solid curve generated from a cubicpolynomial fit. The minimum of this curve is the extracted frequency. 24Chapter 5ConclusionB mixing is a convenient event for measuring some of the CKM matrix elements. Because themass difference between its heavy and light states is directly applicable to finding Vtd and Vtsnumerous experiments are dedicated to improving the statistics for measuring its value. Howeverthe meson?s high rate of oscillation in addition to the difficulty involved in flavor tagging makethe results of these experiments troublesome to analyse.The Template Fitting Method has been developed for the purpose of improving the statisticalerror on the quantity deltam. In section 4.3 it has been shown that the Template Fitting Methodcan find the frequency of oscillation for data sets composed of 35% false points with a residual of0.062 which is well within the error bounds of ?0.14. It is clear that the Template Fitting Methodhas applications in extracting the B mixing frequency, however it has not yet been developed tothe point of completion. The scope of this project has not allowed for the development of a morefine tuned process.In the future there are many improvements to the Template Fitting Method which may makeit a far more accurate tool for data analysis. Using the de-noising capabilities of wavelets toextract out high frequency white noise may make it possible to correctly extract out larger valuesof deltam. Averaging a number of wavelets at the same frequency to create generic templates whichhighlight the noise independent qualities of the data set could reduce the effects of a low signal tonoise ratio. Applying non-rectangular bounding areas that better fit the sinusoidal peaks foundin the coefficient diagrams, may improve the calculation of chitmp.Experiments which are expensive to make and run commonly produce low statistics and con-tain large sources of error. This is a significant problem with many of the enormous experimentsproduced to verify the standard model. Measuring B meson decays requires an immense amountof resources. Despite that, it is possible that the biggest improvement to the accuracy of thesemeasurements will not come from the expansion of the experiment, but from new methods foranalysing the data they produce.25Bibliography[1] Belle Group, Nature 452, 332-335 (2008) .[2] W.-M. Yao et al. (Particle Data Group), J. Phys. G 33, 1 (2006).[3] M. Gell-Mann, A.Pais, Phys. Rev. 97, 1387 (1955) .[4] T.D. Lee and C. N. Yang, Phys, Rev, 104, 254 (1956).[5] J.H. Christenson et al., Phys. Rev. Lett. 13, 138 (1964).[6] D. Griffiths, ?Introduction To Elementary Particles?, John Wiley & Sons (1987).[7] C. Gay, Annu. Rev. Nucl. Part. Sci. 50, 577-641 (2000).[8] J. Rademacker, ?Evaluation of the LHCb RICH Detectors and a Measurement of the CKMAngle?, University of Oxford, 1-28 (2001)[9] L. Wolfenstein, Phys. Rev. D. 31, 2381 (1985).[10] A. Bondar, P. Pakhlov, A. Poluektov, Phys.-Usp. 50, 669-690 (2007).[11] A. Grossmann, J. Morlet, ?Decomposition of Hardy functions into square integrable waveletsof constant shape? SIAM J. Math., 15, 723-736 (1984)[12] Y. Meyer, R. Ryan, ?Wavelets Algorithms & Applications?,SIAM (1993).[13] W. Hardle, G. Kerkyacharian, D. Picard, A. Tsybakov, ?Wavelets, Approximations, andStatistical Applications?, Springer (1998).[14] P. Addison, ?The Illustrated Wavelet Transform Handbook?, IoP (2002).[15] C. Chui, ?An Introduction to Wavelets?, Academic Press Inc. (1992).26Appendix AProgramsAll programs presented here are written in the Matlab scripting language and will run in anystable Matlab Version 7.0 or better environment with the Wavelet Toolbox. The Curve FittingToolbox also available for Matlab is of considerable use when analysing the chi-frequency plots fora minimum with error bounds.A.1 Program to Generate Pu,m(t)function [pu_data,pm_data]=generate_master_data(samples,mxtau,gamma,freq)%Function [pu_data,pm_data]=generate_master_data(samples,mxtau,gamma,freq)%Returns the master datasets (pu_data,pm_data) each with a length samples.%The master curve is sampled anywhere from t=0 to t=mxtau.%Seed the random number generatorrand(?state?,sum(100*clock));%Initializationpu_data=zeros(1,samples);pm_data=zeros(1,samples);pu_peakval=0.5*gamma*exp(-gamma*0/freq)*2;pm_peakval=0.5*gamma*exp(-gamma*0.5/freq)*2;pu_master=inline(?0.5*g.*exp(-g*t).*(1+cos(2*pi*f.*t))?);pm_master=inline(?0.5*g.*exp(-g*t).*(1-cos(2*pi*f.*t))?);%Monte Carlo Method for Data Generationfor i=1:samples,%Generate un-mixed Data.while pu_data(i) == 0,%Select a random sampling point.rndtime=rand*mxtau;%Evaluate the function at the sampling point.master_time=pu_master(freq,gamma,rndtime);%Randomly drop valuesif rand*pu_peakval<master_time,27Appendix A. Programspu_data(i)=rndtime;endend%Generate Mixed Data, done in the same was as un-mixed data.while pm_data(i) == 0,rndtime=rand*mxtau;master_time=pm_master(freq,gamma,rndtime);if rand*pm_peakval<master_time,pm_data(i)=rndtime;endendend28Appendix A. ProgramsA.2 Program that Applies Simulated Detection Errorsfunction [pu_final,pm_final]=sim_error(pu_data,pm_data,masterpnts,...bgapnts,bgbpnts,smearfun,swaprate)%[pu_final,pm_final]=sim_error(pu_data,pm_data,masterpnts,...% bgapnts,bgbpnts,smearfun,swaprate)%run after generating master sample data using gen_master_data%datafile is the name of a .mat file containing at least:%2 variables of min size <1x100000 (double)>, called pm_data and pu_data%Arguements:%masterpnts - <100000, the number of points to use from the pu/pm master sets%bgapnts - <10000, the number of points to take from the background a curve%bgbpnts - <10000, the number of points to take from the background b curve%smearfun - usually ?0.1*t?, the time dependent function for smearing data%swaprate - the % (0<sr<1) of misses that should occur between pu and pm datarand(?state?,sum(100*clock));%Initializeload bga_samp %loads bga_noise_samp with 10000 points, the 1st background curveload bgb_samp %loads bgb_noise_samp with 10000 points, the 2nd background curvebga_master_sz=10000;bgb_master_sz=10000;bga_noise_samp=bga_noise_samp/max(abs(bga_noise_samp));bgb_noise_samp=bgb_noise_samp/max(abs(bgb_noise_samp));pm_master_sz=length(pm_data);pm=zeros(masterpnts+bgapnts+bgbpnts,1);pu=zeros(masterpnts+bgapnts+bgbpnts,1);%select points to use at random from master noisy setspmnormfact=(1/2)*max(abs(pm_data));punormfact=(1/2)*max(abs(pu_data));for i=1:masterpnts,pm(i)=pm_data(ceil(rand*pm_master_sz));pu(i)=pu_data(ceil(rand*pm_master_sz));endfor i=masterpnts+1:masterpnts+bgapnts,pm(i)=bga_noise_samp(ceil(rand*bga_master_sz))*pmnormfact;pu(i)=bga_noise_samp(ceil(rand*bga_master_sz))*punormfact;endfor i=masterpnts+bgapnts+1:masterpnts+bgapnts+bgbpnts,pm(i)=bgb_noise_samp(ceil(rand*bgb_master_sz))*pmnormfact;29Appendix A. Programspu(i)=bgb_noise_samp(ceil(rand*bgb_master_sz))*punormfact;endclear bga_noise_samp;clear bgb_noise_samp;clear pm_data;clear pu_data;%smear the data using the smear function arguementsigmafunc=inline(smearfun);shiftvals=randn(size(pm));shiftvals=shiftvals.*sigmafunc(pm);mxdat=max(pm);for i=1:length(pm)while(pm(i)+shiftvals(i)>mxdat)shiftvals(i)=randn(1);shiftvals(i)=shiftvals(i).*sigmafunc(pm(i));endendpm_smear=pm+shiftvals;clear pm;shiftvals=randn(size(pu));shiftvals=shiftvals.*sigmafunc(pu);mxdat=max(pu);for i=1:length(pu)while(pu(i)+shiftvals(i)>mxdat)shiftvals(i)=randn(1);shiftvals(i)=shiftvals(i).*sigmafunc(pu(i));endendpu_smear=pu+shiftvals;clear pu;%swap the data using miss rate arguementk=0;datsz=size(pm_smear);for i=1:datsz(1),for j=1:datsz(2),if(rand<=swaprate),k=k+1;temp=pu_smear(i,j);pu_smear(i,j)=pm_smear(i,j);30Appendix A. Programspm_smear(i,j)=temp;endendendpm_final=pm_smear;pu_final=pu_smear;31Appendix A. ProgramsA.3 Program to Extract deltam using the Template FittingMethodfunction frq = template_extraction(tfreqs,tbox_pm,tbox_pu,xfreq,dirc,...basename,generate_data,with_plots)%[frq,err] = template_extraction(template_freq,template_box,...% experimental_freq,directory,base_filename,generate_data,with_plots)%Generates or loads experimental data and templates to perform the template%fitting method to extract the oscillation frequency in the experimental data.%%Parameters:%frq,err: The frequency extracted along with the error associated with the% measurement.%tbox: A 3 vector specifying the bounding box to be used in generating the% templates and experimental data. Has the form:% [init_scale,del_scale,last_scale; init_posit,del_posit,last_posit]%tfreqs: A vector specifying the template frequencies to make the templates at%xfreq: The experimental frequency to attempt to find%dirc,basename: Strings. if generate_data==0 then the templates are loaded% from the directory/file ?dirc/basname_[pu|pm]_f<template freq>?% Otherwise this is the base filename used to save the% templates under.%generate_data: An int, if==3 then generates JUST templates & loads exp.% if ==2 then the script generates BOTH the template data and% experimental data required and saves for re-use. if==1, then generates% JUST the experimental data, and loads the template files. if==0, then% BOTH the template and experimental data are loaded.%with_plots: A boolean, if==1 then plots of each template are saved,% along with plots corresponding to the fitting process%Initializationscales_pm=tbox_pm(1,1):tbox_pm(1,2):tbox_pm(1,3);scpm_len=length(scales_pm);posits_pm=tbox_pm(2,1):tbox_pm(2,2):tbox_pm(2,3);pspm_len=length(posits_pm);scales_pu=tbox_pu(1,1):tbox_pu(1,2):tbox_pu(1,3);scpu_len=length(scales_pu);posits_pu=tbox_pu(2,1):tbox_pu(2,2):tbox_pu(2,3);pspu_len=length(posits_pu);prefix=strcat(?./?,dirc,?/?,basename,?_?);%Template Generation or Loadingswitch generate_data32Appendix A. Programscase 3 %gen temp, load exp%Load the experimental datasuffix=sprintf(?xdata_f%i?,dub2int(xfreq));name=strcat(prefix,suffix);load(name);%Generate temapltestemplates=gen_templates(tfreqs,scales_pm,posits_pm,...scales_pu,posits_pu,with_plots,prefix);case 2 %gen temp and exp%Generate experimental dataxdata=gen_xdata(xfreq,scales_pm,posits_pm,...scales_pu,posits_pu,with_plots,prefix);%Generate templatestemplates=gen_templates(tfreqs,scales_pm,posits_pm,...scales_pu,posits_pu,with_plots,prefix);case 1 %load temp, gen exp%Generate experimental dataxdata=gen_xdata(xfreq,scales_pm,posits_pm,...scales_pu,posits_pu,with_plots,prefix);%Load the template Datasuffix=?templatedata?;name=strcat(prefix,suffix);load(name);case 0 %load temp and exp%Load Experimental Datasuffix=sprintf(?xdata_f%i?,dub2int(xfreq));name=strcat(prefix,suffix);load(name);%Load the template Datasuffix=?templatedata?;name=strcat(prefix,suffix);load(name);otherwisereturn;end%Now the experimental data and template data are present, in:% xdata(pum,sc,pos) and templates(pum,fq,sc,pos)%Perfrom the Frequency Extraction using template truth matching.tcoef_vals=zeros(2,length(tfreqs)); %The template coefficientstchisq=zeros(2,length(tfreqs)); %Chi-square template-xdata fit value.%Precalculate the uncertainty for finding chi-squared33Appendix A. Programssigpm=sqrt(squeeze(xdata(1,1:scpm_len,1:pspm_len)*1.25));sigpu=sqrt(squeeze(xdata(2,1:scpu_len,1:pspu_len)*1.25));%options used for minimizing the template coefficient.opts=optimset(?MaxFunEvals?,1e5,?MaxIter?,1e5);cnt=0; %counter used for indexingtmp=zeros(max(scpm_len,scpu_len),max(pspm_len,pspu_len));for i=1:length(tfreqs)cnt=cnt+1;tmp(1:scpm_len,1:pspm_len)=squeeze(templates(1,i,1:scpm_len,1:pspm_len));[tcoef_vals(1,cnt),tchisq(1,cnt)]=fminsearch(@(tcoef)...sqrt(sum(sum(((tmp(1:scpm_len,1:pspm_len)-...tcoef.*squeeze(xdata(1,1:scpm_len,1:pspm_len)))..../sigpm(1:scpm_len,1:pspm_len)).^2))),1,opts);tmp(1:scpu_len,1:pspu_len)=squeeze(templates(2,i,1:scpu_len,1:pspu_len));[tcoef_vals(2,cnt),tchisq(2,cnt)]=fminsearch(@(tcoef)...sqrt(sum(sum(((tmp(1:scpu_len,1:pspu_len)-...tcoef.*squeeze(xdata(2,1:scpu_len,1:pspu_len)))..../sigpu(1:scpu_len,1:pspu_len)).^2))),1,opts);end%Save the template comparison data.suffix=sprintf(?templatechisq_f%i?,dub2int(xfreq));name=strcat(prefix,suffix);save(name,?tchisq?);[val,pm_ind]=min(tchisq(1,:));[val,pu_ind]=min(tchisq(2,:));if(with_plots)mn=0;mx1=max(tchisq(1,:));mx2=max(tchisq(2,:));figure; plot(tfreqs,tchisq(1,:),?o?);hold on; plot([xfreq,xfreq],[mn,mx1],?--?);xlabel(?Frequency (\Delta m) in Oscillations/Lifetime?);ylabel(?Chi-squared?);suffix=?pm_extractchisq_f?;name=strcat(prefix,suffix);saveas(gcf,name,?fig?);close(gcf);figure;plot(tfreqs,tchisq(2,:),?o?);hold on; plot([xfreq,xfreq],[mn,mx2],?--?);xlabel(?Frequency (\Delta m) in Oscillations/Lifetime?);ylabel(?Chi-squared?);34Appendix A. Programssuffix=?pu_extractchisq_f?;name=strcat(prefix,suffix);saveas(gcf,name,?fig?);close(gcf);endfrq(1)=tfreqs(pm_ind);frq(2)=tfreqs(pu_ind);endfunction integer = dub2int(dub)%function integer = dub2int(double)%Converts a double to an integer by multiplying by 10 to get rid of decimal,%and takes the absolute value.while(mod(dub,1)>0 && dub<1e3)dub=dub*10;endinteger=abs(dub);endfunction xdata=gen_xdata(xfreq,scales_pm,posits_pm,...scales_pu,posits_pu,with_plots,prefix)gen_samp=100000;datasetsamp=49000;bgasamp=500;bgbsamp=500;hpnts=2000;mxtau=5;gamma=1;smear=?0.1*t?;swap=0.3;[pu_data,pm_data]=generate_master_data(gen_samp,mxtau,gamma,xfreq);[pu_td,pm_td]=sim_error(pu_data,pm_data,datasetsamp,...bgasamp,bgbsamp,smear,swap);xdata=zeros(2,max(length(scales_pm),length(scales_pu)),...max(length(posits_pm),length(posits_pu)));xdata_hist(1,:)=hist(pm_td,hpnts);xdata_hist(2,:)=hist(pu_td,hpnts);abs(cwt(xdata_hist(1,posits_pm),scales_pm,?morl?));35Appendix A. Programsxdata(1,1:length(scales_pm),1:length(posits_pm))=...abs(cwt(xdata_hist(1,posits_pm),scales_pm,?morl?));xdata(2,1:length(scales_pu),1:length(posits_pu))=...abs(cwt(xdata_hist(2,posits_pu),scales_pu,?morl?));if(with_plots)%Generate P_m plot and save/closefigure;imagesc(squeeze(xdata(1,1:length(scales_pm),1:length(posits_pm))));xlabel(?position (unitless)?);ylabel(?scales?);suffix=sprintf(?xdata_pm_f%d?,dub2int(xfreq));name=strcat(prefix,suffix);saveas(gcf,name,?fig?);close(gcf);%Generate P_u plot and save/closefigure;imagesc(squeeze(xdata(2,1:length(scales_pu),1:length(posits_pu))));xlabel(?position (unitless)?);ylabel(?scales?);suffix=sprintf(?xdata_pu_f%d?,dub2int(xfreq));name=strcat(prefix,suffix);saveas(gcf,name,?fig?);close(gcf);endsuffix=sprintf(?xdata_f%i?,dub2int(xfreq));name=strcat(prefix,suffix);save(name,?xdata_hist?,?xdata?);endfunction templates=gen_templates(tfreqs,scales_pm,posits_pm,...scales_pu,posits_pu,with_plots,prefix)gen_samp=100000;datasetsamp=49000;bgasamp=500;bgbsamp=500;hpnts=2000;mxtau=5;gamma=1;smear=?0.1*t?;swap=0.3;mxscal=max(length(scales_pm),length(scales_pu));mxposi=max(length(scales_pm),length(scales_pm));%Generate the templates36Appendix A. Programstdata_hist=zeros(2,length(tfreqs),hpnts); %[pm or pu, frequency, tau]templates=zeros(2,length(tfreqs),mxscal,mxposi);%[pm|pu,freq,scale,position]%Generate P_u and P_m curves.for i=1:length(tfreqs)[pu_data,pm_data]=generate_master_data(gen_samp,mxtau,gamma,...tfreqs(i));[pu_td,pm_td]=sim_error(pu_data,pm_data,datasetsamp,...bgasamp,bgbsamp,smear,swap);tdata_hist(i,1,:)=hist(pm_td,hpnts);tdata_hist(i,2,:)=hist(pu_td,hpnts);end%Free up some space.clear pu_data pm_data pu_td pm_td;%Generate the templates for each P_u and P_m curve.for i=1:length(tfreqs)templates(1,i,1:length(scales_pm),1:length(posits_pm))=...abs(cwt(tdata_hist(i,1,posits_pm),scales_pm,?morl?));templates(2,i,1:length(scales_pu),1:length(posits_pu))=abs(cwt(tdata_hist(i,2,posits_pu),scales_pu,?morl?));endif(with_plots)for i=1:length(tfreqs)%Generate P_m plot and save/closefigure;imagesc(squeeze(templates(1,i,1:length(scales_pm),...1:length(posits_pm))));xlabel(?position (unitless)?);ylabel(?scales?);suffix=sprintf(?template_pm_f%d?,dub2int(tfreqs(i)));name=strcat(prefix,suffix);saveas(gcf,name,?fig?);close(gcf);%Generate P_u plot and save/closefigure;imagesc(squeeze(templates(2,i,1:length(scales_pu),...1:length(posits_pu))));xlabel(?position (unitless)?);ylabel(?scales?);suffix=sprintf(?template_pu_f%i?,dub2int(tfreqs(i)));37Appendix A. Programsname=strcat(prefix,suffix);saveas(gcf,name,?fig?);close(gcf);endendsuffix=?templatedata?;name=strcat(prefix,suffix);save(name,?templates?);clear tdata_hist;end38

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.52966.1-0085956/manifest

Comment

Related Items