Program Twopoint.f Description

Upgoing Muon Analysis

Program Twopoint.f Description

John G. Learned
9 November 1996

Being updated 12/96.... stand by please.


This program has evolved over a long period, initially written for analyzing upcoming muon data from the IMB experiment, then incorporating data from other underground experiments which have published lists of upcoming muon times and directions, and later with the addition of a calculation of the twopoint correlation function between arrival directions of muons. The program was written in FORTRAN on DEC VAX computers. It has been ported to a UNIX running ALPHA system in November 1996, and upgraded with the addition of the capability to analyze data from SuperKamiokande.

The data from the program presents a number of histograms (about 60) which are printed to a file and to an hbook file for viewing and manipulation with PAW. We shall describe that output in some detail below, but it falls into two categories: the first type of histogram is to check on data (plotting times, detector coordinates and such), and the second is data analyzed in various ways and displayed in potentially interesting coordinates (e.g. in galactic coordinates, solar coordinates, displaying smoothed sky maps, and showing the twopoint correlation function).

Input data can be selected from a variety of available data sets: the CWI experiment (South Africa), the Kolar Gold Fields experiments (early data), the Baksan Experiment, the IMB experiment, the Kamioka experiment, and finally the Super-Kamiokande experiment. MACRO data will be added when it becomes available. The idea is to be able to add these data sets and make comparisons, particularly in searching for point sources. The program contains a list of favorite suspected sources, and will test for concentrations around those directions.

The program twopoint.f is available in the UHHEPG alpha cluster nfs/users/dumand/jgl/superk/upmu. It needs to link with subroutine package timesub.f, has include files and with the common blocks, and with CERNLIB, all as indicated by the Makefile. The controls may be found in a namelist file called twopoint.par, which can be simply edited to control all choices of data, cuts, processing and output. You can ftp the entire package, or look at data from the files therein. This is an evolving program, and I have been busily improving it lately, so no guarantees come with it... it is tool in development not a production program.

A sample output is available in the /post area, as you must have seen. The graphs are in postscript, The two dimensional ones are unavoidably large, sorry. You can also look at the output in old fashioned lineprinter format in the twopoint.prt file, also located in the /post area.

Program Output

We begin with a narrative description of the output which is written to Unit 6 (FORTRAN default) or to the terminal if you do not specify otherwise. This provides a convenient way to describe the program function and controls.

The program reports times at several points during processing. This is because the computational time can become very long for large data sets (600 events takes several minutes), growing as the square of the number due to the background calculation in twopoint. If you want to make a shorter run, consider turning off the twopoint calculation. The fancy smoothed skyplot is also slow and can be turned off as well.

The controls are reported out first, with a report of which of the eight possible experiment data files were selected. The parameter file contains the file names. The controls are as follows:

AMIN, AMAX: min and max of the asphericity, which only applies (now) to the IMB contained event set.

RJDCUTMIN, RJDCUTMAX: permits selecting events from some particular era, as given by the Julian date minimum and maximum, as floating point numbers (real*8), fractional days.

DISTMIN, DISTMAX: not used now.

DEDXMIN, DEDXMAX: Cuts on acceptable band of energy loss rate in units of MeV/cm.

EMIN, EMAX: energy deposit minimum and maximum, in GeV. Mostly appropriate for contained events, but could be used to eliminate showers (though this is perhaps more easily done with the dE/dX cut).

PLCUT: Requires a minimum path length in the detector of the number of centimeters specified.

CENDISTMIN, CENDISTMAX: are cuts on the minimum and maximum number of centimeters for the impact parameter of the track with respect to the detector center.

FIDISTMIN, FIDISTMAX: allow cuts on the minimum and maximum distance of the centers of tracks from the edge of the fiducial volume. This is a good way to get rid of wall hugging events.

SELECT_MUONS: allows for selecting only muons (for contained events).

VERBOSE: turns on the full output with all diagnostic histograms and data file dump.

IDUM: the random number generator seed. Random numbers are used for Monte Carlo runs of course, but also in normal runs for the smearing of event directions in the twopoint background calculations and in making the fancy skyplot.

SAVE_HISTOS: saves histograms to an hbook file for reading to PAW.

MONTE_CARLO_RUN: set up to generate a bootstrap statistics randomized sample by simply reassigning the times with actual coordinates. However, the detector coordinates are randomized a little bit in a gaussian distribution with magnitude given by SMEAR, in degrees. The time is randomized a little too by 15 minutes. In the code one can also enter a signal at a given location and with specified number of signal events, in order to test the sensitivity of signal finding algorithms.

NTRIALS: sets the number of background runs made with scrambled data sets calculating pair distances for use in the twopoint correlation function normalization. The reason for this technique is simply that the data set is not uniformly distributed on the sky, and thus calculating solid angle will not work well for determining the probability of pairs versus angular distance. One hundred trials seems to be fine, perhaps less would be alright. If unlikely correlations are found then runs with more data sets should be used to see if the results are stable. The scrambling works in the same way as for a Monte Carlo data set. The IMB data is handled differently than other data sets in that the smearing near the horizon is treated carefully so that events do not migrate to smaller zenith angles. Will probably fix it so SK data done this way too.

CUT: is used for a cluster search, assessing the number of clusters within an angular region defined by CUT.

SMEAR: is the angle mentioned above, employed to smear out the angles a little bit. Note that the small angle approximation is used so that angles should be less than about 10 degrees. Probably one degree is about right for SK.

SMOOTH: is the angle used to smear the events in making the background plots for the skyplot. The angular resolution of the detector should be used here, about one degree for SK. At the moment it is one size fits all data sets. There is no harm in using too large an angle here, however, as it only dilutes possible signals by spreading them, over too large and area.


Before writing the plots out, the program dumps the event list. The format is a little different for each detector, and it will do this for each data set. The list gives the event identifications, right ascension, declination, time and local coordinates for each event. The time is reported in UT and in Julian day (note 0.5). The detector coordinates are all regularized to x as East and y as North, with the azimuthal angle measured as N of E. It should be noted that various experiments have used different coordinates systems, with azimuth varying from the above to being E of N or even W of S.


Plots 1-12 deal with the twopoint correlation function. You can cut to the chase and just go to plots 9-11.

Plots 16-31 deal with plotting the muon directions in variour coordinate systems.

Plots 32-45 show distributions in times and angles, mostly for diagnostic purposes (has proved very useful in the past, with stuck bits, time recording problems, and such showing up.)

Plots 46-57 present the nice smoothed sky maps in equatorial coordinates, as well as background calculations, and a particularly nice representation for finding clustering (id=210) in excess (or deficit) of data over background in standard deviation units.

We go through the plots in the order printed (regretably not very logical... as it developed), and describe what is plotted. The plots are listed as HBOOK number/ID.

1/10: Pair separation angle distribution. Shows the number of data pairs versus the angular separation in degrees. This has a hump in the middle, as expected because of the solid angle. This is a differential distribution.

2/11: Pair separation cos(angle) distribution. Same as above but with solid angle binning (0-1). This is generally not flat because the detector does not monitor the whole sky, and thus the number of distant pairs falls off with increasing angular distance. This is why we need a Monte Carlo assessment of the pair distribution to assess background.

3/20: MC pair separation angle distribution. So, as in 1/10, this plots the MC data set. The MC data is generated from the real data, bootstrap statistics style, by shuffling the times in the data set, and stirring the angles a little too.

4/21: MC pair cos(angle) distribution. As indicated above, this shows the effects of anisotropic sky observation.

5/30: Pair separation angle distribution, data/background. Now we have divided the real data set distribution by the ensemble of shuffled data sets, bin by bin.

6/31: Pair separation cos(angle), data/background. As above.

7/40: 2Pt correlation Function versus angle. This finally is the twopoint correlation function, in differential form.

8/41: 2Pt correlation of cos(angle). Ditto.

9/50: Cumulative 2 point correlation as a function of angle. This is really more like what we want to see, since we are searching for correlations due to some nearly point source, and the correlation function will maximize within some angle typical of the angular spread of the source as convoluted through neutrino and muon scattering and track reconstruction.

10/51: Cum 2pt CF of cos(angle). Ditto.

11/60: Standard deviation of the 2pt CF versus angle. This is the place to look for significance of any effects of clustering.

12/61: SD of 2pt CF Vs cos(angle). Ditto.

13/120: Number distribution within cut. Not currently filled.

14/121: MC Number distribution within angle cut. Also turned off for now, as this takes much time.

15/122: Obsvd/Exptd within angle cut. Ditto.

16/69: Right Ascension - Declination distribution of close pairs of muons. Ditto.

17/70: RA-DEC distribution of muons. This gives the equatorial coordinate plot of the events.

18/71: RA-sin(DEC) distribution of muons. Ditto, except for using sine of declination.

19/72: RA distribution of events. Right Ascension projection of event directions.

20/73: Declination distribution of events. Declination projection, ditto.

21/74: Sin(Declination) distribution of events. Ditto.

22/76: Galactic Coordinate map of events directions.

23/77: Galactic Longitude projection.

24/78: Galactic Latitude projection.

25/79: Ecliptic Coordinate map of event directions.

26/80: Ecliptic Longitudinal projection.

27/81: Ecliptic Latitudinal projection.

28/82: Local Coordinate map. In azimuth (N of E) versus zenith angle of origin.

29/83: Local Azimuth angle projection.

30/84: Local Zenith angle projection.

31/840: Local cos(zenith angle) Distribution.

32/109: Anisotropy distribution of events. Applicable to contained data only.

33/85: Energy deposited in detector, GeV.

34/110: Path length in detector, in cm.

35/111: dE/dX, in MeV/cm.

36/86: Julian date of events, -2445000.

37/87: Local time of events.

38/101: Local hour distribution of events. Mainly for checking input data.

39/102: Local minute distribution of events. Ditto.

40/103: Local seconds distribution of events. Ditto.

41/104: Decimal GMT distribution of events.

42/105: Decimal local sidereal time distribution of events.

43/106: Day distribution of events.

44/107: Month distribution of events.

45/108: Year distribution of events.

46/200: Weighted sky plot. This gives the events in RA and DEC, with event locations smeared out on the sky according to the parameter SMOOTH. The normalization is such that the peak of an event has value 1.0, and thus the amplitude at any point represents something like the number of events that may have originated in that direction.

47/210: Weighted sky plot in standard deviation units. This plots the events after they have been compared to the MC sky map as discussed below. The plot is shifted up by 2sd, to allow for negative values, and scaled up by 5.0. Hence 10 units corresponds to zero standard deviations. Three SD positive would be 25.

48/220: Distribution of values in the SD plot. This give the histogram of the values from all bins, so there is much overcounting (depending upon the value of SMOOTH compared to the bin size).

49/205: MC declination projection. This shows the declination distribution used to generate the background events. The MC for this purpose, assumes that the data is uniformly distributed in RA.

50/206: Smoothed MC declination distribution. The data used in the MC is smoothed to help with low statistics.

51/207: MC zenith angle distribution.

52/ 208: MC azimuth angle distribution.

53/209: MC Local coordinate distribution. This is mainly for checking the program.

54/201: The events plotted in RA and DEC. Larger scale than in id =70, and same as 210 and such for making overlays.

55/202: The MC events in RA and DEC. Ditto.

56/202: Galactic plane for overlay.

57/203: Some favorite sources, also for overlay.