A Brief Note of Analyzing and Plotting ν μ Disappearance in SBN Detector under ROOT Framework
AA Brief Note of Analyzing and Plotting ν µ Disappearance inSBN Detector under ROOT Framework
Castaly FanJanuary 2020
Abstract
This is a brief technical note of analyzing the ν µ disappearance in SBN detector. We hereprovide a kind of method of plotting the histograms and the heat map plot. We will explorethe properties via ROOT framework supported by CERN. The main language we used is C++. Contents χ values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118.4 Appendix 4: Heat map plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 a r X i v : . [ phy s i c s . d a t a - a n ] O c t Background
For searching sterile neutrinos, one way is to look over the disappearence (reduction) of muonneutrinos ( ν µ ) in a muon neutrino beam. The data of detector properties we have is shown asbelow: Table 1: Data of detector properties[1]Detector Baseline Active LAr mass POT (protons on target)SBND 110m 112 tons 6 . × MicroBooNE 470m 89 tons 13 . × ICARUS 600m 476 tons 6 . × The cross section σ ( E ν ) for CCQE ν mu interactions on Ar is shown as the left figure below.And the initial neutrino flux φ ( E ν ) without oscillations at MicroBooNE is shown as the right figurebelow. (a) cross section (b) flux at MicroBooNE Figure 1: The figure of cross section and flux.Given the number of argon nuclei N Ar , we can calculate the event rate N ev ( E ν ) = σ ( E ν ) × N Ar × [ φ ( E ν ) × ∆ E ν ] × P OT (1)Armed with the equation 1 and the data of Table 1, we can start the task step by step.
Equation 1 shows that the event rate depends on the number of 4 Ar , which can be expressed as N Ar = M × N A W (2)where M is the mass in gram, W is the atomic weight, and N A is the number of nuclei per mole(namely, 6 . × ). Now, let us consider N Ar of each detectors. Given the LAr mass of Table 1,we obtain: 2 . × for SBND • . × for MicroBooNE • . × for ICARUSThen, we can create the event rate’s histogram according to Equation 1.Back to the Figure 1, the cross section for each detector is the same, whereas the curve of theflux would be scaled by the factor r since the flux falls off with the distance r from the target.For instance, given that the flux at MicroBooNE depends on the distance 470 m, the flux at SBNDwould be scaled by (470 / m and the flux at ICARUS would be scaled by (470 / m.Given that Equation 1, the unit should be considered in order to create the event rate’s his-togram. So we can simply write down the event rate equation in terms of the units:[ N ] = [ cm Ar ] × [ Ar ] × [ 1 M eV · P OT · m ] × [ M eV ] × [ P OT ] (3)Figure 2: Event rate’s histogram. The curves from up to down are respectively corresponding toSBND, ICARUS, and MicroBooNE.Here, [ N ] is denoted by the number of event. Then we have to notice the unit should benormalized, in particular, [ cm ] is for cross section whereas [ m ] is for flux.As a result (See Appendix 1), the histogram of event rate for each detector can be created asFigure 2. 3y calculating the total number of events according to the Figure2, we can obtain the value3 . × , which is almost matched with the number of ν µ events during ν µ n → µ − p processbased on the SBN proposal (namely, 3,122,600). To explore the disappearance of ν µ , we have to know the probability where ν µ remain, denotedby P ν µ → ν µ . In general, if P ν µ → ν x means the probability that ν µ turn into ν e , ν τ , or ν s (sterileneutrinos), then P ν µ → ν µ would be reasonably equal to 1 − P ν µ → ν x . Specifically, we have a formuladescribing ν µ disappearance as below[2]: P ν µ → ν µ = 1 − sin (2 θ µµ ) sin (cid:18) ∆ m L E ν (cid:19) (4)The equation above is taken by natural units. Thus, we can write down the form in SI units: P ν µ → ν µ = 1 − sin (2 θ µµ ) sin (cid:18) .
27 ∆ m LE ν (cid:19) (5)where θ is the mixing angle, L is the oscillation distance in kilometers, and E is the neutrinoenergy in GeV.Figure 3: Event rate’s histogram with oscillations. The curves of SBND, ICARUS, andMicroBooNE are same as Fig 2, whereas have subtle differences between both diagrams.4y arbitrarily setting sin (2 θ µµ ) = 0 . m = 1 eV , the factor of probability can be simplywritten as P ν µ → ν µ = 1 − . (cid:18) . LE ν (cid:19) (6), which shows the probability of oscillations depends on the distance L and the neutrinos energy E ν .To create the event rate histogram with neutrino oscillations, we can consider the factor ofequation 6. The distance L is corresponding to each detector’s baseline (see Table 1), as for theenergy E ν , it is given by the flux histogram as Figure 1. Now, we can obtain the event ratehistogram with oscillations, which is shown as Figure 3 (See Appendix 2).By looking at the Figure 3 roughly, it is hard to see obvious differences between that with theone without oscillations (Figure 2). However, there are actually few subtle differences between bothhistograms. Now what we have to do is to analyze the details of both curves via statistical ways. Statistically, to do data fitting work, we need to use least squares in order to minimize the sum ofsquared errors and residuals[3]. Consider a function y i , i = 1 , ..., N as a function of another variable x i without error. Each y i has a different unknown mean λ i and a known variance σ i . Then, supposethe true value λ is the function of x and depends on unknown parameters θ = ( θ , ..., θ m ), that is, λ = λ ( x ; θ ). Now, consider the χ ( θ ) quantity: χ ( θ ) = N (cid:88) i =1 ( y i − λ ( x i ; θ )) σ i (7)In other words, y , ..., y N are measured with errors σ , ..., σ N at each x ( x , ..., x N ) value withouterrors. The true value λ i = λ ( x i ; θ ), where θ is adjusted to minimize the value of χ given by theequation above.Assume un-oscillated event rate histogram is y i and oscillated one is another function λ . Forthe variance σ i , consider it as λ , since we regard the case as a Poisson distribution (i.e. the meanequal to the variance). Thus, the chi-squared value of one detector can be simply obtained from: χ = N (cid:88) i =1 ( y i − λ i ) λ i (8)As for the total chi-squared value for three detectors can be written as the sum χ total = χ a + χ b + χ c (9)where the index a, b, and c can be represented as the detector SBND, MicroBooNE, and ICARUSrespectively.Then we can print out the χ value for each detector (Appendix 3). we get: • SBND: χ a = 0 . • MicorBooNE: χ b = 1 . ICARUS: χ c = 3 . χ value is χ total = χ a + χ b + χ c = 5 . Given that the χ value for each detector, we can make a heat map plot to show whether theexperiment can be consistent with or rule out statistically. To realize it, we must make a 2-D plot,that is, TH2D in ROOT framework (Appendix 4).Ultimately, the plot on canvas is shown as Figure 4. The amazing result is important in testingthe sensitivity for each detector. Note that the axis is set as logarithmic value. For the formal heatmap plot, the curves can be shown more obviously, which is shown as Figure 5.Figure 4: Heat map plot given from each detector’s χ value. The contour curves can be seen asthe sensitivity prediction for the detectors.6igure 5: The formal heat map plot in SBN proposal (arXiv:1903.04608). Compared with bothplots, the curves are almost similar so that the result would be correct[2]. To create the 2-D heat map plot, the first step is to calculate event rate according to cross section’sand flux’s data. Then, considering neutrino oscillations, we can see the slightly different betweenboth event rate’s histograms. To know the difference between null hypothesis (non-oscillated) withthe oscillated one, we must figure out the χ value. Finally, compiling the result into the originalloop, we can attain the heat map plot with the curves. It is noteworthy that physicists can lookover the sensitivity from such contour plot. In this way, it can let us know the experimental detailsof ν µ disappearance in SBN detector. I am grateful to Prof. Andrew Mastbaum for giving me a precious opportunity to be engaged inSBN program’s research. It is a wonderful tour for me to learn plenty of technological skills as wellas data analysis work, which is crucial to my physics career in the near future.7 eferences [1] Pedro A. N. Machado, Ornella Palamara, David W. Schmitz.
The Short-Baseline NeutrinoProgram at Fermilab . (2019) arXiv:1903.04608.[2] R. Acciarri, C. Adams, R. An, C. Andreopoulos, A.M. Ankowski, et al . A Proposal for a ThreeDetector Short-Baseline Neutrino Oscillation Program in the Fermilab Booster Neutrino Beam .(2015) arXiv:1503.01520.[3] Glen Cowan.
Statistical Data Analysis . (1998) New York: Oxford University Press.8
Appendices
Here is the collection of some related codes with respect to ROOT framework.
To create the plot for event rate histograms, we can log into the ROOT framework and run thefollowing command (some ”//” means the note for the command line below): void draw_histogram ( TH1D * h, TGraph * cross_section , double parameters ) { std :: cout << parameters << std :: endl ; // for every bin in x- axis ( GeV ) for ( int i = 1; i <= h-> GetNbinsX (); ++i) { float binCenter = h-> GetBinCenter (i); // value of cross section float y1 = cross_section -> Eval ( binCenter ); // value of flux float y2 = h-> GetBinContent (i); // eventrate formula float N = y2 * y1 * parameters ; // assign the y value back to the plot h-> SetBinContent (i,N); // debug output // std :: cout << N << std :: endl ; } h-> GetYaxis () -> SetTitle (" Number of Events / GeV "); h-> Draw (" SAME "); std :: cout << h-> Integral (" width ") << std :: endl ; } void eventrate () { // find the directory where cross section file is TFile * cross_section_file = TFile :: Open (" cross_section . root "); cross_section_file ->ls (); // read the cross section file TGraph * cross_section = ( TGraph *) cross_section_file -> Get (" qel_cc_n ;1"); // draw the cross section figure // cross_section -> Draw (); // find the directory where flux file is TFile * flux = TFile :: Open (" flux . root "); flux ->ls (); // read the flux file TH1D * h = ( TH1D *) flux -> Get (" numu_CV_AV_TPC "); // create three clones for three cases TH1D * h_sbnd = ( TH1D *) h-> Clone (" h_sbnd "); TH1D * h_microboone = ( TH1D *) h-> Clone (" h_microboone "); TH1D * h_icarus = ( TH1D *) h-> Clone (" h_icarus "); draw_histogram ( h_sbnd , cross_section ,
1E -38/1 E4 *1.6856 E30 /0.05*6.2 E20 /1 E6 *(470.0/110) *(470.0/110) ); draw_histogram ( h_microboone , cross_section ,
1E -38/1 E4 *1.3395 E30 /0.05*13.2 E20 /1 E6); draw_histogram ( h_icarus , cross_section ,
1E -38/1 E4 *7.1638 E30 /0.05*6.2 E20 /1 E6 *(470.0/600) *(470.0/600) ); } The command here is slightly different from the first one, that is, we have introduced ”baseline” asa new parameter: void draw_histogram ( TH1D * h, TGraph * cross_section , double parameters , double baseline ) { std :: cout << parameters << std :: endl ; for ( int i = 1; i <= h-> GetNbinsX (); ++i) { double binCenter = h-> GetBinCenter (i); double y1 = cross_section -> Eval ( binCenter ); double y2 = h-> GetBinContent (i); // eventrate formula double S = sin (1.27/1 E3 * baseline / binCenter ); // std :: cout << i << ": " << y2 << std :: endl ; double N = (1.0 - (0.1) * S * S )* y1 * y2 * parameters ; h-> SetBinContent (i,N); // std :: cout << N << std :: endl ; } h-> GetYaxis () -> SetTitle (" Number of Events / GeV "); h-> Draw (" SAME "); std :: cout << h-> Integral (" width ") << std :: endl ; } void eventrate2 () { TFile * cross_section_file = TFile :: Open (" cross_section . root "); cross_section_file ->ls (); TGraph * cross_section = ( TGraph *) cross_section_file -> Get (" qel_cc_n ;1"); // cross_section -> Draw (); TFile * flux = TFile :: Open (" flux . root "); flux ->ls (); TH1D * h = ( TH1D *) flux -> Get (" numu_CV_AV_TPC ");
TH1D * h_sbnd = ( TH1D *) h-> Clone (" h_sbnd "); TH1D * h_microboone = ( TH1D *) h-> Clone (" h_microboone "); TH1D * h_icarus = ( TH1D *) h-> Clone (" h_icarus "); draw_histogram ( h_sbnd , cross_section ,
1E -38/1 E4 *1.6856 E30 /0.05*6.2 E20 /1 E6 *(470.0/110.0) *(470.0/110.0) , draw_histogram ( h_microboone , cross_section ,
1E -38/1 E4 *1.3395 E30 /0.05*13.2 E20 /1E6 , draw_histogram ( h_icarus , cross_section ,
1E -38/1 E4 *7.1638 E30 /0.05*6.2 E20 /1 E6 *(470.0/600) *(470.0/600) , } χ values Here, we need to introduce the parameter s2 ( sin θ ) and m (mass): void draw_histogram ( TH1D * h, TGraph * cross_section , double parameters , double baseline , double s2 , double m) { std :: cout << parameters << std :: endl ; for ( int i = 1; i <= h-> GetNbinsX (); ++i) { double binCenter = h-> GetBinCenter (i); double y1 = cross_section -> Eval ( binCenter ); double y2 = h-> GetBinContent (i); double S = sin (1.27/1 E3 * m * baseline / binCenter ); // std :: cout << i << ": " << y2 << std :: endl ; double y = y1 * y2 * parameters ; double N = (1.0 - s2 * S * S )* y; h-> SetBinContent (i,N); // std :: cout << N << std :: endl ; } h-> GetYaxis () -> SetTitle (" Number of Events / GeV "); h-> Draw (" SAME "); std :: cout << h-> Integral (" width ") << std :: endl ; } double chi2 ( TH1D * h, TH1D * g) { double chi_sq = 0.0; // for ( int i = 1; i <= h-> GetNbinsX (); ++i) for ( int i = 1; i <= 135; ++i) { double y = h-> GetBinContent (i); double lambda = g-> GetBinContent (i); if ( lambda == 0) continue ; double chi_sq_i = (y- lambda )*(y- lambda )/ lambda ; std :: cout << " << lambda << std :: endl ; chi_sq = chi_sq + chi_sq_i ; } return chi_sq ; } void eventrate_chi () { TFile * cross_section_file = TFile :: Open (" cross_section . root "); cross_section_file ->ls (); TGraph * cross_section = ( TGraph *) cross_section_file -> Get (" qel_cc_n ;1"); // cross_section -> Draw (); TFile * flux = TFile :: Open (" flux . root "); flux ->ls (); TH1D * h = ( TH1D *) flux -> Get (" numu_CV_AV_TPC ");
TH1D * h_sbnd = ( TH1D *) h-> Clone (" h_sbnd "); TH1D * h_microboone = ( TH1D *) h-> Clone (" h_microboone "); TH1D * h_icarus = ( TH1D *) h-> Clone (" h_icarus "); draw_histogram ( h_sbnd , cross_section ,
1E -38/1 E4 *1.6856 E30 /0.05*6.2 E20 /1 E6 *(470.0/110.0) *(470.0/110.0) , draw_histogram ( h_microboone , cross_section ,
1E -38/1 E4 *1.3395 E30 /0.05*13.2 E20 /1E6 , draw_histogram ( h_icarus , cross_section ,
1E -38/1 E4 *7.1638 E30 /0.05*6.2 E20 /1 E6 *(470.0/600) *(470.0/600) , TH1D * h_sbnd2 = ( TH1D *) h-> Clone (" h_sbnd2 "); TH1D * h_microboone2 = ( TH1D *) h-> Clone (" h_microboone2 "); TH1D * h_icarus2 = ( TH1D *) h-> Clone (" h_icarus2 "); draw_histogram ( h_sbnd2 , cross_section ,
1E -38/1 E4 *1.6856 E30 /0.05*6.2 E20 /1 E6 *(470.0/110.0) *(470.0/110.0) , draw_histogram ( h_microboone2 , cross_section ,
1E -38/1 E4 *1.3395 E30 /0.05*13.2 E20 /1E6 , draw_histogram ( h_icarus2 , cross_section ,
1E -38/1 E4 *7.1638 E30 /0.05*6.2 E20 /1 E6 *(470.0/600) *(470.0/600) , double chi2_sbnd = chi2 ( h_sbnd , h_sbnd2 ); double chi2_microboone = chi2 ( h_microboone , h_microboone2 ); double chi2_icarus = chi2 ( h_icarus , h_icarus2 ); std :: cout << " chi2_sbnd : " << chi2_sbnd << std :: endl ; std :: cout << " chi2_microboone : " << chi2_microboone << std :: endl ; std :: cout << " chi2_icarus : " << chi2_icarus << std :: endl ; } Armed with the event rate histograms as well as the χ values, we can write down the loop in ouroriginal function. The entire command lines are shown as below: void draw_histogram ( TH1D * h, TGraph * cross_section , double parameters , double baseline , double s2 , double m) { // std :: cout << parameters << std :: endl ; for ( int i = 1; i <= h-> GetNbinsX (); ++i) { double binCenter = h-> GetBinCenter (i); double y1 = cross_section -> Eval ( binCenter ); double y2 = h-> GetBinContent (i); double S = sin (1.27/1 E3 * m * baseline / binCenter ); // std :: cout << i << ": " << y2 << std :: endl ; double y = y1 * y2 * parameters ; double N = (1.0 - s2 * S * S )* y; h-> SetBinContent (i,N); // std :: cout << N << std :: endl ; } h-> GetYaxis () -> SetTitle (" Number of Events / GeV "); h-> Draw (" SAME "); // std :: cout << h-> Integral (" width ") << std :: endl ; } double chi2 ( TH1D * h, TH1D * g) { double chi_sq = 0.0; // for ( int i = 1; i <= h-> GetNbinsX (); ++i) for ( int i = 1; i <= 135; ++i) { double y = h-> GetBinContent (i); double lambda = g-> GetBinContent (i); if ( lambda == 0) continue ; double chi_sq_i = (y- lambda )*(y- lambda )/ lambda ; std :: cout << " << lambda << std :: endl ; chi_sq = chi_sq + chi_sq_i ; } return chi_sq ; } void eventrate_heatmap () { TFile * cross_section_file = TFile :: Open (" cross_section . root "); cross_section_file ->ls (); TGraph * cross_section = ( TGraph *) cross_section_file -> Get (" qel_cc_n ;1"); // cross_section -> Draw (); TFile * flux = TFile :: Open (" flux . root "); flux ->ls (); TH1D * h = ( TH1D *) flux -> Get (" numu_CV_AV_TPC ");
TH1D * h_sbnd = ( TH1D *) h-> Clone (" h_sbnd "); TH1D * h_microboone = ( TH1D *) h-> Clone (" h_microboone "); TH1D * h_icarus = ( TH1D *) h-> Clone (" h_icarus "); draw_histogram ( h_sbnd , cross_section ,
1E -38/1 E4 *1.6856 E30 /0.05*6.2 E20 /1 E6 *(470.0/110.0) *(470.0/110.0) , draw_histogram ( h_microboone , cross_section ,
1E -38/1 E4 *1.3395 E30 /0.05*13.2 E20 /1E6 , draw_histogram ( h_icarus , cross_section ,
1E -38/1 E4 *7.1638 E30 /0.05*6.2 E20 /1 E6 *(470.0/600) *(470.0/600) , auto heatmap = new TH2D (" heatmap "," heatmap ; log ( sin ^2 theta ); log ( Delta_m ^2) ",
100 , -3 ,0 ,
100 , -2 ,2) ; for ( int i = 1; i < heatmap -> GetNbinsX () +1; i ++) { for ( int j = 1; j < heatmap -> GetNbinsY () +1; j ++) { float x = pow (10 , heatmap -> GetXaxis () -> GetBinCenter (i)); float y = pow (10 , heatmap -> GetYaxis () -> GetBinCenter (j)); TH1D * h_sbnd2 = ( TH1D *) h-> Clone (" h_sbnd2 "); TH1D * h_microboone2 = ( TH1D *) h-> Clone (" h_microboone2 "); TH1D * h_icarus2 = ( TH1D *) h-> Clone (" h_icarus2 "); draw_histogram ( h_sbnd2 , cross_section ,
1E -38/1 E4 *1.6856 E30 /0.05*6.2 E20 /1 E6 *(470.0/110.0) *(470.0/110.0) , draw_histogram ( h_microboone2 , cross_section ,
1E -38/1 E4 *1.3395 E30 /0.05*13.2 E20 /1E6 , draw_histogram ( h_icarus2 , cross_section ,
1E -38/1 E4 *7.1638 E30 /0.05*6.2 E20 /1 E6 *(470.0/600) *(470.0/600) , double chi2_sbnd = chi2 ( h_sbnd , h_sbnd2 ); double chi2_microboone = chi2 ( h_microboone , h_microboone2 ); double chi2_icarus = chi2 ( h_icarus , h_icarus2 ); double chi2_total = chi2_sbnd + chi2_microboone + chi2_icarus ; heatmap -> SetBinContent (i,j, chi2_total ); std :: cout << " chi2_sbnd : " << chi2_sbnd << std :: endl ; std :: cout << " chi2_microboone : " << chi2_microboone << std :: endl ; std :: cout << " chi2_icarus : " << chi2_icarus << std :: endl ; std :: cout << "x" << x << std :: endl ; std :: cout << "y" << y << std :: endl ; delete h_sbnd2 ; delete h_microboone2 ; delete h_icarus2 ; } } heatmap -> Draw (" colz "); gPad -> SetLogz (); }}