RICH PID

Page maintainer: Rebecca

This page is considered done. It been reviewed by Yoshi. There may be missing elements, but they are all flagged and the text has no errors.

Intro

The Ring Imaging CHerenkov (RICH) detector is used to identify pions, kaons and protons at HERMES. It is described in Bootcamp. While the track identity is in the uDST tables, to be more precise you should do "RICH unfolding", which gives each track a probabilistic pion, kaon, and protons weight. See below for more details.

In 1995-1997 the RICH did not yet exist, and the Cherenkov Detector was used.

Essentials

A beautiful (but slightly outdated) description of the RICH detector, algorithms, and unfolding methods can be found on the RICH web page under Using the RICH Detector for Physics Analysis . HOWEVER, please note that the e-, h-, and c-tune Pmatrices are no longer used to produce systematic errors AND all references to the Q matrix should be ignored, only the inverted Pmatrix, $P^{-1}$ matrix is used (although in some places references to the Q matrix in reality meaning the Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle P^{-1}} matrix. Just be careful.)

Methods

The RICH information is processed with three different algorithms, IRT, DRT, and (for multiple tracks in a detector half) EVT. Additionally, the RICH PID Scheduler tries to choose the BEST method between IRT and DRT. This is described in more detail in Bootcamp.

Currently, the recommended method is EVT for multiple tracks in a detector half and DRT for single tracks.

The uDST tables

The g1Track table links to the smRICH table via four different links, IRT, DRT, EVT, and BEST. To use the recommended method, EVT+DRT, one can use:

if (COUTAB(smRICH) > 0) { /* Guard against case where table is empty */
  NATREL(g1Track, g1Track.EVT, smRICH, ok);
  if (ok && smRICH.iMethod==3) {
    /* If this is really EVT (method=3) then take this information   */
    itype=smRICH.iType;
  } else {
    /* if there is no EVT information (single track OR EVT failed),  */
    /*   get DRT information                                         */
    NATREL(g1Track, g1Track.DRT, smRICH, ok);
    if (ok) {
      itype = smRICH.iType;
    } else {
      itype = -999;
    }
  }
} else {
  itype = -999;
}

For some productions (and all future productions) the EVT link is set to DRT for a single track in the detector half. See the Status of EVT Productions section below for details. For these productions, only the first line and the check on the "ok" variable is needed.

The smRICH table then contains the relevant infomation for the chosen track and method. The most important entries are

iType - the identified particle type with the code
- 0 - rQp of 0
- 3 - pion
- 4 - kaon
- 5 - proton
rQp - the log of the ratio of the two highest probabilities (slightly more complicated for EVT). If this is zero it means two particle type hypotheses are equally probable and it is not possible to decide on the particle type. Typically user make a cut rQp>0

RICH Unfolding (Pmatrices)

To further improve the RICH identification, an unfolding is done. Pmatrices give the probability that a track that is truly of type A is identified as type B, where A and B range over pions, kaons, and protons. They are constructed using Monte Carlo where both the true and identified types are know, and are a function of track multiplicity and momentum. By inverting the Pmatrix to get the Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle P^{-1}} matrix one can then unfold individual tracks by giving them a pion "weight", a kaon "weight" and a proton "weight".

The current Pmatrices are available at

/group01/richgrp/Pmatrix_v4.0

PLEASE be sure you are reading in the Pmatrices/ inverted Pmatrices correctly by recreating the example plot found below in the Example Codes section.

Choosing Your Method

There are a few decision you need to make before you begin your unfolding.

Geometry - Based on the data set you are unfolding, choose between
- 1999 geometry, for 99-05 data
- 2006 geometry, for 06-07 data

Method
- EVT - recommended for more analyzers
- IRT- provided only for historical reasons!!!

Using this information you can navigate to the correct directory from /group01/richgrp/Pmatrix_v4.0

Track counting

When determining which Pmatrix you should use, the number of tracks in the detector half should be counted before any cuts at all are made. All long tracks (check g1Track.iSelect to be sure you do not have magnet tracks) should be counted.

Valid Momentum Range

The Pmatrices are provided for 1-16GeV, but in most cases analyzers should only use the RICH in the range 2GeV - 15GeV. There are a few exceptions:

If you are only interested in pions, it is possible to go down to 1Gev.
The systematic errors for kaons and especially for protons are large at low momentum. For analyses that are RICH systematics limited (and not statistically limited) a minimum momentum cut of 4GeV should be considered for protons and possibly for kaons.

If you are using a subset of the available hadron types (pions, kaons, protons) then you should be careful to do the RICH unfolding correctly. In such a case you can make, for example, a "pion only" Pmatrix, where the columns/rows are Pion and Not-a-Pion. These can be constructed from the "pmatyields" files that are provided by adding together the appropriate rows and columns and then normalizing by the total in each column (including the unidentified category).

How to create "pion only" inverted Pmatrix

/group01/richgrp/Pmatrix_v4.0/Geom99/EVTall/center/pmatyields_EVT_all.asc

(or whatever file is appropriate for your geometry and which p-matrix version you want if you're using this for a systematic).

Look at the first four non-comment lines in this file you're see the raw yields for the 1-2GeV bin for a single track event. The columns give the number of true pions / kaon / protons and the rows are the number of tracks identified as pions (1), kaons (2), protons (3), unidentified (4).

#tracks Pbin Pmin Pmax h(id)=[1...4] pi(true) k(true) p(true)
1  1     1.0     2.0    1     156363.000        522.000       3782.000
1  1     1.0     2.0    2          3.000         19.000         16.000
1  1     1.0     2.0    3        100.000        188.000       1195.000
1  1     1.0     2.0    4       5223.000       8531.000      68170.000

The pmatrix version of this file is pmat_EVT_all.asc and the corresponding lines are

#tracks Pbin Pmin Pmax h(id)=[1...4] pi(true) k(true) p(true)
1  1     1.0     2.0    1       0.967060       0.056371       0.051693
1  1     1.0     2.0    2       0.000019       0.002052       0.000219
1  1     1.0     2.0    3       0.000618       0.020302       0.016333
1  1     1.0     2.0    4       0.032303       0.921274       0.931755

All that's been done is each column is normalized by the column sum (just over these four rows, not the whole file). For example, 522 / (522+19+188+8531) = 0.056371.

So, what you need to do is

1 - convert pmatyields_EVT_all.asc into a pion / not-pion format, so our small example it would have one less column and two less rows. The first few lines would look like this:

# #tracks Pbin Pmin Pmax h(id)=[1...4] pi(true) not-pi(true)
1  1     1.0     2.0    1     156363.000        4304.000
1  1     1.0     2.0    2     5326.000         78119.000

2 - calculate a p-matrix from your new pion / non-pion file, in the same way as pmatyields_EVT_all.asc was converted into pmat_EVT_all.asc. The first few lines would be:

# #tracks Pbin Pmin Pmax h(id)=[1...4] pi(true) not-pi(true)
1  1     1.0     2.0    1     0.9670        0.0522
1  1     1.0     2.0    2     0.0329        0.9478

3 - invert this matrix, to get something of the form of invpmat_EVT_all.asc.

Pmatrix file formats

Pmatrix

The main Pmatrix files are the pmat_METHOD_all.asc

Where METHOD = IRT or EVT

Within each file you have something like this

# P matrix
# Version: 4.0
# RICH PID method: EVT
# Charge sep: all
# MC source: disNG 99 with own BKG
# 
# h(id): 1:pi 2:k 3:p 4:unidentified
# Format:
# #tracks Pbin Pmin Pmax h(id)=[1...4] pi(true) k(true) p(true)
1  1     1.0     2.0    1       0.967060       0.056371       0.051693
1  1     1.0     2.0    2       0.000019       0.002052       0.000219
1  1     1.0     2.0    3       0.000618       0.020302       0.016333
1  1     1.0     2.0    4       0.032303       0.921274       0.931755
1  2     2.0     3.0    1       0.962514       0.022842       0.039187
1  2     2.0     3.0    2       0.004946       0.785202       0.042091
1  2     2.0     3.0    3       0.032535       0.191957       0.918722
1  2     2.0     3.0    4       0.000004       0.000000       0.000000

The header specifies the P matrix type (Version, PID method, charge separation, type of MC production the P matrices were produced from).

The columns are:

Track multiplicity
Momentum bin
Minimum momentum
Maximum momentum
IDENTIFIED as H type (H is: 1=pion, 2=kaon, 3=proton, 4=unidentified particle)
Fraction of TRUE pions identified as H
Fraction of TRUE kaons identified as H
Fraction of TRUE protons identified as H

inverted Pmatrix

The inverted matrices also exist in the invpmat_METHOD_CHARGESEP.asc files. These have the format:

# Inverted P matrices
# Version: 4.0
# RICH PID method: EVT
# Charge sep: all
# MC source: disNG 99 with own BKG
# 
# Format:
# #tracks Pbin Pmin Pmax h(true)=[1...3] pi(id) k(id) p(id)
1    1     1.0     2.0    1    1.0361067    4.5894997   -3.3407559
1    1     1.0     2.0    2   -0.0062370  561.8389659   -7.5136427
1    1     1.0     2.0    3   -0.0314511 -698.5422761   70.6916403
1    2     2.0     3.0    1    1.0405166   -0.0196392   -0.0434822
1    2     2.0     3.0    2   -0.0046308    1.2880708   -0.0588151
1    2     2.0     3.0    3   -0.0358806   -0.2684329    1.1022972

The header specifies again the inverted P matrix type (Version, PID method, charge separation, type of MC production the P matrices were produced from).

The columns are:

Track multiplicity
Momentum bin
Minimum momentum
Maximum momentum
TRUE type H (H is: 1=pion, 2=kaon, 3=proton)
H weight for IDed pions
H weight for IDed kaons
H weight for IDed protons

Pmatrix hadron yields

In addition, the hadron yields which were used to produce the P matrices are available in the pmatyields_METHOD_all.asc files. The output format is identical to the pmat_METHOD_all.asc files, but instead of probabilities they show the total number of hadrons of type X identified as Y in the respective P matix bin.

Systematic Error

In order to provide a systematic error on the RICH unfolding several different Pmatrices are provided for each of the combinations above (geometry, and method). Previously this is where the c-, h-, and e- tunes were used. Analyzers should use all the provided Pmatrices and use the greatest difference as a systematic error.

Generators - Since MC does not perfectly reproduce our data we use both Pythia and disNG to produce Pmatrices, each using a background file from the production itself.

Background files - The DRT and EVT algorithms include a value for the probabilistic number of "background" "hits" in a given PMT. This values is contained in the "background file" which can be extracted from an set of HRC files. Because the background files extracted from different data years and different MC production vary, Pmatrices were produced with disNG using a variety of different backgrounds:
- disNG - the background from the disNG production itself, this is our best guess at the "correct" Pmatrix
- pythia - the background from a pythia production of the same geometry
- data - the background extracted from data of a year with similar geometry (2000 data for 99geometry, 2006 data for 06geometry)

There is a link in each Pmatrix directory to "center", which is the Pmatrix that should be used as the central value, the disNG production with its own (disNG) background file.

NOTE: If your analysis depends an a variable that is integrated over in the Pmatrices (ie. Fehler beim Parsen (MathML mit SVG- oder PNG-Rückgriff (empfohlen für moderne Browser und Barrierefreiheitswerkzeuge): Ungültige Antwort („Math extension cannot connect to Restbase.“) von Server „https://wikimedia.org/api/rest_v1/“:): {\displaystyle \phi} ) you should consider doing your own PEPSI challenge (aka MC self-test) to determine the systematic error due to the RICH.

More Info

The old RICH webpage has a wealth of information, some of which is outdated.

Pmatrix versions

The current Pmatrix version is v4.0

v4.0
- First Pmatrix set containing the new EVT PID method (which is recommended)
- The different tunes (e/h/center-tune) are discarded
- IRT method is only provided for "historic" reasons

Cuts on the Pmatrices

The pmatrices were produced with the following cuts:

Event level cuts
- require beam lepton
Track level cuts
- RICH EVT link exists (long track)
- fiducial volume cuts (see the Standard Cuts page)
  - z vertex cut
  - calo position
  - field clamp acceptance cuts

Status of EVT Productions

last updated by Rlamb 22:58, 21 July 2010 (UTC)

Papers and Notes

IRT -- The IRT algorithm is described in Ralf Kaiser's internal note

DRT -- The DRT algorithm is described in E Cisbani's internal note

EVT -- Lots and lots and lots of details on the work done to create EVT, including how to use EVT, systematic studies done, and details on the various Pmatrices can be found in Rebecca's internal note

BEST -- The RICH PID Scheduler is described in A. Mass' internal note

Pmatrices
- Work done to tune the MC (and the source of the e-tune (electron tune) mirror roughness parameter) is documented in Brecht Hommez's RICH note
- Old study on Pmatrices with decay particles (source of the h-tune (hadron tune) mirror roughness parameter) is documented in the Hasagawa and Miyachi internal note
Many many many other papers, if they are not in the HERMES documents database, can be found on the RICH Internal notes page

Code Repository

Test your RICH unfolding by producing the Example plot! Here you will find details on how to produce the plot and the final plot you should reproduce.

NEEDEXPERT Can someone please xcheck my plot!!! (I already xchecked it with my analysis code vs. Hanna++) Rlamb 12:07, 22 May 2008 (CEST)

For an example on how to do RICH unfolding with Hanna++ see the Hanna++ Example Codes

For a basic example code on reading in the Pmatrix have a look at this example code and the accompanying HOWTO file