Muon reconstruction
1. Instructions
Script macros/rec/stage1_mc_NEW.C is most up-to-date script.
Usage
Install BARS as described in this twiki. Script is launched from ROOT with the following syntax:
.x stage1_mc_NEW.C+("<WOUT file with MC events>","<output ROOT file with reconstructed BARS objects>", 0, <noise suppression method id>, <threshold of hit selection>)
Noise suppression method is one of 3 methods described in chapter 2:
- Causality and muon criterion
- Hotspot selection and string group building
- Causality and residual cleaning
At the output 3 files are produced:
- file with the event mask monitoring histograms (default: mask_monitor.root)
- file with reconstruction monitoring histograms (default: reco_monitor.root)
- file with output BARS object and other data (name is given in the launch command)
The output BARS fomat is described in chapter 4
2. Noise suppression
Large fraction of impulses falling into the 5mks event frame are due to noise. With main contribution at low amplitudes being PMT dark current and lake noise (amplitudes of order of 1 p.e.). At higher amplitudes (which?) contributes noise due to discharge in OM electronics. Sample of signal impulses in the event should be selected with as high as possible purity while keeping the detector sensitivity at acceptable level. There are thee methods of signal impulse selection implemented currently
- Causality and muon criterion.
- Hotspot and string group building. In this method pair of adjacent OM's with impulses with amplitudes above A1 and A2 which fall into causality window are found at each string. Around such pair a group of impulses which follow the cerenkov signal propagation pattern along the same string are selected. At the next stage groups of impulses at different strings are merged into global group. The group which shows the best quality/nhits ratio is seected for further analysis. Configuration block contains the following parameters (corresponds to rec/stage1_mc_NEW.C):
hotspotHitGroups.SetOutputMaskName(selectedHitsMaskName);
this parameter defines the event mask name which will be used by next stages of processing,
hotspotHitGroups.SetHotSpotThreshold1(thresh);
hotspotHitGroups.SetHotSpotThreshold2(thresh);
here thresholds for hotspot selection are set,
hotspotHitGroups.SetSafetyWindow(30);
here the time (ns) margin for causality criteria for impulses at adjacent OM's is set
hotspotHitGroups.SetSafetyWindowHotSpot(30);
time (ns) margin for assembling the group at the same string
hotspotHitGroups.SetLinkWindow(30);
time (ns) margin for linking of groupd at different strings
hotspotHitGroups.AssembleGroups(kFALSE);
if this parameter is set to kFALSE groups string groups are not constructed but instead just two impulses from hotspot are linked with another string
hotspotHitGroups.IsData(kFALSE);
this should be set to kTRUE if real data is processed. Described method requires signal hits at adjacent OM's. About 90% of tracks fail this topological criterium, so that method results in low sensitivity to the neutrino signal. Although it works well for the atmospheric muon flux reconstruction. 3. Causality and residual cleaning. In this method all loop over all impulses is done causality-connected group is constructed starting from each impulse and passed to further stges of selection afterwards the group with the best quality/nhits ratio is selected. After the causality group is constructed the estimation of trajectory is performed as follows. The direction is obtained using sum of vectors defined by all possible pairwise combinations of impulses at different strings ordered in time. The reference point is defined by position and time of impulse with largest amplitude. Using preliminary trajectory estimation cleaning of residuals is performed in two stages as follows. At the first stage cleaning of time and distance residuals is done. It is done in 9 iterations with gradual tightening of cuts on residual time and distance (rho). Values of cuts change from 1000 ns to 300ns for time and from 300m to 60m for rho. ...This should be added to config... Next step is cleaning of residuals in phit variable. It is done in 5 iterations with tightest value of phit being 10e-4. This method provides purity of order of 90% for threshold of 1.5 p.e. and 10 times better efficiency than method 2. Configuration parameters are as follows:
hitGroups.SetOutputMaskName(selectedHitsMaskName);
event mask name
hitGroups.SetHotSpotThreshold1(thresh);
hitGroups.SetHotSpotThreshold2(thresh);
thresholds for initial impulse in causality group and threshold for the rest of impulses
hitGroups.SetLinkWindow(30);
time safety margin for causality criterion
hitGroups.ExcludeTimeRhoOutliers(kTRUE);
do the time and rho residual cleanin
hitGroups.ExcludePhitOutliers(kTRUE);
do the phit residual cleaning
hitGroups.IsData(kFALSE);
true if processing of real data is performed.
3. Muon trajectory reconstruction
Method for the reconstruction
Muon direction is reconstructed in two distinct stages: fast initial estimation of trajectory and minimisation of quality function to find optimal trajectory. INitial estimation is done as described in chapter 2: for the direction estimation sum of vectors is used which is defined by all possible pairwise combinations of impulses at different strings ordered in time, the reference point is defined by position and time of impulse with largest amplitude. Next stage is minimisation of the quality function. There are 9 functions implemented currently, although performance results are so far obtained just for function 1:
1
Q=\Sigma_{1}^{Nhits} \frac {(t_{i}^{exp}-t_{i}^{theor})^2}{\sigma_{i}} + w*A*D
quality function inspired by ANTARES paper ... A and D are proportional to distance and time, w - reltive weight of second member. Presently set to 0.2.
2
Q=\Sigma_{1}^{Nhits} \frac{(t_{i}^{exp}-t_{i}^{theor})^2}{\sigma_{i}}
regular time chi2
3
Q=\Sigma_{1}^{Nhits} \frac{(t_{i}^{exp}-t_{i}^{theor})^2}{\sigma_{i}} + w*D
time chi2 and distance member
4
likelyhood function under development
5
M-estimator used by ANTARES (https://arxiv.org/abs/1708.03649)
Q=\Sigma_{1}^{Nhits} q_{i}\sqrt{1+(f(t_{i}^{exp}-t_{i}^{theor}))^2/2}
where q_{i} is amplitude in i-th channel, adjustable f - factor, set to 1 presently.
6
modified M-estimator
Q=\Sigma_{1}^{Nhits} q_{i}\frac{1}{1+1/(f(t_{i}^{exp}-t_{i}^{theor}))^2}
7
maximisation of phit product
Q=-\Sigma_{1}^{Nhits}\log_{10}p_{hit}(r)
8
combination of chi2 and phit
Q=\Sigma_{1}^{Nhits}(w\frac{(t_{i}^{exp}-t_{i}^{theor})^2}{\sigma_{i}}-\log_{10}p_{hit}(r))
9
robust distance function with parameter mu=40
Q=\Sigma_{1}^{Nhits}f_{i}
f_{i} = \rho_{i}^2, \,\,\, \rho_{i} < \mu
f_{i} = \mu (2 \rho_{i} - \mu), \,\,\, \rho_{i} \geq \mu
Program flow
Class BMuonTrajectory hold trajectory parameters and contains all methods for the reconstruction: initial trajectory estimation, minimisation, residual cleaning, etc..
Class BMuonTrajectoryProducer is child class of MTask class. Int its BMuonTrajectoryProducer::Process() method it creates the BMuonTrajectory object using the event mask obtained from earlier stages of processing, calls methods for trajectory reconstruction and fills the BRecParameters container.
4. Output of the reconstruction
RecParameters contains the following recontsructed quantities:
| type | name | meaning |
|---|---|---|
| Double_t | fFuncValue; | value of minimisation function |
| Int_t | fQual; | NOT FILLED |
| Int_t | fNumPar; | NOT FILLED |
| Int_t | *fIsFixed; | [fNumPar] NOT FILLED |
| Double_t | *fInitial; | [fNumPar] NOT FILLED |
| Double_t | *fMinLimit; | [fNumPar] NOT FILLED |
| Double_t | *fMaxLimit; | [fNumPar] NOT FILLED |
| Double_t | *fStep; | [fNumPar] NOT FILLED |
| Double_t | *fErrors; | [fNumPar] errors of trajectory parameter estimation |
| Int_t | fNhit; | number of impulses used for the reconstruction |
| Double_t | *fTres; | [fNhit] redsidual times for channels: t_exp-t_theor |
| Int_t | *fNchGeom; | //[fNhit] ID's of channels |
| Int_t | fNhitMC; | // NOT FILLED |
| Double_t | *fTresMC; | //[fNhitMC] NOT FILLED |
| Int_t | *fNchGeomMC; | //[fNhitMC] NOT FILLED |
| Double_t | fThetaRec; | // reconstructed polar angle in degrees in (0-180). Zero angle corresponds to the vertical upward direction of particle velocity. |
| Double_t | fThetaMC; | // MC zenith angle in degrees in (0-180), the same definition of angle. |
| Double_t | fPhiRec; | // reconstructed azimuth angle in degrees in (0-180). Angle is calculated wrt. X-axis direction (...) |
| Double_t | fPhiMC; | // MC azimuth angle in degrees in (0-180), same definition (???) |
| Double_t | fX0Rec; | // reconstructed x0 in meters. X0 is the coordinate of trajectory point in the plane Z=0. |
| Double_t | fX0MC; | // NOT FILLED |
| Double_t | fY0Rec; | // reconstructed y0 in meters. Y0 is the coordinate of trajectory point in the plane Z=0. |
| Double_t | fY0MC; | // NOT FILLED |
| Double_t | fTimeRec; | // Rec time in ns |
| Double_t | fTimeMC; | // MC time in ns |
| double | fP; | // NOT FILLED |
| double | fPMC; | // NOT FILLED |
| Double_t | fFuncValueMC; | //NOT FILLED |
| double | fE; | //NOT FILLED |
| double | fZDist; | //ZDist - largest distance between OM projections on the track |
| double | fZDistMC; | //NOT FILLED |
| int | *fImpulseNumber; | //[fNhit] numbers of impulses used in the reconstruction |
| double* | fRho; | //[fNhit] distances from trajectory to OMs used in the reconstruction |
| int | fNcalls; | // number of function calls during minimization |
| TVector3 | fInitialDirection; | //vector of initial direction estimation |
| TVector3 | fInitialPoint; | //vector of reference point estimation |
| Double_t | fInitialT; | //time at reference point |
| Double_t | fInitQuality; | //value of minimization function for the initial approximation |
| Double_t | fPHit; | //product of Phit - probability to give active hits |
| Double_t | fPnonHit; | //probability to hav given configuration of empty OM'sat the distance r<30m from track |
| Double_t | fPHitSum; | //sum of Phit |
| Double_t | fTResMin; | //mimimal residual time among all impulses |
| Double_t | fTResMax; | //maximal residual time |
| Double_t | fRhoMax; | //maximal rho |
| Double_t | fRhoMin; | //minimal rho |
| Int_t | fNnonHit; | //number of OM in pnon-hit caculation (r<30m from the track) |
5. Reconstruction monitoring
Modules for the monitoring of hit selection and track reconstruction are run from the same program. At the output two additional files are produced which contain the monitoring results.
6. Action items (20/06/18)
- Reconstruction and monitoring automation
- Certified data processing chain
- Optimisation of noise suppression: adjust intrinsic algorithm parameters