Cascade reconstruction version 2

The cascades are one of the two detection channels of neutrinos in the neutrino telescopes. Recently, the most energetic events in the IceCube are cascades. The goal of the cascade reconstruction is to identify their distinctive pattern in the measured data and reconstruct from the measured hits the position and time of the cascade and also its energy and direction. The first two parameters are obtained from the Chi^2 fit. The latter ones are calculated with the log likelihood technique.

1. Instructions

New BARS program for cascade reconstruction in $BARSSYS/programs/reconstruct-cascade/. Right now it works only with real data. No MC simulation of cascades in BARS format has been provided yet. Moreover, only the first step of the cascades reconstruction (position and time) is implemented for the time being. For precise results, the time calibration (offsets and timecal_dzh) and real dynamic geometry (geometry.export-linear.root) is applied based on the parameters of the program.

2. Inputs

  • XeventID.joint.events.root - Joint BExtractedImpulseTel files
  • geometry.export-linear.root - dynamic geometry files
  • $BARSSYS/build/config/year/calib/clusterY/offset/offsets - offset calibration files
  • $BARSSYS/build/config/year/calib/clusterY/time/timecalib_dzh - time calibration files

3. Outputs

In the same folder where XeventID.joint.events.root is located are produced three files:

  • recCasc_hist.root - Summary histograms for better understanding of the processed run
  • recCasc_vis.root - Visualizations of events that passed all the used filters/cuts
  • recCasc_nTuple.root - Important parameters of the reconstructed cascades

recCasc_hist.root

The list of basic summary histograms for the whole run as well as for sub-sets of filtered cascades on different levels of filtering.

  • h_nHits - number of all hits per event
  • h_chargePerHit - charge distribution of all measured pulses for all OMs
  • h_nHitsAfterQ - number of hits that passed Q filter per event
  • h_chi2AfterQ - the resulting chi2 after the first fit
  • h__nHitsAfterT - number of hits that passed T filter per event
  • h_chi2AfterT - the resulting chi2 after the second fit
  • h_timeRes - time residuals of all pulses from events that passed charge chi2 filter
  • h_nHitsPerOM - number of hits on the given OM. This histogram is used to decide if the OM is working or "dead".
  • h_XYposition - X and Y position of cascades that passed Z filter/cut

The list of histograms will be definitely extended in the future based on the needs of users. Do you have suggestion for new histogram, please, let us now.

recCasc_vis.root

For every reconstructed cascade that passed all the filters/cuts the TCanvas with 9 TPads is created. The first eight TPads show hits (black points) and expected times of detection (region defined by two green lines) on individual lines in Time-Z space. The red point represents the fitted time and position of the cascade. The ninth pad shows the distance vs. charge of all the pulses in the expected time region (green lines). This last pad helps to better distinguish between near to horizontal atmospheric muon bundles from cascades. In the case of cascade, the charge decreases with distance almost exponentially and no high charge should be present at larger distances.

recCasc_nTuple.root

Summary of the parameters stored in the NTuples:

  • runID - cascade run ID
  • EventID - cascade event ID
  • NPulses - number of pulses in the event
  • NPulsesT - number of pulses produced by cascade
  • QRatio - ratio of the charge of the pulses within and without +-50 ns time window around expected time
  • CloseHits - number of close OMs above the cascade position that were hit by cascade (maximum = 15)
  • X - reconstructed x position
  • Y - reconstructed y position
  • Z - reconstructed z position
  • Time - time of cascade production (calibrated one)

4. Program flow

  1. Existence of all input files is checked and user is informed in the case of missing ones

  2. joint.events.root

  3. offsets
  4. timecalib_dzh

  5. geometry.export-linear.root

  6. The loop over all the events in the run is started

  7. Batch of filters/cuts inter-spaced by 2 fits is applied to every event in run in the following order:

  8. N Filter - only events with number of pulses bigger than NCut (default value = 70) pass

  9. Q Filter - only events with more than QCutHits pulses (default value = 6) with charge larger than QCut (default value = 2000 FADC channels) pass
  10. Initial cascade parameter estimation - initial time and position of the cascade is estimated as the time of the pulse with the largest charge and the position of OM that detected this pulse
  11. First fit - the first chi2 fit of only pulses with charge larger than QCut (default value = 2000 FADC channels) is performed
  12. Chi2 Filter - Chi2 of the fit has to be smaller than QCutChi2 (default value = 50)
  13. T Filter - only events with more than TCutHits pulses (default value = 20) in the +- TCutTimeWindowNs window (default value = 50 ns) around the expected time pass
  14. Second fit - the second chi2 fit of only pulses within the +- TCutTimeWindowNs window (default value = 50 ns) around the expected time is performed
  15. Chi2 Filter - Chi2 of the fit has to be smaller than TCutChi2 (default value = 50)
  16. Z Filter - only events with cascade reconstructed Z position less than ZCut (default value = 600) pass
  17. T Delay Filter - time delay between cascade time and time of the first hit is calculated and this value has to be smaller than TDelayCut (default value = 400 ns)
  18. Q Ratio Filter - sum of the charge of all the pulses within the +- TCutTimeWindowNs window around the expected time is calculated and it is divided by the overall charge of all the pulses in the event and this ratio has to be larger than QRatioCut in percents (default value = 85 percent)
  19. Branch Filter - More than half of the pulses within the +- TCutTimeWindowNs window around the expected time has to be above the cascade Z position
  20. Close hits Filter - Pulses have to be found on more than 10 out of 15 WORKING closest OMs that are above the reconstructed cascade Z position

  21. Event visualization - the visualizations of the events that passed by all the above mentioned filters are rendered

  22. NTuples, histograms and visualizations are saved

Obligatory format of time calibration files

Time calibration as well as offsets are automatically read from $BARSSYS/build/config/ folders. Both calibration files have to have fixed format to enable automated reading. 

1. # first line with comment
2. empty line
3. # String 1
4. 12 values separated by spaces or tabs for the first section
5. 12 values separated by spaces or tabs for the second section
6. 12 values separated by spaces or tabs for the third section
7. empty line
8. # String 2
9. 12 values separated by spaces or tabs for the first section
10. ...

Obligatory file structure of the input files

After BARS installation described here, set in your .bashrc following environment variables according to your system. 

export BARSSYS=/home/fajtak/work/bars2/build
export PATH=$BARSSYS/bin:/$HOME/.local/bin/:$PATH
export LD_LIBRARY_PATH=$BARSSYS/lib:$LD_LIBRARY_PATH

export CEPH_MNT="/media/fajtak/Data/BaikalData"
export BARS_DATA_PROCESSED=$CEPH_MNT
export BARS_DATA_JOINT=$CEPH_MNT
export BARS_DATA_RAW=$CEPH_MNT/raw
export BARS_DATA_GEOMETRY=$CEPH_MNT/geometry

All the joint events should be saved in the CEPH_MNT folder. All geometry files should be saved in the BARS_DATA_GEOMETRY folder. The structure of the folders in the BARS_DATA_GEOMETRY has to follow the structure of geometry folder at Baikalsnr. In the CEPH_MNT create folders exp19_barsv051, exp18_barsv051, exp17_barsv051 and exp16_barsv051 and in every folder create folders cluster0, cluster1, ..., clusterN according to the year. In these folders there are folders of runs (for example 0025) and it this folder you can put h0025.joint.events.root. With this folder structure the whole processing will be automatic and you don't need to set any input and output paths.

5. Usage

Instalation

Go to the bars program folder

cd bars/programs/

Download the reconstruct-cascade program from git

git clone https://github.com/fajtak/reconstruct-cascade.git

Go to the build folder

cd reconstruct-cascade/build

Create a MakeFile and compile

cmake ..
make

Install it

make install

and start it

./reconstruct-cascade -s 2018 -c 2 -r 25
./reconstruct-cascade -s 2018 -c 2 -r 25 -x ../src/config.rc -n 25000

The first example above shows the minimal working example. You always have to specify season (-s), cluster (-c) numbered from zero, and run (-r). Other optional parameters are number of processed events (-n) and your own configuration file, where you can change the values of cuts without recompiling the whole program.

Configuration parameters

The meaning of the parameters is explained in Program Flow above. If you want to use your own values, you can change the values in src/config.rc file and then start program with parameter -x ../src/config.rc.

Default values are following:

NCut 70
QCut 2000
QCutHits 6
QCutChi2 50
TCutTimeWindowNs 50
TCutHits 20
TCutChi2 50
ZCut 600
TDelayCut 400
QRatioCut 85

Helpful bash scripts

The program is accompanied by the set of bash and root scripts that are located in the build folder.

  • changeDataSource
  • downloadJoint
  • downloadJointAllRuns
  • processAllRuns

  • deleteAllResults

  • studyReconstructedCascades.C

changeDataSource

How to use:

source ./changeDataSource 0
source ./changeDataSource 1

This script can be used to change the data storage without rewriting the .bashrc. Just change the CEPH_MNT="" in the script to your second data storage (in my case external HDD drive). Then when you start a new terminal where CEPH_MNT is set to your usual local datastorage you can call the script in the above mentioned way and the data storage will be changed. You can define two different options for parameter 0 and 1 or add another ones.

downloadJoint

How to use:

./downloadJoint 18 2 25

This script downloads the joint.events.root from given year, cluster and run to your data storage to appropriate folder. To specify year use just last two digits. Clusters are numbered from zero. To be able to do that you have to mount cephfs to your computer with this command:

sudo sshfs baikal@baikalsnr.jinr.ru:/mnt/cephfs /mnt/cephfs -o reconnect -o allow_other

downloadJointAllRuns

How to use:

./downloadJointAllRuns 18 2

This script does the same what the previously described one but it downloads all runs from the given year and cluster. To specify year use just last two digits. Clusters are numbered from zero.

processAllRuns

How to use:

./processAllRuns 18 2

This script calls ./reconstruct-cascade for all runs from the given year and cluster. To specify year use just last two digits. Clusters are numbered from zero.

deleteAllResults

How to use:

./deleteAllResults 18 2

This script deletes all output files (recCasc_*.root) from all runs from the given year and cluster. To specify year use just last two digits. Clusters are numbered from zero.

studyReconstructedCascades.C

How to use:

root -l
.x studyReconstructedCascadesC+(18,2)

This script takes results from all recCasc_nTuple.root from given year and cluster and lets you to do your own analysis of the well reconstructed cascades. In the script there are examples how to use different variables and create histograms and graphs.