Smearing

From EIC
Revision as of 13:45, 16 September 2013 by Tpb (talk | contribs) (→‎Generic Tracker: Defined N in generic tracking formula)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

The eic-smear package provides a means for applying detector smearing to Monte Carlo (MC) data, whereby kinematic properties from the MC particles are modified according to a description of detector performance. The package contains a series of classes and functions to facilitate the smearing of ROOT trees created from Monte Carlo output using BuildTree. It is intended to be extremely flexible, allowing you to change what type of smearing you want to simulate quickly and painlessly.

Overview

Smearing is designed to allow rapid, approximate estimates of the effect of detector acceptance and performance on physics observables. As implemented in eic-smear, it is much faster to run than a full detector simulation (e.g. Geant), both in terms of user time spent specifying the detector properties and in CPU time spent processing events. For example, the smearing code typically processes hundreds to thousands of events per second on a modern CPU, which is comparable to the generation rate of many of the event generators supported. Compare this to the performance of a detailed Geant simulation, where specifying the detector properties can be a non-trivial task, and events may each take many seconds to be processed. That said, the smearing code is "not" intended as a replacement for full detector simulations. Rather, it is a complimentary tool for roughly assessing the impact of changes in detector performance on a physics measurement, in a way that is much more rapid, but less detailed, than a full detector simulation.

Code structure

Smearing is managed via a Detector object, which is composed of one or more objects inheriting from the Smearer class. A number of classes inheriting from Smearer are implemented (Device, ParticleID etc), and the user is free to define their own as needed. Each Smearer defines some detector behaviour, optionally with an acceptance in one or more particle properties, for example angle, momentum or energy. A Smearer should not necessarily be considered to be a direct analogue of a single real-world detector subsystem (though in some cases it may be). Rather, it is a more abstract representation of some element of the whole detector's performance. For example, while a user may define a single Device object for smearing momentum, this Device may represent the combined tracking performance and acceptance of a number of physical detector subsystems.

The Detector applies the smearing of all its component Smearers to an input MC event, and generates a corresponding smeared event. File I/O is handled by the function SmearTree.

The sections below describe the main functions and classes the user needs to work with to perform smearing. All classes described below are contained within the Smear namespace to prevent name clashes with other code. Class documentation can be found here.

Extensibility

The smearing code is designed to interface directly with the output from Monte Carlo events built via the BuildTree routine in eic-smear. However, users can apply smearing to their own event structures, so long as they inherit from the necessary eic-smear base classes. Developers should consult the class documentation for details for the eic-smear API.

The SmearTree function

SmearTree drives processing of an input MC tree and generates a corresponding output smeared tree. Note that, unlike the classes detailed below, the SmearTree function is not a member of the Smear namespace. It takes up to four arguments:

  1. A Detector object, defining the detector behaviour.
  2. The name of an input MC tree to smear.
  3. The name of the output smeared file.
  4. The maximum number of events.

The last two are optional. By default all input events are processed and the output name is determined automatically: if the original file name is "filename.root", the smeared file name is "filename.smear.root". Hence the SmearTree function can be called via:

SmearTree(detector, "filename.root");                        // Use default output name and all events
SmearTree(detector, "filename.root", "outname.root", 10000); // Generate outname.root with up to 10,000 events

The Detector Class

The Detector class is a container for Smearer objects. When it's Smear() method is called on an MC particle, it performs the smearing of all its component Smearers on that particle. The results are stored in a smeared particle object. The user should assemble a Detector defining the smearing they wish to perform, and pass it to the SmearTree function. Smearer objects are added to a Detector via AddDevice().

Importantly, it only makes sense to add one smearer for each type of kinematic variable in a given acceptance range; multiple Smearers for a given variable whose acceptance overlap do not stack - only one will be applied, and which one will be applied is undefined. Note that it is perfectly acceptable to have multiple Smearers affecting the same variable so long as their acceptances do not overlap. Also note that here "acceptance" refers not just to geometrical acceptance, but to such properties as particle energy, charge and species i.e. anything that can be specified to restrict the particles to which a Smearer is applied.

The Detector also calculates smeared event-wise kinematics for DIS events (namely x, y, Q2 and W2), using up to 3 different methods of calculation.

  • The null momentum (NM) method uses the scattered electron to calculate event-wise variables in the null momentum (m → 0) approximation.
  • The Jacquet-Blondel (JB) calculates kinematics using the particles in the hadronic final state.
  • The double-angle (DA) method combines information from the scattered electron and the hadronic final state.

These are selected via:

mydetector.SetEventKinematicsCalculator("NM");       // Selects only one
mydetector.SetEventKinematicsCalculator("JB DA NM"); // Selects all three, order does not matter

By default no kinematic calculations are performed on smeared events.

The Device Class

This is a simple class that smears exactly one kinematic variable, x, according to a formula describing the Gaussian resolution, σ(x), in that variable. Note that the formula should provide the absolute resolution, not the relative resolution σ(x)/x. Supported variables for smearing are: energy, momentum (total, transverse or longitudinal) and angle (polar or azimuthal). Which one to smear is passed to the Device constructor via a (case-sensitive) string:

"E", "P", "pT", "pZ", "theta", "phi"

or one of the following enumerations:

enum Smear::KinType { kE, kP, kTheta, kPhi, kPt, kPz };

The resolution formula is defined via a string, containing some combination of mathematical operations and/or variables chosen from the strings listed above. Up to four variables can be included in the formula defining resolution. For example, to smear energy according to a resolution of σ(E) = 12% sqrt(E) + 0.05 (i.e. σ(E)/E = 12% / sqrt(E) + 0.05 / E):

Smear::Device energy("E", "0.12 * sqrt(E) + 0.05");       // or equivalently…
Smear::Device energy(Smear::kE, "0.12 * sqrt(E) + 0.05");

To smear polar angle according to σ(P) = 1% / (P * sqrt(sin(θ))):

Smear::Device theta(Smear::kTheta, "0.01 / (P * sqrt(sin(theta)))");

An optional third argument defines the "genre" of particle to smear:

  • 0: All particles (default).
  • 1: Electromagnetic (photon and electron/positron only).
  • 2: Hadronic (pion, kaon, proton etc).

For example:

Smear::Device emcal(Smear::kE, "0.12 * sqrt(E) + 0.05", 1); // Only smears photon, electrons, positrons

Generic Tracker

The Tracker class implements a generic function for momentum resolution in a solenoidal field. The total resolution contains two contributions: the intrinsic detector resolution, and a multiple scattering component. The intrinsic resolution (due to the detector point resolution) is given by

[math]\displaystyle{ \left\vert \frac{\sigma(p)}{p} \right\vert _{intrinsic} = \frac{p\sigma}{0.3BL_{T}^{2}}\sqrt{\frac{720}{N+4}} }[/math]

for a track of momentum [math]\displaystyle{ p }[/math] following a transverse path length [math]\displaystyle{ L_{T} }[/math] through a magnetic field [math]\displaystyle{ B }[/math] in a detector with point resolution [math]\displaystyle{ \sigma }[/math]. [math]\displaystyle{ N }[/math] is the number of measurement points along the track. The factor 720 is replaced by 320 in the case that there is a vertex constraint. The resolution due to multiple scattering is

[math]\displaystyle{ \left\vert \frac{\sigma(p)}{p} \right\vert_{scattering} = \frac{0.016\sqrt{R}}{0.3BL_{T}\beta} }[/math]

for a track passing through [math]\displaystyle{ R }[/math] radiation lengths of material. The two contributions are added in quadrature. Strictly, these formulae are derived for transverse (not total) momentum, for tracks perpendicular to the beam. However we assume that the effects of angular resolution are small and that the formulae apply to total momentum in that case. See the references at the end of this section for further reading.

The Tracker describes a hollow cylindrical volume, defined via inner and outer radii and the z positions of the two end planes. These variables, along with magnetic field strength, material budget (in number of radiation lengths), number of fit points and the intrinsic point resolution, are all settable via the constructor. The angular acceptance is automatically determined via the geometrical parameters, but can also be set manually if the user wishes to restrict acceptance to a region smaller than the physical dimensions..

Users do not directly create a Tracker object. Instead the code provides for two types of tracking volumes: radial, via RadialTracker, and planar, via PlanerTracker. A schematic representation of each is shown below.

Tracking types.png

Both volumes perform the smearing on momentum with contributions from intrinsic resolution and multiple scattering, but they differ in the way the number of fit points n is calculated for each track that lies within the acceptance. The radial tracker mimics the behaviour of a drift chamber or TPC, and n is assigned as a fraction of the maximum possible fit points based on the track's transverse path through volume, rounded to the nearest integer. The planar tracker emulates vertical layers, equally-spaced along the z-axis, and n is the number of planes the particle would intersect.

For example, to simulate a TPC with R = 20 cm to R = 40 cm, with total length of 1 m (extending from z = -50 cm to z = 50 cm), in a 2 Tesla magnetic field with 0.03 radiation lengths, 80 micron point resolution and a maximum of 40 fit points:

Smear::RadialTracker tracker1(0.2, 0.4, -0.5, 0.5, 2., 0.03, 80.e-6, 40);

Similarly, to declare a region of 6 equally spaced tracking layers with inner radius R = 2 cm to accommodate the beam pipe and outer radius R = 1 m, extending from z = 50 cm to z = 100 cm:

Smear::PlanarTracker tracker2(0.02, 1.0, 0.5, 1.0, 2., 0.03, 80.e-6, 6);

If not specified separately, the angular acceptance of a tracking volume is automatically calculated from geometric parameters. If both of these volumes were created in the same detector script, their acceptance regions would overlap. As it stands, the code does not allow for a particular kinematic variable (say, p) to be smeared by more than one smearing device. To assign each tracker to a specific region in angular acceptance:

tracker1.Accept.AddZone(central);
tracker2.Accept.AddZone(forward);

where "central" and "forward" are acceptance zones defined previously in the script (see section on Acceptance). With these angular acceptance zones implemented, the example layout described above ("tracker1" +"tracker2") functions effectively like the drawing below. Some results from this example layout are shown here.

Example.png

References:

The ParticleID Class

The ParticleID class is a Smearer used to select particle species using a pre-calculated probability matrix stored in a plain text file. It does not smear any kinematic properties, but instead can be used to describe particle (mis)identification probabilities by a detector. It reads a table of identification matrices that describe the probability for a particle of some type to be identified as another type. By default the ParticleID will use smeared values of momentum to make its determinations, but the original MC values can be used instead instead by using

ParticleID::SetPIDUseMC(true);

The eic-smear distribution contains an example matrix table, "PIDMatrix.dat", in the data/ subdirectory, which describes the performance of the HERMES RICH detector. The name of the file should be passed to the ParticleID constructor. For example:

Smear::ParticleID hermes("/path/to/PIDMatrix.dat");

The user can define their own custom (mis)identification matrices by following the format of PIDMatrix.dat.

PID Matrix File Format

The file contains one matrix for each of n momentum bins. The file must contain the following 3 formatting lines, before any data it contains:

!T true1 true2 true3
!F false1 false2 false3 false4
!P Nbins

Here "true1", "true2" and "true3" are the PDG codes of 3 particles which the ParticleID will attempt to identify. "false1" to "false4" are the PDG codes of 4 particles which the ParticleID may identify one of "trueX" as. These lists should contain positive numbers; the ParticleID will identify both that particle and the associated anti-particle if it exists. Lastly, Nbins is the number of momentum bins. If you do not want to use one of the true or false slots, you should enter a bogus (invalid as a PDG code) value such as -999. For example, a typical formatting header would be

!T 211 321 -999
!F 211 321 2212 0
!P 50

A ParticleID with the above formatting will attempt to identify pions and kaons. Each pion and kaon it sees may be identified as a pion, a kaon, a proton or a 0, which the user can interpret as a failure to identify the particle. It also has 50 momentum bins. Lines of data in the P matrix file should take the form

1 pbinNumber pMin pMax FalseIndex element1 element2 element3

Every data line must begin with the number 1. pbinNumber is the number of the momentum bin, which cannot be less than 1 or more than Nbins. pMin is the lower bound of the momentum bin, pMax is the upper bound of the momentum bin (should be the same for every line with the same number bin), FalseID is the PDG code of the particle of which the identity will be assigned. elementX are the probability to identify true1, true2, true3 as FalseID. For example a couple of lines of a typical P matrix data file would be

1 4 4.0 5.0 2 0.05 0.94 0.
1 4 4.0 5.0 3 0.04 0.06 0.

With the example formatting lines above, the above data lines specify that a pion with a momentum 4 < p <5 GeV has a 0.05 probability of being identified as a kaon, and a 0.04 probability of being identified as a proton. A kaon with a momentum 4 < p < 5 GeV has a probability of 0.94 of being identified as a kaon, and a probability of 0.06 of being identified as a proton. The 0's are placeholders to make sure the reader doesn't get confused (three probability columns are required). The Smear::Detector can contain arbitrarily many ParticleID instances, so if you want to identify more than 3 particles, you should use more than one ParticleID instance, with the appropriate acceptances. For example you may have one which identifies hadrons, and another which identifies leptons.

PerfectID

By default smeared particles have no species assignment (PDG code). PerfectID copies the ID of an input particle in its acceptance to the smeared counterpart without modification.

Acceptance

Every Smearer has an instance of the Acceptance class, named Accept, describing a region of (E, p, θ, φ) associated with it. It will only smear variables from the particles which fall within the device acceptance. By default the acceptance allows all energies, momenta and angles, so the user need not explicitly set the acceptance if they wish a Smearer to accept everything. An Acceptance contains one or more Zones, where each Zone describes a continuous (E, p, θ, φ) region. The Acceptance can describe (possibly non-adjacent) regions of (E, p, θ, φ) by adding multiple Zones:

Smear::Acceptance::Zone central(TMath::Pi() / 4, 3 * TMath::Pi() / 4); // Defines the angular range [pi/4, 3pi/4]
Smear::Device emcal(Smear::kE, "0.12 * sqrt(E) + 0.05", 1);
emcal.Accept.AddZone(central);

The syntax reflects the fact that Zone is a nested class of Acceptance. By default all stable particles will be smeared. You might want to have a highly specialized device, which only smears certain particles. For example, for a Device called dev:

dev.Accept.AddParticle(211);  // now detects only pi+
dev.Accept.AddParticle(-211); // now detects pi+ and pi-
dev.Accept.AddParticle(321);  // now detects pi+, pi- and K+

For more information, see the Acceptance class documentation.

Data Format

When you use the SmearTree function a new root file will be generated. The file contains a TTree, named "Smeared", with a single branch, named "eventS". Each entry in the branch is an instance of an event class named EventSmear, which is essentially a stripped down version of the Monte Carlo events produced by BuildTree. EventSmear contains a list of particles of class ParticleMCS ("particle Monte Carlo smeared"), which is essentially a stripped down version of ParticleMC used by BuildTree. A different class is used to save disk space, as the Monte Carlo particles contain a large amount of information which is not required (or appropriate) for the smeared particles.

Smeared data is stored using a simple scheme which matches smeared data to MC data by index. The nth MC event in a tree corresponds to the nth smeared event, and the nth particle in an MC event corresponds to the nth particle in the corresponding smeared event. This is to avoid the need for a complicated indexing scheme matching MC and smeared particles. Many particles stored in the MC tree are not smeared, due either to not being final-state particles, or not falling within the acceptance of any Smearer. Such particles are stored in the smeared tree as NULL pointers and hence occupy almost zero space.

To simultaneously navigate the MC and smeared trees, one can use the ROOT "friendship" mechanism for trees:

root[1] TTree tree("EICTree", "filename.root")
root[2] tree.AddFriend("Smeared", "filename.smear.root")

One can then use commands such as tree.Draw("Smeared.particles.p:EICTree.particles.p") to plot smeared momentum versus momentum. One can also use

tree.SetBranchAddress("event", &event);
tree.SetBranchAddress("eventS", &eventS);
tree.GetEntry(0);

Each call to GetEntry() will then read the events from both trees simultaneously.

Example Scripts

Below is an example of a script assembling a Detector:

Smear::Detector BuildExample() {
   gSystem->Load("libeicsmear");
   // Create an electromagnetic calorimeter with sigma(E) = 20% * sqrt(E)
   // Genre == 1 (third argument) means only photons and electrons are smeared
   Smear::Device emcal(Smear::kE,  "0.2*sqrt(E)", 1);

   // Create our tracking capabilities, by a combination of mometum, theta and phi Devices.

   // Create a smearer for momentum
   Smear::Device momentum(Smear::kP, "0.01 * P");
   Smear::Device theta(Smear::kTheta, "0.05 * P");
   Smear::Device phi(Smear::kPhi, "0"); // "0" indicates perfect performance i.e. sigma(phi) = 0

   // Create a based on Hermes RICH.
   Smear::ParticleID rich("PIDMatrix.dat");

   // Our detector only covers the central region
   Smear::Acceptance::Zone central(TMath::Pi() / 4., 3. * TMath::Pi() / 4.);

   emcal.Accept.AddZone(central);
   momentum.Accept.AddZone(central);
   theta.Accept.AddZone(central);
   phi.Accept.AddZone(central);
   rich.Accept.AddZone(central);

   // Create a detector and add the devices
   Smear::Detector det;
   det.AddDevice(emcal);
   det.AddDevice(momentum);
   det.AddDevice(theta);
   det.AddDevice(phi);
   det.AddDevice(rich);
   det.SetEventKinematicsCalculator("NM JB DA"); // The detector will calculate event kinematics from smeared values

   return det;
}

To run this in ROOT (assuming libeicsmear is already loaded):

root[0] .L ExampleScript.cpp
root[1] SmearTree(BuildExample(), "mcFile.root");

For those so inclined, it can be more convenient to make use of the ROOT Python bindings to do this in a Python script:

import ROOT
import math

def build_example():
    """Example for building a simple Smearing Detector."""
    ROOT.gSystem.Load('libeicsmear')

    Device = ROOT.Smear.Device

    # Create the devices, which all have the same acceptance in this example.
    # - an electromagnetic calorimeter with sigma(E) = 20% * sqrt(E).
    # - tracking capabilities, by a combination of mometum, theta and phi.
    # - a PID detector based on Hermes RICH.
    devices = [
        Device(ROOT.Smear.kE, '0.2*sqrt(E)', 1), # Genre 1 -> photon, electron
        Device(ROOT.Smear.kP, '0.01 * P'),
        Device(ROOT.Smear.kTheta, '0.05 * P'),
        Device(ROOT.Smear.kPhi, '0'), # '0' indicates perfect performance
        ROOT.Smear.ParticleID('PIDMatrix.dat')
    ]

    # Our detector only covers the central region
    central = ROOT.Smear.Acceptance.Zone(math.pi / 4., 3. * math.pi / 4.)

    det = ROOT.Smear.Detector()
    det.SetEventKinematicsCalculator('NM JB DA')
    for i in devices:
        i.Accept.AddZone(central)
        det.AddDevice(i)

    return det

Information Welcome!

We are in the process of gathering information for testing and implementation of smearing, so if you would like to see some aspect of your group's detector included, we will be pleased to add it into one of the available layouts. You should contact Thomas Burton or Elke Aschenauer to provide parametrizations, or whatever else.