1. Part B: Here you will apply what you learned to compute the volume of an n-dimensional sphere, and its error, essentially computing integrals using the MC acceptance/rejection method.

• Bevington & Robison: Testing the Fit, you will be fitting your data in part three to a Gaussian function for a discussion on this topic check chapters 9-11

Part B. Calculate the Volume of N-dimensional spheres

Part C: Generate Poisson deviates, histogram, and plot the distribution

1. Part A: Here you will first employ the Monte Carlo acceptance/rejection method to compute the value of π and its statistical uncertainty. You may use any programming language you want as long as it is C, C++, Python, MATLAB, or FORTRAN. Since all have deployed an IDE and have been encouraged to use it to develop python scripts and Dr. Boeglin's excellent set of python wrapper tools in his <a data-wikiword="LabTools3" href="http://wanda.fiu.edu/boeglinw/LabTools3/" target="_blank" title="LabTools3"> LabTools3 </a> encourage you to continue to use that for this and all of the computing projects assigned in this class.

• Bevington & Robison: Probability and Distributions, Binomial, Poisson, and Gaussian distributions all of chapter 2 pgs 17-31.

GradingRubric

1. Part A: Here you will first employ the Monte Carlo acceptance/rejection method to compute the value of π and its statistical uncertainty. You may use any programming language you want as long as it is C, C++, Python, MATLAB, or FORTRAN. Since all have deployed an IDE and have been encouraged to use it to develop python scripts and Dr. Boeglin's excellent set of python wrapper tools in his LabTools3 encourage you to continue to use that for this and all of the computing projects assigned in this class.

Part C: Generate Poisson deviates, histogram and plot the distribution

In Part C the assignment is to generate a sample of Poisson deviates using Monte Carlo methods, specifically the transformation method. The method depends on the conservation of probabilities which stipulates that a draw from a uniform random distribution is equivalent to a draw from any other distribution, including the Poisson. Symbolically,

where is the uniform distribution defined on the interval and is, in this case, the Poisson distribution function, a discrete distribution,

with x an integer and mean μ.

Here is how to do it:

1. draw a random number r from a uniform distribution by calling the random() function
2. now compute the Poisson probability, the argument of the summation above, for integers x starting at 0 and continuing on to some n. Each time add terms to form the summation as required by the equality above.
3. Compare the summation as the integers increase and compare to the random number r drawn from the uniform distribution above.
4. Once you have a summation that is greater than the r, which corresponds to the LHS, write the integer in the previous iteration x into a table and that will be your Poisson deviate that corresponds to r,

Repeat this process N = 10,000 times to obtain 10,000 Poisson deviates.

You will need as an input parameter the mean of the distribution. You will then histogram the deviates and plot the histograms, include these in your report. Remember that histogram binning is important. If you don't get this right the plots will look funny, either there will be spaces between bins of some bins will be much larger than their neighbors, indicating overflow. Given that you are plotting Poisson distributions that are discrete the bin width is = 1.Please repeat the generation of Poisson deviates 3 times: one with

• Bevington & Robison: Probability and Distributions, Binomial, Poisson, and Gaussian distributions all of chapter 2 pgs 17-31.
• Bevington & Robison: Monte Carlo Techniques, The acceptance-rejection method you'll use, and most aspects of this lab are spelled out in chapter 5 pgs 75-94.

• The computer program can be written in any language you want as long as it is FORTRAN, C, C++, Python, or MATLAB.
• If you write a C, C++, or FORTRAN program you will need to compile your code. You do this by issuing the appropriate compiler command "g++" for a C or C++ program, "gfortran" for FORTRAN.
• An outline of the procedure is described below. The exact prescription is detailed in my places, online as well as in Bevington & Robison. Basically, the idea is to compare the area of a circle to the area of a square. This ratio is proportional to π and we use this fact to get π from an MC estimate of the ratio. The pseudo-code that will accomplish this is:

1. 1. Next check if x,y point falls within a circle by the condition if (x**2+y**2 < 1) then increment a counter (integer) by one.
   2. Compute π by taking the ratio of your counter to the number of iterations. As this is proportional to π you have completed the computation of the central value of π.
   3. Finally, calculate the statistical uncertainty of π given that the process used to determine π is a binomial one the uncertainty will thus be proportional but not equal to the binomial uncertainty.

• Also part of the assignment is to find the number of events needed to obtain statistical uncertainties of 1%, 0.01%, and 0.00001%. One way to do this is the brute force way and use the program changing the number of iterations until you get the desired result. This could take a while at least for the smallest uncertainties. There is an analytic method that involves solving for N using the value of pi from your calculator. Make sure you include these estimates in your paper.

Part B. Calculate the Volume of an N-dimensional spheres

Modify your program, or write a new one, to determine the volumes of N-dimensional spheres (N=1,2,3,4,5 but you've already done the 2D one) with radius=1.0. Use 10,000 tries and make sure you quote your uncertainties. You should also compare your results with analytic values that you can compute by some other means. Please include the means used to compute the analytic values.

Part D. Generate and Fit a Gaussian

In this exercise, you employ the Monte Carlo technique to generate data distributed as a normal or Gaussian distribution. Then you will fit the generated data sample to a Gaussian Probability Distribution Function (PDF), determine the fit parameters and their errors, and discuss the goodness of the fit via the reduced chi2.

There are many ways to generate a normally distributed set of data. One way is to start off with a uniform (flat) distribution of data points, generated between 0 to 1, and transform them to a normally distributed set of data via the Box-Muller 2D Gaussian transformation described in Bevington and Robinson or here<a href:http://pdg.lbl.gov/2011/reviews/rpp2011-rev-monte-carlo-techniques.pdf > http://pdg.lbl.gov/2011/reviews/rpp2011-rev-monte-carlo-techniques.pdf </a> This clever algorithm is both efficient and easy to implement. There are many sources you can read about it, there is even a Wikipedia article on it. There is an easier way but one I do not want you to use. It turns out that given the usefulness of having random numbers drawn from a Gaussian are, every programming language has normal or Gaussian deviate generators in standard libraries. For example, numerical Python has a built-in normal random number generator in NumPy and probably all other math libraries. Again please do not use a built-in method or any functions that invoke a normal random number generator. The assignment calls for you "coding up" the Box-Muller procedure yourself.

Here is how to do it:

• Use uniform random number generators to generate random numbers drawn from the Gaussian centered at 0 with a variance of 1.0 using the Box-Muller transform as described in the references above.
• Bin the transformed random number data into 40 bins from -4.0 to + 4.0. 1.0.
• Fit the histogram generated above with a Gaussian function, or Gaussian PDF.

• Find the parameters of the fit and quote their uncertainties. Do you know how these were determined? Feel free to explore how fitting package work in general, and how they find fit parameters. Quote the parameters and their uncertainties prominently in your lab report. You should discuss how these were determined by the fitter. You should look through Bevington and Robinson for background and feel free to explore how the functions used in Python actually do the fits (this may take a bit of digging around Dr. Boeglin's LT.box code since those tools wrap standard Python curve fitting libraries found in scipy and maybe matplotlib. Learning about this is something everyone do even if it is done at a superficial level.

• Determine the goodness of fit. You can use the reduced chi2 if your fitter uses the least-squares method to determine the fit parameters. Also, discuss in your lab report the use of the reduced chi2, what it means and why can we use it as a means by which to estimate the "goodness of fit".

GradingRubric

1. You will first use the Monte Carlo acceptance/rejection method to compute the value of π and its statistical uncertainty. You may any programming language listed below: C, C++, python, matlab, or FORTRAN.

• The computer program can be written in any language you want as long as it is FORTRAN, C, C++, python, or Matlab.

• Bevington & Robison: Probablity and Distributions, Binomial, Poisson and Gaussian distributions all of chapter 2 pgs 17-31.

• An outline of the procedure is described below. The exact perscription is detailed in my places, online as well as in Bevinton & Robison. The basically the idea is to compare the area of a circle to the area of a square. This ratio is proportional to π and we use this fact to get π from a MC estimate of the ratio. The pseudo code that will accomplish this is:

1. 1. Compute π by taking the ratio of your counter to the number of iterations. As this is proportional to π you have completed the computation of the central value of %pi;.
   2. Finally calculate the statistical uncertainty of π given that the process used to determine π is a binomial one the uncertainty will thus be proportional but not equal to the binomial uncertainty.

• Also part of the assignment is to find the number of events needed to obtain statistical uncertainties of 1%, 0.01% and 0.00001%. One way to do this is the brute force way and use the program changing the number of iterations until you get the desired result. This could take a while at least for the smallest uncertainties. There is analytic method that involves solving for N using the value of pi from your calculator. Make sure you include these estimates in your paper.

• Generate data distributed normally, start with rand() and transform to Gaussian or use the Box-Muller transforms.

GradingRubric

• You can create your program either on your PC or Mac.
• The computer program can be written in any language you want as long as it is FORTRAN, C, C++ or python.

Suggested Reading for this lab:

• You can create your program either on your PC or Mac or on the Linux server medianoche.hep.fiu.edu. An account on medianoche was or will be created for you during class.
• The computer program can be written in any language you want as long as it is FORTRAN, C, C++ or python. I've placed "starter" code files in each account's home directory to get you started. The files are called sample.cc, sample.py etc. Please not that account content will be deleted as soon as the course ends. So please move off anything you want to keep before then. Even if you decide to write your code on your PC/Mac you will still be required to copy the code over to the Linux server and make sure it runs there. Please note that the python environment is different on the Linux server so your python code may not function properly. If so see me and we'll work to resolve this individually. To gain access to the Linux server see the page HowToAccessComputing.
• If you write a C, C++ or FORTRAN program you will need to compile your code. You do this by issuing the appropriate compiler command "g++" for a C or C++ program, "gfortran" for FORTRAN.
• An outline of the procedure for using MC methods to compute the value of pi is described below and available from many sources on the web. Basically the idea is to compare the area of a circle to the area of a square. Since this ratio is proportional to &pi and can be easily computed its value numerically via this MC method. Here is pseudo code that will accomplish this:
  1. Throw (generate) two random numbers "x" and "y" with a random number generator. In c++ for example this can be done by calling the function rand(). It is available in the standard C++ library stndlib. Make sure the number returned in the interval of [0,1].
  2. Next check if x,y point falls within a circle by the condition if (x**2+y**2 < 1) then increment an counter (integer) by one.
  3. Compute &pi by taking the ratio your counter to the number of iterations since this is proportional to pi you are done.
  4. Finally calculate the statistical uncertainty of &pi given that the process to find is binomial then uncertainty in pi is proportional to the binomial uncertainty.
• Also part of the assignment is to find the number of events needed to obtain statistical uncertainties of 1%, 0.01% and 0.00001%. One way to do this is the brute force way and use the program changing the number of iterations until you get the desired result. This could take a while at least for the smallest uncertainties. There is analytic method that involves solving for N using the value of pi from your calculator. Make sure you include these estimates in your paper.

There are many ways to generate a normally distributed set of data. One way is start off with a uniform (flat) distribution of data points, generated between 0 to 1 and transform them to a normally distributed set of data via the Fox-Muller transform, described here http://pdg.lbl.gov/2011/reviews/rpp2011-rev-monte-carlo-techniques.pdf. This method is probably the easiest to understand but not the simplest to implement. That would be to use the builtin normal random number generators available from libraries that come with your programming environment. For example numerical python has builtin normal random number generator and so does root. I am not sure if the stnd libraries in C/C++ do but I am sure they exist. I would prefer you use the easier to understand method using rand() but the you can use the simpler ones too.

Here is how to proceed:

• Generate data distributed normally, start with rand() and transform to Gaussian or use the functions in numpy or root.
• Bin the data into 40 bin from -4.0 to + 4.0, if you chose your normally distributed data to have a mean of 0 and width of 1.0, the default.
• Fit the histogram generated above with a Gaussian PDF
• Find the parameters of the fit and quote their uncertainties. Do you know how these were determined? Feel free to explore how your fitting package works to determine these parameters. Also quote their uncertainties and discuss how these were determined by the fitter. For useful background refer to Bevington and Robinson. Learning about this is something everyone do even if it is done at a superficial level.
• Determine the goodness of fit. You can use the reduced chi2 if your fitter uses the least squares method to determine the fit parameters. Also discuss what value of the reduced chi2 means and write in your report if your chi2/dof is reasonable or not.

GradingRubric

Suggested Reading for this lab:

• Bevington & Robison: Probablity and Distributions, Binomial, Poisson and Gaussian distributions all of chapter 2 pgs 17-31.
• Bevington & Robison: Monte Carlo Techniques, The acceptance rejection method you'll use and most aspects of this lab are spelled out in chapter 5 pgs 75-94.
• Bevington & Robison: Testing the Fit, you will be fitting your data in part three to a Gaussian function for a discussion on this topic check chapter 9-11

• The computer program can be written in any language you want: (FORTRAN, C/C++, java, python...). I've placed code in each account to get you started. There is sample.* file for c++ and python. If you want to use FORTRAN I can provide you with a sample.F file too. At the beginning of the semester you were given an account on one of our HEP analysis machines, "medianoche.hep.fiu.edu" please feel free to use that server as you wish for the remainder of the semester. Accounts are typically deleted as soon as you complete the course. I encourage you to make use of the linux server, it will come in handy later in your studies or work life. The server is configured with a gnu C/C++ and FORTRAN compiler plus one or more python interpreter. Instructions on how to access our server from either a PC or a Mac are available in the HowToAccessComputing page. To compile C++/FORTRAN code issue the command "g++" of "gfortran" from the terminal window. As I informed you, you may use your own computer to work instead but you'll need to deposite your raw code on the medianoche server for me to inspect the code and to see if your program runs and produces the proper output, it dosen't have to be pretty but it must work. If you'd like you can insert your code in your paper as an appendix.
• An outline of the procedure for using MC methods to compute the value of pi is described below and available from many sources on the web. Basically the idea is to compare the area of a circle to the area of a square whose. This ratio is simply the value of pi. The area is constructed by generating a pair of random numbers representing x and y pairs and then simply counting the number of times these randomly generated points fall within a circle and comparing it to the number which fall in a square, just the number of tries. Here is pseudo code that will accomplish it:
  1. Throw (generate) two random numbers with a random number generator. In c++ for example this can be done by calling the function rand(). Its in one of the standard C++ libraries stndlib. Make sure the number returned in the interval of [0,1].

• You need to do a proper treatment of uncertainties in this lab. First consider where are the uncertainties in this process, which come about from the acceptance/rejection method which we use. This is just a binomial process. As mentioned in the reading material above Bevington and Robinson have a nice writeup on binomial distributions also Chapter 10 in Melissinos and Napolitano discusses this. But many sources are available on the web. You should calculate this statistical errors in your program and print out the result. That way you can get directly from the program no matter how many iterations use.
• Part of the assignment is to find the number of events needed to obtain statistical uncertainties of order 1%, 0.01% and 0.00001%. One way to do this it to use the program you wrote to estimate the number of events. However, you may find it difficult to do this for a relative error of 0.00001%. There is fact an analytic way to do it (Hint: for this step you can use the known value of pi).
• Also, elaborate on sources of systematic uncertainties if any. Systematic uncertainties result from sources of errors that do not arise from random fluctuations. In this experiment I can think of only one source of systematic, namely the random number generator used. A way to test this is to change the random number generator, by either creating your own or finding an alternative function to call. Then repeat the experiment and see if this significantly changes your result. Describe what you observe and quote the discrepancy as an additional uncertainty. Systematic uncertainties such as these are added in quadrature with statistical uncertainties.

• Write a program that uses the acceptance/rejection Monte Carlo Method to calculate the value of π. Use your program to find the value to an accuracy of 0.1% and quote that as your main result of part A.
• The computer program can be written in any language you want: (FORTRAN, C/C++, java, python...). There are many sources on the web that will guide you on how to do this complete with code snippets etc., At the beginning of the semester you were given an account on one of our HEP analysis machines, "medianoche.hep.fiu.edu" please logon that guy to write, compile and run your program. The server has a gnu C/C++ and FORTRAN compiler plus pretty much any thing else you may need. There are instructions on how to access the server if you are using a PC or Mac here HowToAccessComputing. To compile C++/FORTRAN code issue the command "g++" or "gfortran" from a terminal on the server. If you insist on using your own computer to do your work please consult with me first. You will be required to turn in your code as part of the write up. Please include it as an appendix with your lab report.
• An outline of the procedure for using MC methods to compute the value of pi is described below and available from many sources on the web. Basically though the idea is to compare the area of a circle to the area of a square whose ratio is simply the value of pi. The area is formed by generating random numbers in two dimensions, ie., x and y separately and counting the number of times the randomly generated x and y number falls within the possible x y space. Here is some pseudo code that will accomplish this:

1. 1. Increment an integer, initial set to 0, each time you loop through your block of code increasing the counter if and only if the x or y point falls within the square.

• You will also be required to do a proper treatment of uncertainties in this lab. In particular the statistical uncertainties associated with Monte Carlo processes are, ie., the acceptance/rejection method being used is a binomial process. There are many good references on the web on statistics on the web, Wikipedia for example, but I find Bevington and Robinson and/or Chapter 10 in Melissinos and Napolitano particularly good introductory sources. The statistical errors should included in your program as mentioned in step 5. of the pseudo code above.
• Part of the assignment is to find the number of events needed to obtain statistical errors of order 1%, 0.01% and 0.00001%. One way to do this it to use the program you wrote to estimate the number of events. However, you may find it difficult to do this for the relative error of 0.00001%. There is fact an analytic way to do it (Hint: for this step you can use the known value of pi).
• Also, elaborate on sources of systematic uncertainties if any. Systematic uncertainties result from sources of errors that are not statistical in nature. In this experiment I can think of only one source of systematic, namely the random number generator used. A way to test this is to change the random number generator, by either creating your own or finding an alternative function to call from C++. Then repeat the experiment and see if this significantly changes your result. Describe what you observe and quote the discrepancy as an additional uncertainty. Systematic uncertainties such as these are added in quadrature with statistical uncertainties.

• In this exercise you employ the Monte Carlo technique to generate a normal or Gaussian distribution. Then you will fit this distribution to a Gaussian Probability Distribution Function (PDF). This is another common use of the Monte Carlo method, especially in particle and nuclear physics. Basically, we take the uniform (flat) distribution from 0 to 1 provided to you by the built in function rand() and use a transform to convert that into a distribution of random numbers now distributed as a Gaussian. In fact, this use of the Monte Carlo method can be used to generate any kind of distribution but some are more difficult than other, see Bevington and Robinson for a useful description of the use of transforms in Monte Carlo techniques.

• To generate your Gaussian data distribution you will again use a random number generator over interval 0 to 1. The set of number, uniformly distributed from 0 to 1 will be "transformed" into a new set of number now distributed as Gaussian. There are different algorithms available to do the transformation from the central limit theorem to transforms or algorithms such as the Box-Muller transform. The Particle Data Group Math Review section (http://pdg.lbl.gov/2011/reviews/rpp2011-rev-monte-carlo-techniques.pdf) describes a couple of ways. If you want to research alternatives please feel free to do so.
• Once written the code to generate the distribution of Gaussian data I would like you to histogram into 40 bins from -4.0 to 4.0. The histograming can in principle b part of your code or you can use some other histograming package to do it for you. One way is to call hbook histograming functions from the FORTRAN libraries available on the server or you can cut and paste the output into excel. If you are using python I am pretty sure there are histograming functions or methods available as well. There are many other ways too. Feel free to consult with your instructor. I would recommend though to do it yourself, gives you a little more practice programming...

• Fit your distribution to a Gaussian Probability Distribution Function (PDF). As stated earlier you can use whatever package you want to do this but I recommend you use the python packages you used in Modern Lab or MNFIT, an older but more sophisticated analysis package similar to what we use in HEP root. If you are familiar with root you can use that but the learning curve is rather steep so instead I'll describe briefly below how to use the MNFIT data analysis package which is as powerful as root but easier, in my opinion, to use. .
• To use mn_fit to fit your histogram converted the histogram into an ASCII formatted file that can be read by mn_fit. There are example formatted ASCII histograms in a file on the medianoche server in /home/jorge/PHY4810L/dat_test.mnd. Information on how to access the medianoche server are detailed in here The format of the file includes an entry for each bin, each row corresponds to one bin, the value of X (the centroid of the bin) and the value of Y (the number of entries per bin). In the dat_test.mnd file there are three examples on how you can enter data. Histogram 1 has errors in Y (the sqrt of the number of entries), the 2 histogram specifies no errors (mnfit uses the sqrt(Y) for the error) and Histogram 3 has errors for each variable which are allowed to be asymmetric. For this lab use histograms of type 2 or 1 if you want to calculate your own errors.
• MN_FIT has been installed on the medianoche.hep.fiu.edu server. To access medianoche please check the HowToAccessComputing web

• • You are now on the medianoche server indicated by [jorge@medianoche ~/]$ prompt

    [jorge@medianoche ~/]$ mn_fit 

    hit return or type 2 to display the X11 graphics window

    MN_CMD> dat_fetch histogramDATA.mnd 0 

    Now fetch the name of your ASCII or hbook histogram (there is a sample data .mnd file attached below)

    MN_CMD>index 

These commands will first list all histograms read into mnfit then plot histogram 10&0 (the &0 indicates a secondary identifier) and add a Gaussian function as the probability distribution function to be used in a fit to your data. The fun add Gaussian will expect arguments, use sigma type, and enter reasonable values for the parameters when prompted. By reasonable I mean values close to what you see for your plot.

Now do the following:

MN_CMD> fit 1

The fit 1 command will fit the histogram number 1 with the function loaded into the fitter. You will have to select what type of fit to use. You can use the chi2 fit for this exercise. The "minimize" command will put you into the minimization and error analysis package called MINUIT. MINUIT is a set of functions algorithms and a complete package from CERN. Its been around for a very long time so has lots of features and is used by many people not only in high energy/nuclear physics. It may be of interest to you to spend some time exploring the documentation from CERN. The "display" command will displays your histogram with the fitted function (the function with the best parameters obtained by the fitting procedure overlay on the plot complete with statistical errors. The numbers are also displayed in tabular form below. The command "minos" will recalculate the errors using a more sophisticated approach where the minimizer doesn't assume that your errors are parabolic. For more information check the CERN documentation

Please note that while you will work closely with your lab partner you are expected to write your own code and present to me on the server as part of your grade.

ALSO: Include in your report the number of events needed to get a statistical uncertainty of 1%, 0.1% and 0.00001%.

Part C. Generate and Fit a Gaussian

Using Monte Carlo Methods generate a data distributed as a Gaussian, bin the data into a histogram and then fit the histogram with a Gaussian probability distribution function (PDF).

• Write a program that uses the acceptance/rejection Monte Carlo Method to calculate the value of π. Use your program to find the value to an accruacy of 0.1% and quote that as your main result of part A.

• You will also be required to do a proper treatment of uncertainties in this lab. In particular the statistical uncertainties associated with Monte Carlo processes are, ie., the acceptance/rejection method being used is a binomial process. There are many good references on the web on statistics on the web, wikipedia for example, but I find Bevinton and Robinson and/or Chapter 10 in Mellission and Napolitano particularly good introductory sources. The statistical errors should included in your program as mentioned in step 5. of the pseudo code above.
• Part of the assignment is to find the number of events needed to obtain statistical errors of order 1%, 0.01% and 0.00001%. One way to do this it to use the program you wrote to estimate the number of events. However, you may find it difficult to do this for the relative error of 0.00001%. There is fact an analytical way to do it (Hint: for this step you can use the known value of pi).
• Also, elaborate on sources of systematic uncertainties if any. Systematic uncertainties result from sources of errors that are not statistical in nature. In this experiment I can think of only one source of systematic, namely the random number generator used. A way to test this is to change the random number generator, by either creating your own or finding an alternative function to call from C++. Then repeat the experiment and see if this significantly changes your result. Describe what you observe and quote the discrepency as an additional uncertainty. Systematic uncertainties such as these are added in quadrature with statistical uncertainties.

Assignment: Use Monte Carlo Methods to determine the volume of an N-dimensional sphere of unit radius when N=1,2,3,4 and 5.

• In this exercise you employ the Monte Carlo technique to generate a normal or Gaussian distribution. Then you will fit this distribution to a Gaussian Probability Distribution Function (PDF). This is another common use of the Monte Carlo method, especially in particle and nuclear physics. Basically, we take the uniform (flat) distribution from 0 to 1 provided to you by the built in function rand() and use a transform to convert that into a distribution of random numbers now distributed as a Gaussian. In fact, this use of the Monte Carlo method can be used to generate any kind of distribution but some are mored difficult than other, see Bevington and Robinson for a useful description of the use of transforms in Monte Carlo techniques.
• In this lab you will first create the aforementioned set of data, Then use a fitting package, (an application that allows you to plot data, fit functions to data sets and display them and the fit parameters for publication). In an real experiments, particularly when the physical processes are difficult to extract or the data analysis is particularly sophisticated, the ability to simulate the reaction of the experiment to simulated data is a critical part of the data analysis.
• To generate your Gaussian data distribution you will again use a random number generator over interval 0 to 1. The set of number, uniformly distributed from 0 to 1 will be "transformed" into a new set of number now distributed as Gaussian. There are different algorithms available to do the transformation from the central limit theorm to transforms or algorithms such as the Box-Muller transform. The Particle Data Group Math Review section (http://pdg.lbl.gov/2011/reviews/rpp2011-rev-monte-carlo-techniques.pdf) describes a couple of ways. If you want to research alternatives please feel free to do so.
• Once written the code to generate the distribution of Gaussian data I would like you to histogram into 40 bins from -4.0 to 4.0. The histogramming can in principle b part of your code or you can use some other historgraming package to do it for you. One way is to call hbook histogramming functions from the FORTRAN libraries available on the server or you can cut and paste the output into excel. If you are using python I am pretty sure there are histograming functions or methods available as well. There are many other ways too. Feel free to consult with your instructor. I would recommend though to do it yourself, gives you a little more practice programming...
• Once you have your histogram you will next need to fit the simulated data with whatever fitting package you choose. Dr. Boeglin's Modern Physics Lab provided you with an interface to the python fitting libraries that will be more than suitable for this exercise. In fact, the entire lab can be done in python scripts. Dr. Boeglin has done a nice job of documenting procedures on how to use the very powerful python libraries. All of Dr. Boeglin's software is available on the medianoche server as well. Additionally, I make available to you a particularly simple but powerful package, called MINFIT, that I used in as a graduate student and served me quite well. This package or application is also installed on the medianoche server. There are instructions below on how to use it and an online manual is available.

• Write a program that uses the acceptance/rejection Monte Carlo Method to calculate the value of π.
• The computer program can be writen in any language you want: (FORTRAN, C/C++, java, python...). There are many sources on the web that will guide you on how to do this complete with code snipets etc., At the begining of the semester you were given an account on one of our HEP analysis machines, "medianoche.hep.fiu.edu" please logon that guy to write, compile and run your program. The server has a gnu C/C++ and FORTRAN compiler plus pretty much any thing else you may need. There are instructions on how to access the server if you are using a PC or Mac here HowToAccessComputing. To compile c++/fortran code issue the command "g++" or "gfortran" from a terminal on the server. If you insist on using your own computer to do your work please consult with me first. You will be required to turn in your code as part of the write up. Please include it as an appendix with your lab report.
• An outline of the procedure for using MC methods to compute the value of pi is described below and available from many sources on the web. Basically though the idea is to comapre the area of a circle to the area of a square whose ratio is simply the value of pi. The area is formed by generating random numbers in two dimensions, ie., x and y seperately and counting the number of times the randomly generated x and y number falls within the possible x y space. Here is some psuedo code that will accomplish this:

• Part of the assignment is to find the number of events needed to obtain statistical errors of order 1%, 0.01% and 0.00001%. You can use your program to estimate the number of events experimentally but you may find it difficult to do that with the 0.00001% guy. Feel free to discuss this with your instructor. There may be an easier way

• • You are now on the medianoche server indicated by [jorge@medianoche ~/]$ prompt

    [jorge@medianoche ~/]$ mn_fit 

    hit return or type 2 to display the X11 graphics window

    MN_CMD> dat_fetch histogramDATA.mnd 0 

    Now fetch the name of your ascii or hbook histogram (there is a sample data .mnd file attached below)

    MN_CMD>index 
    MN_CMD>plot 10&0 
    MN_CMD>fun add gaussian 

    These commands will first list all histograms read into mnfit then plot higtogram 10&0 (the &0 indicates a secondary identifier) and add a Gaussian function as the probability distribution function to be used in a fit to your data. The fun add Gaussian will expect arguments, use sigma type, and enter reasonable values for the parameters when prompted. By reasonable I mean values close to what you see for your plot.

    Now do the following:

    MN_CMD> fit 1
    MINUIT> minimize 
    MINUIT> display 
    MINUIT> minos 

    The fit 1 command will fit the histogram number 1 with the function loaded into the fitter. You will have to select what type of fit to use. You can use the chi2 fit for this exercise. The "minimize" command will put you into the minimization and error analysis package called MINUIT. MINUIT is a set of functions algorithms and a complete package from CERN. Its been around for a very long time so has lots of features and is used by many people not only in high energy/nuclear physics. It may be of interest to you to spend some time exploring the documentation from CERN. The "display" command will displays your histogram with the fitted function (the function with the best parameters obtained by the fitting procedure overlayed on the plot complete with statistical errors. The numbers are also displayed in tabular form below. The command "minos" will recalculate the errors using a more sophisticated approach where the minimizer doesn't assume that your errors are parabolic. For more information check the CERN documentation

Lab Assignment: Monte Carlo Techniques

This is the first lab that you will complete. Please prepare a lab report written in the form of a Journal of Physics Article that has the form illustrated in Dr. Raue's paper linked to from PHY4810LSyllabus page. This lab report should be done independently, ie everyone will turn one in. For the lab report you may consult each other and myself but are expected to work somewhat independently, ie prepare your own program and turn in a unique report. Site all references properly.

Part A. Calculate the value of π (pi) using Monte Carlo methods

Assignment: Calculate the value of π and find the number of events needed to get a statistical uncertainty of 1%, 0.1% and 0.00001%.

• Use Acceptance/Rejection Monte Carlo Method to calculate the value of π.
• Do this by writing a computer program, in any language you want: (FORTRAN, C/C++, java, python...). There are many sources on the web that will guide you on how to do this complete with code snipets etc., You can use our HEP analysis machine, "medianoche.hep.fiu.edu" which has a functioning gnu C/C++ and FORTRAN compiler. To compile c++/fortran code issue the command "g++" or "gfortran" from a terminal on the server. Java, python and variety of other scripting languages common on the Linux platform are also available on medianoche. Of course you can use your own computer if you wish. Free programing languages are available for windows based PCs and the MAC already has a c++/c compiler and scripting languages such as python. I will require you to turn in your code as part of the write up. Please include it as an appendix with your lab report.
• An outline of the procedure for using MC methods to compute the value of pi is described below but the basic idea is to comapre the area of a circle to the area of a square whose ratio is simply the value of pi. The area is formed by generating random numbers in two dimensions, ie., x and y seperately and counting the number of times the randomly generated x and y number falls within the possible x y space. Here is some psuedo code that will accomplish this:
  1. Throw (generate) two random numbers with your favorite random number generator. This is done by calling a random number generating function or subroutine from the standard libraries (for example in C++ the function is called rand()) in the interval of [0,1].
  2. Increment an integer, initialy set to 0, each time you loop through your block of code increasing the counter if and only if the x or y point falls within the square.
  3. Next check if x,y point falls within a circle. If it does increment another counter, also set to 0 initially, by one.
  4. Now calculate PI by simply taking the ratio of the number of times (counts) in the square (your first counter) vs. number of counts in circle (your second counter). This is nothing more than the ratio of areas of your circle to your square.
  5. Calculate the statistical error for the binomial statistical process that results from the acceptance/rejection method.
• You will also be required to do a proper treatment of errors in this lab. In particular the statistical errors associated with the Monte Carlo need to be included in your report. The aforementioned MC process which uses a acceptance/rejection method is inherently a binomial process. You will need to research this to learn how to incorperate this into your results. A good source is in Chapter 10 of Mellissinos but many good sources are available on the web. In step 5. of the pseudo code above you will need to compute the statistical error as you your program runs through the counts.
• Part of the assignment is to find the number of events needed to obtain statistical errors of order 1%, 0.01% and 0.00001%.
• Also, elaborate on sources of systematic uncertainties if any. Systematic uncertainties result from sources of errors that are not statistical in nature. In this experiment I can think of only one source of systematic, namely the random number generator used. A way to test this is to change the random number generator or the way it is envoked and then examine the difference in the results. Describe what you observe and quote an additional error even if its smaller than the statistical error. Please note that systematic uncertainties are added in quadrature with statistical uncertainties.

Part B. Calculate the Volume of an N-dimensional sphere (N=1-5)

Assignment: Use Monte Carlo Methods to determine the volume of N-dimensional sphere of unit radius from N=1 to N=5.

For this lab use 10,000 tries to determine the central value f the volume. In your results quote that and the associated statistical error. Also compare your result with theoretical value determined analytically.

C. Generate and Fit a Gaussian

Assignment: Using Monte Carlo Methods generate a data distributed as a Gaussian, bin the data into a histogram and then fit the histogram with a Gaussian probability distribution function (PDF).

In this exercise you will be using a random number generator to create a sample distribution which will then be fitted to a PDF of the known form. This is another common use of the Monte Carlo method, in data analysis whereby distributions are created to simulate a phenomena of interest. In this lab you will first create the data and then use a fitting package or application to fit a known PDF to the data generated. Typically, in an experiment the data will come from the experiment instead of a computer but the data processing, data fitting and parameter estimation procedure can be learned from the simulated data we'll create here.

There are essentially two parts to this lab. First you will generate and then process the a set of simulated data which when histogrammed will be distributed as Gaussian. To generate your simulated data distribution you will once again use a random number generator used previously to create a uniform distribution from 0 to 1. The uniformly distributed output will be applied in an algorithm to create the Gaussian distributed data. There are many different algorithms available to do this and a variety of sources exist to aid you, for example: there are a couple of ways described in the Particle Data Group Math Review section. Alternatives that are faster more or less accurate exists and you are encouraged to explore on the on the web.

Once you've coded the MC algorithm you will need to histogram the output, remember we are ultimately going to fit the histogram to a Gaussian PDF as part of this exercise. The histogramming process can be done in the code itself or you can feed the raw output data to histogramming package. In High Energy Particle physics we use HBOOK to histogram or ROOT or a variety of other applications. In fact, MS excel also has a histogramming feature built into the data analysis package. The easiest thing to do is to create your histogram directly in your Monte Carlo code. Once you have your histogram you will need to fit the simulated data you will need to fit with a fitting package. In principle you could write your own fitting package but this would be a considerable amount of work. Your best bet is to use an application with fitting built in. In this class we have made available a package that is relatively easy to use and very powerful. The package or application is deployed on the HEP analysis cluster and is called mn_fit. There are instructions below on how to use it and a online manual is available. Alternatively you can use the fitting packages in MATLAB, or Matehmatica if you so desired. I believe even MS excel can be coxed to do curve fitting, at least a simple Gaussian fit as we will do.

Some helpful hints on how to proceed Generating the simulated data distributed as a Gaussian

• Throw random numbers between [0,1] and use the algorithm defined in say the PDG link above to create your Gaussian distribution. I would generate at least 1000 to 10,000 sample points, enough so that each bin below has no less than 5 entries. If you go too far the fitting procedure will have difficulty with dealing with small errors.
• Bin the output variable (your independent variable (x) ) in a histogram with 40 bins. Given that your Gaussian data will likely have a mean of 0 and width of 1.0 I would range your independent variable between -4.0 and +4.0 As mentioned earlier binning data can be done within your MC program by simply counting the number of times the output dependent variable falls within a range of values defined by your bin size. For example if you use 40 bins and your (x) range is from -5.0 to +5.0 then the bin size is 10/40 or .25. The independent variable (plotted on the horizontal axis) is then the centroid of your bin and the dependent variable (plotted on the vertical axis) is the number of entries that fall within the value of x in a particular bin.
• Fit your distribution to a Gaussian Probability Distribution Function (PDF). As stated earlier you can use whatever package you want to do this, MATLAB, Mathematica, Mapel, ROOT, PAW to name a few. In what follows I will describe how to use my favorite fitting package called mn_fit.
• To use mn_fit to fit your histogram converted the histogram into an ASCII formated file that can be read by mn_fit. There are example formatted ascii histograms in a file on the medianoche server in /home/jorge/PHY4810L/dat_test.mnd. Information on how to access the medianoche server are detailed in hereThe format of the file includes an entry for each bin, each row corresponds to one bin, the value of X (the centroid of the bin) and the value of Y (the number of entries per bin). In the dat_test.mnd file there are three examples on how you can enter data. Histogram 1 has errors in Y (the sqrt of the number of entries), the 2 histogram specifies no errors (mnfit uses the sqrt(Y) for the error) and Histogram 3 has errors for each variable which are allowed to be asymmetric. For this lab use histograms of type 2 or 1 if you want to calculate your own errors.
• MN_FIT has been installed on the medianoche.hep.fiu.edu server. To access medianoch please check the HowToAccessComputing web
• Here is sample session on how to get onto medianoche and use mn_fit to fit your histogram:
  • First you will need to login to into medianoche via a ssh (secure shell) either through putty or the ssh client on your Mac or linux box. (check HowToAccessComputing).
  • You are now on the medianoche server indicated by [jorge@medianoche ~/]$ prompt

[jorge@medianoche ~/]$ mn_fit // hit return or type 2 to display the X11 graphics window
MN_CMD> dat_fetch histogramDATA.mnd 0 // fetch the name of your ascii or hbook histogram (there is a sample data .mnd file attached below)
MN_CMD> index // will list all histograms read into mnfit
MN_CMD>plot 10&0 //to draw histogram 10 (histogram id=10 and secondary ID = 0)
MN_CMD> fun add gaussian // add the Gaussian function as the fitting function, use the sigma type
1st prompt = sigma, 2nd = area, 3rd, mean, 4th = sigme, enter in reasonable starting points for these parameters to start the minimizer
MN_CMD> fit 1 // fit the histogram number 1 with the function loaded into the fitter. You will have to select what type of fit to use. You can use the chi2 fit for this exercise.
MINUIT> minimize // The fit command above puts you into the minimization and error analysis package called MINUIT. MINUIT is a set of functions algorithms and a complete package that's beenaround for ever. It may be of interest to you to spend some time exploring the documentation from CERN)
MINUIT> display // Displays your histogram with the function overlayed and statistical and error analysis results
MINUIT> minos // Recalculates the errors using a more sophisticated approach where the minimizer doesn't assume parabolic errors, for more info see documentation

• To get an output from mn_fit you can either ask mn_fit to output a encapsulated postscript file and paste that into your document or you can just use your laptop's "Print Screen" key (on windows its Alt + Prnt Scrn, on Mac Command + Shift + 4) and then paste onto document. If you opt to use mn_fit's output functions beware that its only usable outputs is postscript or encapsulated postscript. Either will like format will be understood by a printer or a printer driver that reads postscript. You will not be able to display on screen unless you embed a bitmap image in the postscript.
• When you have a fit ready to include in your lab report come and see me to discuss what it means and how to interpret the results. You should include what you learn, with regards to the goodness of fit" the meaning of the chi2 in your reported. This information is reported at the top of the fit graphics displayed by mn_fit. In your write up describe what these are and report on what they mean. Check what happens when you have histograms with lots of events verses those with few. Lots here are of order 10,000 few is of order 100. Are the errors of the number of events per bin computed correctly?

GradingRubric