In this report the results from a 10 week internship are presented.  The goal of the internship was to apply different unfolding approaches to conduct measurements of energy spectra from data aquired by FACT, the First G-APD Cherenkov Telescope. FACT is the first operational telescope of its kind, employing a camera equipped with silicon photo multipliers (G-APD aka SiPM) to primarily detect gamma rays. Improving the unfolding method can help with better interpretation of the data and more accurate physics results without the need for new equipment or more observations. The approaches tested during this internship range from simplistic matrix inversion to an improvement over of the previous standard (TRUEE)

Bieker, Jacob

Brügge, Kai

Börner, Mathis

Nöthe, Maximillian

Eldorado - Ressourcen aus und für Lehre, Studium und Forschung

RISE Germany InternshipUnfolding FACT DataTechnicalReport Jacob Bieker, Mathis Börner, KaiBrügge, Maximillian Nöthe09/2017technische universitätdortmundPart of the work on this technical report has been supported by Deutsche Forschungs-gemeinschaft (DFG) within the Collaborative Research Center SFB 876 "ProvidingInformation by Resource-Constrained Analysis", project C3.Speaker: Prof. Dr. Katharina MorikAddress: TU Dortmund UniversityJoseph-von-Fraunhofer-Str. 23D-44227 DortmundWeb: http://sfb876.tu-dortmund.de1 IntroductionIn this report the results from a 10 week internship are presented. The goal of theinternship was to apply different unfolding approaches to conduct measurements of energyspectra from data aquired by FACT, the First G-APD Cherenkov Telescope. FACT isthe first operational telescope of its kind, employing a camera equipped with silicon photomultipliers (G-APD aka SiPM) to primarily detect gamma rays. Improving the unfoldingmethod can help with better interpretation of the data and more accurate physics resultswithout the need for new equipment or more observations. The approaches tested duringthis internship range from simplistic matrix inversion to an improvement over of theprevious standard (TRUEE [7]).2 UnfoldingThe goal of unfolding is to determine which input distribution f (y) of a desired quantityy leads to the measured distribution g (x) of the observable x. For continuous distributionthis can be written as an integral equationg (x) =∫A (x, y) f (y) dy. (1)In this equation A (x, y) is the so-called detector response function and returns the condi-tional probability to measure a value x given a value for y. In most applications discretedistributions f and g are used and (1) transforms to~g = A~f. (2)The discretization for f is often motivated by the physics that is being investigated and isusually a simple equidistant binning in y or a spline representation of the distribution (e.g.in TRUEE). For an equidistant binning ~f contains the number of events in the differentbins and for a spline representation it the coefficient vector for the splines. In (2) theformer detector response function A (x, y) is replaced by the detector response matrix A.Analog to A (x, y) the detector response matrix contains the conditional probabilities tomeasure an event in a bin i of the observable distribution ~g given the event is in bin j ofthe desired distribution ~f .As shown later the challenges of solving (2) is driven by the detector response matrix A.The characteristics of this matrix mainly depends on the binning chosen for the observablespace. In subsection 2.1 the different combinations of unfolding and binning approachesare presented and used to illustrate the evolution from the most basic approach to thethe new approach that improves upon the current standard (TRUEE) used for FACTdata [8].2.1 Different Unfolding ApproachesIn this section the different approaches and the challenges associated with them areexplained.12.1.1 Inverse MatrixThe simplest way to unfold data is to invert the detector response matrix. The invertedmatrix can then be multiplied with the measured signal ~g to unfold back to the realspectrum ~f using ~f = A−1~g. Depending on the uncertainty in the measured distribution,this will lead to large fluctuations in the solution. The cause for those fluctuations can beshown by transforming (2) into a basis in which A is a diagonal matrix [1] (A = UDUᵀ):~g = A~f → Uᵀ~g = DUᵀ ~f → ~c = D~bMatrix D is a diagonal matrix with A’s eigenvalues on the diagonal1. Therefore, theunfolding via inverting the matrix in this space is simply:bi =1λici. (3)In (3) the eigenvalues of matrix A are inverted and therefore small eigenvalues becomelarge factors for their corresponding coefficients ci. With increasing i the coefficientsbecome less and less significant. This is what shows up as large fluctuations in the solutionand is sometimes referenced as fitting to noise. This characteristic of the unfoldingemerges from ambiguous connections between a bin in the observable space and a binof the sought-after quantity. The so-called condition number κ of matrix A is a usefulmeasure of the uncertainty of the measurement and how strong the solution will bedominated by noise. It is defined by the ratio between the biggest and the smallesteigenvalues κ = λ0/λN .To tackle the problem of fluctuating solutions multiple options are available: Increasestatistics to get more significant coefficients, ignore small eigenvalues (regularization2) orfind a way to increase the eigenvalues of matrix A. In most analysis all available data isused, so using regularization or increasing the condition of matrix A are the only feasibleoptions.2.1.2 Singular Value DecompositionSingular Value Decomposition (SVD) [6] improves on inverting the matrix by allowingfor A’s that are not square. Without that restriction, more information can be utilized,generally leading to a better condition3 of matrix A.2.1.3 Poissonian LikelihoodIt is known that the underlying distribution in each observable bin gi is a Poisson dis-tribution. This knowledge is not taken into account by neither the matrix inversion nor1We assume that the eigenvalues are ordered by their size2As shown later for this unfolding no regularization is needed, it is not covered in full detail in thisreport.3For rectangular matrices the condition is defined with the singular values instead of the eigenvalues.2by the SVD approach. This can lead to biased solutions when many bins contain small(≤ 10) numbers of events.The poissonian distribution can be taken into account when the problem is formulatedas a likelihood function:L(~g∣∣~f) = m∏i=1(ĝgiigi!exp (−ĝi))· 1√(2π)n det (τ1)exp(−12τ ~fᵀCᵀ1C~f)︸ ︷︷ ︸Tikhonov Regularization. (4)To find the maximum of the likelihood function (4) two different methods were used.Gradient DescentThe gradient descent (Newton method) approach implemented in this study analyticallycalculates the negative log likelihood, L, gradient ~h of the negative log likelihood, andHessian matrix H. The Newton method is an iterative method starting with ~f0. Foreach iteration the gradient and Hessian are calculated to determine next position ~fi+1 =~fi −H−1 ·~h · γ 4. We run this iteration for a fixed (large) number of steps. The negativelog likelihood S, its gradient ~h and Hessian H used are:For minimization:S(~g∣∣~f) = m∑i=1((A~f)i− gi ln((A~f)i))+12τ ~fᵀC′ ~fGradient:dSdfk= hk =m∑i=1Ai,k − giAi,kn∑j=1Ai,jfj+ τ n∑i=1fiC′k,iHesse-Matrix:d2Sdfkdfl= Hk,l =m∑i=1giAi,kAi,l(n∑j=1Ai,jfj)2+ τC ′k,lThe minimum of the negative likelihood provides the most likely spectrum leading tothe measured distribution in the observable space. To obtain the uncertainties herethe inverse Hessian (parabolic approximation of the likelihood space) is used. This is areasonable approximation for bins in f with large statistics, but might break down forbins closer to zero, because we can’t have negative bin counts and the likelihood spacestarts to become asymmetric.4The factor γ is used to increase/decrease the step size.3Markov Chain Monte CarloThe final unfolding method tested during this study is Markov Chain Monte Carlo(MCMC) [5]. The idea of MCMC is to perform a random walk in the likelihood space.For each step of the walk is randomly drawn via rejection sampling based on the like-lihood function. This leads to a sample of the posterior probability distribution of thelikelihood function.The sample can be used to determine the best fit and different uncertainty levels. Oneadvantage of the MCMC approach above the gradient descent is that no approximationsare made about the likelihood space. This leads to reasonable error contours for smallentries in ~f . Another great advantage is that the sample of the likelihood also provides astraightforward way to compare the compatibility between a given spectrum ~fTest and theunfolded result. To calculate the compatibility one simply has to calculate the likelihoodvalue for the test spectrum and compare with the likelihood values for each example inthe MCMC sample. This provides a p-value for the hypothesis that the test spectrumand the unfolded spectrum have the same underlying distribution. For all MCMC resultsin this report the emcee [4] package was used.2.2 Binning MethodsAs mentioned in section 2 the binning of the observable space is crucial to achieve a goodcondition for matrix A. There are no limitations or constraints on the bins’ shapes or inthe dimensionality of the binning. TRUEE for example allows up to three observables.In all these observables the number of bins is set and equidistant binning is used for eachobservable. The danger of increased dimensionality is that the binning can create a lotof empty or sparsely populated bins. This can lead to two major problems. Bins withfew entries give very little constraint on the fit, and too few simulated events in a binintroduce a large uncertainty in the entries of matrix A due to limited statistics.2.2.1 EquidistantOne way to deal with the previously mentioned problems of sparsely populated bins is bymerging bins in a equidistant binning. There are various ways to do this. Ones used hereto test the effectiveness of this were merging the closest bins with low counts until themerged bins contained a set minimum number of events, and merging the bins with thelowest counts together until at least the minimum number of events is contained withinthem. Example of equidistant binning in two dimensions is shown in Figure 1.2.2.2 Decision Tree BinningA more sophisticated way [2] to obtain a binning for the problem utilizes a decisiontree [3]. The idea is to train a decision tree classifier to classify events into the bins of ~f .The training of a decision tree is an recursive process. In each step a cut is performedto achieve the best separation of the different classes. For the two new datasets created4(a) Equidistant binning. (b) Merging closest bins to have at least 10events in each bin.Figure 1: Equidistant without and with merging closest bins until a threshold is reached.from the cut, an new cut is performed and so on. The recursion stops either when thepart of the dataset is purely one class or an external condition is fulfilled. The stoppingpoints are so-called leafs of the tree and the path to the leafs is a sequence of cuts inthe observables. Each leaf represents a rectangular space in the observable space and thewhole space is covered by the leafs.This approach has multiple advantages over the equidistant binning. The biggest advan-tage is that the binning is explicitly optimized for the problem. The decision builds themodel to differ between the different bins in ~f . This should directly lead to an optimizedcondition of the matrix A.Another advantage is that the training process can be controlled with external conditions.E.g. the recursion of the training process can be stopped when only k events are left inthe dataset. This prevents the binning from having empty/weakly populated bins. Avisualization of the binning with two observables is shown in Figure 23 ResultsThe various different methods of unfolding the data are presented below. All methodswere applied to simulations for the FACT telescope. The wanted quantity is the energyof the γ-particles. For every unfolding the binning in the energy was equidistant inlog10(Eγ/ [GeV]) between 2.4 and 4.2.Used DataWhen using the FACT Monte Carlo data, the initial dataset consisted of 391 904 puresignal events with 253 different attributes. In this study, only simulated data without5Figure 2: Equidistant without and with merging closest bins until a threshold is reached.background was used. In addition, only 18 of the attributes were used to train the decisiontree classifier for the binning. To simulate an actual measurement run, the events arerandomly split into three groups. One group of roughly 5000 events is used as the testingdata. All the other data is split into twenty percent for training the decision tree classifier,and eighty percent to build the detector response matrix. For all observable binnings itis ensured to have at least 10 test events in each bin.3.1 SVD UnfoldingAs discussed in subsubsection 2.1.1 the condition of detector response matrix A providesdeep insight into the characteristics of the unfolding. Figure 3 shows the singular valuesand the corresponding condition number for five different binnings. The worst condition(κ = 90 223.0) is obtained with an equidistant binning in one observable with low correla-tion to the energy Eγ . Using an observable with a known high correlation to the energythe condition is improved by more than two orders of magnitude (κ = 169.0). Whentwo observables (visualization Figure 1a) are used the condition can be further improved(κ = 84.6). The full potential of the decision tree-based binning can be utilized whenall 18 observables are used (κ = 23.0). Using the same two observables as used in theequidistant binning gives a condition number of κ = 67.0 for the tree based binning.6Figure 3: Normed singular values and conditions for different binnings.In Figure 4 the results of the SVD unfolding for the binnings are shown. The differencesof the condition can directly be seen in the fluctuations of the results. Only the unfoldingusing the tree binning with all observables shows no obvious fluctuations. This is aindication that such a binning removes or reduces greatly the regularization needed tounfold FACT data.3.2 Poissonian LikelihoodFigure 5 shows the path of the gradient descent approach. The algorithm is capableof finding the maximum of the likelihood function. Looking at the likelihood samplefrom the MCMC in the background of Figure 5 indicates that for the very first bin aparabolic approximation of the likelihood space might not be justified. The distributionis asymmetric because it is close to zero.In the figures 6, 7, 8, 9 the MCMC samples and best fits for the same and for differentbinning as shown in Figure 3 are depicted.Similar to the SVD unfolding, all binnings, except the tree binning using all 18 observ-ables, show no clear maximum in the likelihood space and a large variance in the samplefor all bins.7Figure 4: Unfolding results from SVD unfolding for the different detector response ma-trices obtained with the different binning approaches.84 ConclusionIn the first chapters the benefits and drawbacks of four different unfolding approacheswere discussed. In the application of those approaches it clearly showed that using aMCMC sampler with poissonian likelihood for the unfolding provides the best results.The difference between the SVD and the gradient descent unfolding (both with 18 ob-servable tree binning) were the uncertainties on the unfolded spectra. The first bin in theunfolding illustrated why an approach that doesn’t rely on a parabolic approximation toobtain the uncertainties should be favored. Even the best fit in the bin was ≈ 80 resultsclose to 0 were not excluded. For results close to 0 the parabolic approximation can beexpected to break down.Another insight gained from the studies is that using the decision tree based binningapproach (subsubsection 2.2.2) with all 18 observables in the training improves the con-dition of the problem up to a point where no regularization might be needed. This isa promising finding, but needs to be validated in more extensive studies. A first checkwould be to repeat the unfolding of simulated events multiple times (O(1000)) and in-vestigate if the result is unbiased. To do that for each unfolding, the p-value for thedistribution of the simulation has to be calculated as described in subsection 3.2. For anunbiased solution the p-values should have a uniform distribution.After showing that the unfolding with is unbiased the next logical step is to apply theapproach on real data.References[1] Volker Blobel. “An unfolding method for high energy physics experiments”. In: arXivpreprint hep-ex/0208022 (2002).[2] M. Börner et al. “Measurement/Simulation Mismatches and Multivariate Data Dis-cretization in the Machine Learning Era”. In: ADASS XXVII. Ed. by TBD. Vol. TBD.ASP Conf. Ser. San Francisco: ASP, 2018, TBD.[3] L. Breiman et al. Classification and Regression Trees. Monterey, CA: Wadsworthand Brooks, 1984.[4] Daniel Foreman-Mackey et al. “emcee: the MCMC hammer”. In: Publications of theAstronomical Society of the Pacific 125.925 (2013), p. 306.[5] Charles J Geyer. “Markov chain Monte Carlo maximum likelihood”. In: (1991).[6] Andreas Hoecker and Vakhtang Kartvelishvili. “SVD approach to data unfolding”.In: arXiv preprint hep-ph/9509307 (1995).[7] Natalie Milke et al. “Solving inverse problems with the unfolding program TRUEE:Examples in astroparticle physics”. In: Nuclear Instruments and Methods in PhysicsResearch Section A: Accelerators, Spectrometers, Detectors and Associated Equip-ment 697 (2013), pp. 133–147.9[8] Maximilian Nöthe et al. “FACT Performance of the First Cherenkov Telescope Ob-serving with SiPMs”. In: Proceedings of the 35th International Cosmic Ray Confer-ence. PoS(ICRC2017)791. Proceedings of Science, 2017. url: https://pos.sissa.it/301/791/pdf.10Figure 5: Visualization of the path in the gradient descent minimization. The blue pointsof every 10-th step of the 500 iterations. In the background the likelihood sample obtainedwith the MCMC approach is shown.11Figure 6: Result of the MCMC unfolding using equidistant binning with a single observ-able. The plots on the diagonal show the distribution of the samples for the differentbins in ~f . The non-diagonal plots show the two dimensional distribution of two bins. Allaxes show number of events in bin. The marked points are the best fit.12Figure 7: Result of the MCMC unfolding using equidistant binning with two observables.The plots on the diagonal show the distribution of the samples for the different bins in~f . The non-diagonal plots show the two dimensional distribution of two bins. All axesshow number of events in bin. The marked points are the best fit.13Figure 8: Result of the MCMC unfolding using decision tree based binning with twoobservables. The plots on the diagonal show the distribution of the samples for thedifferent bins in ~f . The non-diagonal plots show the two dimensional distribution of twobins. All axes show number of events in bin. The marked points are the best fit.14Figure 9: Result of the MCMC unfolding using decision tree based binning with all (18)observables. The plots on the diagonal show the distribution of the samples for thedifferent bins in ~f . The non-diagonal plots show the two dimensional distribution of twobins. All axes show number of events in bin. The marked points are the best fit.15

RISE Germany Internship: Unfolding FACT Data

Abstract

Similar works

Full text

Available Versions

Eldorado - Ressourcen aus und für Lehre, Studium und Forschung