Second Tracking Tutorial

Goal of this exercise

This exercise follows up on SWGuideCMSDataAnalysisSchoolIntroductoryTrackingExercise by making the connection between fundamental track distributions and how to use them in a physics analysis. At the end of the exercise, we will all make invariant mass distributions from real 2010 data and search for resonance peaks.

What I hope you will learn from this exercise are:

• how to learn what quantities are available in a dataset for your analysis ("shopping for things to plot");
• starting points to learn more about using tracks;
• working with the geometry of tracks to build new information (we will be calculating vertices from pairs of tracks);
• and finally, how to plot invariant mass distributions and search for resonances.

Contacts

Please contact Jim Pivarski for any questions, comments or suggestions about the tutorial.

Preliminaries

# enable CMSSW commands, if you haven't already
source /uscmst1/prod/sw/cms/cshrc prod

# create a CMSSW project area (version 3_9_7)
scramv1 project CMSSW CMSSW_3_9_7
cd CMSSW_3_9_7/src/
cmsenv

Exploring tracking data

In this section, I will walk you through a very common first task: finding out what data records are available for your analysis. Once upon a time, this was done through a time-consuming cycle of (1) searching documentation, (2) writing code, (3) compiling code, (4) running code, (5) looking at the error message and wondering what went wrong (back to 1). The CMSSW framework provides alternate access methods to the data that are faster for the user but slower for the computer--- allowing us to quickly figure out what we'll be able to do before we write our large-scale analysis code. That is, we sketch before we paint.

We'll start by running a Python/FWLite script in your CMSSW project directory. Copy the following into a file called explore.py:

# Contents of the explore.py script
import ROOT
from DataFormats.FWLite import Handle, Events

# point to data file
events = Events("dcap://cmsdca3.fnal.gov:24144/pnfs/fnal.gov/usr/cms/WAX/11/store/user/pivarski/TracksAndMuons/2010B/events_1_1_7kI.root")

# define (reco) track collection type and (pat) muon collection type
tracks = Handle("std::vector<reco::Track>")
muons = Handle("std::vector<pat::Muon>")

# start to loop over events, but break during the first event
for event in events:
    # get data from the event, using their labels
    event.getByLabel("generalTracks", tracks)
    event.getByLabel("cleanPatMuons", muons)

    # print out a lot of information
    for track in tracks.product():
        print "track", track.pt(), track.eta(), track.phi()
    for muon in muons.product():
        print "muon", muon.pt(), muon.eta(), muon.phi()

    # break from the loop (so that we're still in the first
    #     event when Python gives you interactive access)
    break

and run the script with

python -i explore.py

It's important to run cmsenv (see Preliminaries) before Python, and to include the -i option (so that Python stays alive and gives you a prompt after it's done with explore.py).

The script reads in one event from the data file, loads all tracks and muons from that event, and prints out the pT, eta, and phi of each. You'll notice that there are a lot more "tracks" than "muons."

Interacting with tracks

Now, at the Python prompt (>>>), type the following to find out what kinds of objects we're looking at.

>>> tracks
<DataFormats.FWLite.Handle instance at 0xf59de94c>   # a Handle--- not surprising; everything delivered by CMSSW is in a Handle of some sort
>>> tracks.product()
<ROOT.vector<reco::Track> object at 0xc3fa590>   # a vector (list) of reco::Track objects
>>> tracks.product().size()
166L

>>> track = tracks.product()[0]     # get the first track (could just as easily get index 165)
>>> track
<ROOT.reco::Track object at 0xc4cc940>   # a reco::Track--- what do reco::Tracks have for us to look at?

>>> track.pt()
0.70336815331205838

>>> track.numberOfValidHits()
20

>>> track.dxy(), track.dxyErr()      # oops, I mis-spelled the variable name for the dxy uncertainty
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'reco::Track' object has no attribute 'dxyErr'

>>> track.dxy(), track.dxyError()
(-0.068770636200316079, 0.019096558158866263)     # better

How did I know what to type?

You've already been introduced to the main sources of CMSSW documentation:

• LXR (for finding an object's definition and examples of its use)
• CVS repository (for tracking changes between versions)
• doxygen (for following relationships between superclasses and subclasses)
• Workbook (human-written tutorials)
• SWGuide (human-written reference)

Information about tracks (specifically, reco::Track objects) can be found here:

• TrackBase.h
• Track.h

(reco::Track is a subclass of reco::TrackBase, but most of the variables we frequently use are in reco::TrackBase.)

Since the code knows what data are available, we could just ask it directly:

>>> dir(track)

The above should give you a list of all the member data and member functions in reco::Track objects, though without any information about what they do. (Cross-reference this output with the documentation above.)

Let's try asking for some detailed information, like the x residual (difference between the expected track position and the actual hit position) for the first hit:

>>> track.residualX(0)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
RuntimeError: (file "", line 0) ---- ProductNotFound BEGIN
RefCore: A request to resolve a reference to a product of type: St6vectorIN4reco10TrackExtraESaIS1_EEwith ProductID 2:769
can not be satisfied because the product cannot be found.
Probably the branch containing the product is not stored in the input file.
---- ProductNotFound END
 (C++ exception)

Why was there an error? In this case, it's because the data sample that I prepared for today's exercise has only the minimum information for our objective: making invariant mass plots. The residualX call references a separate data structure, reco::TrackExtra, which contains information about the track's trajectory at different points along the path. This extra information was deliberately put into a separate object so that it can be voluntarily included or excluded from custom datasets. At some point in your analysis, you may be asked whether you need information in RECO (the full set of reconstructed quantities) or just AOD (a subset designed for physics analyses). You can use the above to quickly check for your favorite quantities!

Interacting with muons

The data sample that we're using for this exercise contains only tracks and muons. As a particle signature, a muon is primarily just a track (as opposed to an electron, for which the calorimeter signature is essential). One could think of muons as being special tracks that have hits in the muon chambers as well as the silicon tracker, but muons are represented by a class in CMSSW, much like electrons, photons, jets, missing energy, etc.

In this exercise we will be using pat::Muons, rather than reco::Muons, though there is no difference for the user. Briefly, PAT is a collection of tools for working with high-level physics objects, but any call that can be made from a reco::Muon can also be made from a pat::Muon. The documentation can be found in

• DataFormats/MuonReco/interface/Muon.h
• DataFormats/PatCandidates/interface/Muon.h (wrapper around reco::Muon)
• DataFormats/PatCandidates/interface/PATObject.h (superclass of pat::Muon)

Much like the track example above, we can ask muons about the data they store:

>>> muon1 = muons.product()[0]
>>> muon2 = muons.product()[1]
>>> muon1.pt(), muon2.pt()
(11.123780250549316, 4.9288649559020996)

One particularly important feature of a muon is that it contains up to three tracks: an inner track (tracker-track), an outer, muon-system only track (standAloneMuon), and a combined track (GlobalMuon). The three tracks represent the same particle, but using different subdetectors with different resolutions, and therefore yield slightly different results.

>>> muon1.innerTrack().pt(), muon1.outerTrack().pt(), muon1.globalTrack().pt()
(11.123780617918049, 12.702820606346632, 11.143488366078506)

Depending on how the muon was constructed, not all of these will be available in every muon, so be sure to check isAvailable() in your code.

>>> muon1.innerTrack().isAvailable()
1
>>> muon1.outerTrack().isAvailable()
1
>>> muon1.globalTrack().isAvailable()
1

Though I said that a muon is essentially a track, the muon object contains more useful information than just the tracks, such as the muon's small energy deposits in the calorimeters. See the documentation and tinker with it to learn more.

Incidentally, we are now at a stage where we can calculate our first invariant mass:

>>> from math import sqrt
>>> sqrt((muon1.energy() + muon2.energy())**2 - (muon1.px() + muon2.px())**2 - \
    (muon1.py() + muon2.py())**2 - (muon1.pz() + muon2.pz())**2)
17.844641987474031

Getting familiar with the dataset

The data sample comes from the "Mu" primary dataset (/Mu/Run2010B-Dec22ReReco_v1/RECO), which means that these events were collected by muon triggers. To make this skim, I applied the following cuts:

process.selectedPatMuons.cut = cms.string("pt > 3.0 && \
    ((isTrackerMuon() && numberOfMatches() >= 2 && -2.4 < eta() && eta() < 2.4) || \
    isGlobalMuon())")
process.countPatMuons.minNumber = cms.uint32(2)

All "cleanPatMuons" in the final sample must have p_T > 3 GeV/c, and they must either be clean TrackerMuons (identified by extrapolating the tracker's tracks to the muon system) or GlobalMuons (identified by extrapolating the muon system's tracks to the tracker). Also, all events in the final sample must have at least two "cleanPatMuons".

To see this for yourself, try replacing the loops and break statement in explore.py with a print-out of the number of muons in each event

    print muons.product().size(),
#     for track in tracks.product():
#         print "track", track.pt(), track.eta(), track.phi()
#     for muon in muons.product():
#         print "muon", muon.pt(), muon.eta(), muon.phi()
#    break

and run it from a new instance of Python. You should see a lot of "2"s, a few "3"s, and even more rarely, "4".

To make this observation explicit, re-arrange explore.py to look like the following:

from DataFormats.FWLite import Handle, Events
from math import sqrt
import ROOT    # for plotting histograms

events = Events("dcap://cmsdca3.fnal.gov:24144/pnfs/fnal.gov/usr/cms/WAX/11/store/user/pivarski/TracksAndMuons/2010B/events_1_1_7kI.root")
muons = Handle("std::vector<pat::Muon>")

histogram_number = ROOT.TH1F("histogram_number", "", 6, -0.5, 5.5)
histogram_pt = ROOT.TH1F("histogram_pt", "", 100, 0., 50.)
histogram_mass = ROOT.TH1F("histogram_mass", "", 100, 0., 5.)

for event in events:
    event.getByLabel("cleanPatMuons", muons)
    histogram_number.Fill(muons.product().size())

    for muon in muons.product():
        histogram_pt.Fill(muon.pt())

    if muons.product().size() == 2:
        muon1 = muons.product()[0]
        muon2 = muons.product()[1]

        histogram_mass.Fill(sqrt((muon1.energy() + muon2.energy())**2 - (muon1.px() + muon2.px())**2 - \
                                 (muon1.py() + muon2.py())**2 - (muon1.pz() + muon2.pz())**2))

Run it (with the -i option), and then call the following on the Python prompt.

>>> histogram_number.Draw()
>>> histogram_pt.Draw()
>>> histogram_mass.Draw()

You should see the following (click to enlarge):

We can already see the resonance at an invariant mass of 3.1 GeV/c² (without compiling anything)! We could continue like this, except for two problems:

• the Python/FWLite interface is too slow to process large datasets;
• some operations are not possible in FWLite, such as calculating intersections (vertices) of tracks.

Calculating vertices

A proton-proton collision produces about a hundred tracks (though most of them have p_T < 1.5 GeV/c). If we tried to calculate the invariant masses of all two-track pairs, we'd produce about ten thousand candidates per event, and the vast majority of them would not be real resonances. Our plots would be swamped with backgrounds. In any bump-hunt, one must reduce the number of combinations in an intelligent way: the example above used the fact that very few tracks are muons; in this section, we'll use the fact that tracks are one-dimensional curves in a three-dimensional space--- in general, not guaranteed to intersect.

Vertex-calculations cannot be done in FWLite because they require access to edm::EventSetup (constants, geometry, magnetic field), so we'll build a module in C++. Although we need to switch languages, all of the class structures, member data, member functions, etc. are exactly as they were in the Python examples above. The biggest difference is that we need to start worrying about when to use "." and when to use "->" (pointers in C++).

Start by creating a blank CMSSW analyzer in a Test directory:

mkdir Test
cd Test/
mkedanlzr MyAnalyzer
cd ..
emacs Test/MyAnalyzer/src/MyAnalyzer.cc &      # (or use your favorite editor)

and make the following modifications.

(1) Add these include files directly below "#include "FWCore/ParameterSet/interface/ParameterSet.h":

// for tracks
#include "DataFormats/TrackReco/interface/Track.h"

// for vertexing
#include "FWCore/Framework/interface/ESHandle.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"
#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"

// for making histograms
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "TH1F.h"
#include "TH2F.h"

(How did I know which ones? A series of syntax and runtime errors, each error message telling me what's missing.)

(2) Add the following packages at the top of your Test/MyAnalyzer/BuildFile:

<use name=DataFormats/TrackReco>
<use name=TrackingTools/TransientTrack>
<use name=TrackingTools/Records>
<use name=RecoVertex/KalmanVertexFit>
<use name=FWCore/ServiceRegistry>
<use name=CommonTools/UtilAlgos>
<use name=root>

(Did you notice that these are all just the first two words of each #include (except for TH1F.h and TH2F.h, which were converted to root)? Converting includes to BuildFile entries like this is always sufficient to get an analyzer working.)

(3) Delcare the histograms directly below the "// ----------member data" comment:

TH2F *m_histogram_xy;
TH1F *m_histogram_z;
TH1F *m_histogram_mass_all;
TH1F *m_histogram_mass_displaced;

(4) Book the histograms in the MyAnalyzer::MyAnalyzer constructor:

{
//now do what ever initialization is needed
edm::Service<TFileService> tfile;
m_histogram_xy = tfile->make<TH2F>("histogram_xy", "", 300, -3., 3., 300, -3., 3.);
m_histogram_z = tfile->make<TH1F>("histogram_z", "", 30, -30., 30.);
m_histogram_mass_all = tfile->make<TH1F>("histogram_mass_all", "", 100, 0.4, 0.6);
m_histogram_mass_displaced = tfile->make<TH1F>("histogram_mass_displaced", "", 100, 0.4, 0.6);
}

(5) Replace the contents of the MyAnalyzer::analyze method with the following:

{
   edm::Handle<std::vector<reco::Track> > tracks;
   iEvent.getByLabel("generalTracks", tracks);

   // this wraps tracks with additional methods that are used in vertex-calculation
   edm::ESHandle<TransientTrackBuilder> transientTrackBuilder;
   iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", transientTrackBuilder);

   // loop over unique pairs of opposite-sign tracks, "one" and "two"
   for (std::vector<reco::Track>::const_iterator one = tracks->begin();  one != tracks->end();  ++one) {
      for (std::vector<reco::Track>::const_iterator two = one;  two != tracks->end();  ++two) {
         if (two != one  &&  one->charge() != two->charge()) {
            // make a pair of TransientTracks to feed the vertexer
            std::vector<reco::TransientTrack> tracksToVertex;
            tracksToVertex.push_back(transientTrackBuilder->build(*one));
            tracksToVertex.push_back(transientTrackBuilder->build(*two));

            // try to fit these two tracks to a common vertex
            KalmanVertexFitter vertexFitter;
            CachingVertex<5> fittedVertex = vertexFitter.vertex(tracksToVertex);

            // some poor fits will simply fail to find a common vertex
            if (fittedVertex.isValid()  &&  fittedVertex.totalChiSquared() >= 0.  &&  fittedVertex.degreesOfFreedom() > 0) {
               // others we can exclude by their poor fit
               if (fittedVertex.totalChiSquared() / fittedVertex.degreesOfFreedom() < 3.) {

                  // this is a good vertex: fill histograms
                  m_histogram_xy->Fill(fittedVertex.position().x(), fittedVertex.position().y());
                  m_histogram_z->Fill(fittedVertex.position().z());

                  // important! evaluate momentum vectors AT THE VERTEX
                  TrajectoryStateClosestToPoint one_TSCP = tracksToVertex[0].trajectoryStateClosestToPoint(fittedVertex.position());
                  TrajectoryStateClosestToPoint two_TSCP = tracksToVertex[1].trajectoryStateClosestToPoint(fittedVertex.position());
                  GlobalVector one_momentum = one_TSCP.momentum();
                  GlobalVector two_momentum = two_TSCP.momentum();

                  // calculate mass
                  double total_energy = sqrt(one_momentum.mag2() + 0.14*0.14) + sqrt(two_momentum.mag2() + 0.14*0.14);
                  double total_px = one_momentum.x() + two_momentum.x();
                  double total_py = one_momentum.y() + two_momentum.y();
                  double total_pz = one_momentum.z() + two_momentum.z();
                  double mass = sqrt(pow(total_energy, 2) - pow(total_px, 2) - pow(total_py, 2) - pow(total_pz, 2));

                  // look at masses everywhere and at least 2 cm from the beamline
                  m_histogram_mass_all->Fill(mass);

                  double distance_perp = sqrt(pow(fittedVertex.position().x(), 2) + pow(fittedVertex.position().y(), 2));
                  if (distance_perp > 2.) m_histogram_mass_displaced->Fill(mass);
               }
            }
         }
      }
   }
}

Take a moment now to read through the code that you have just inserted. After extracting the tracks and TransientTrackBuilder from the event record, it performs a loop over all unique pairs of oppositely charged tracks, calculates the vertex of each, and then the invariant mass at the vertex.

(6) Create a configuration file named MyAnalyzer_cfg.py with the following contents:

import FWCore.ParameterSet.Config as cms

# set up the process and the input file list
process = cms.Process("MAKENTUPLE")
process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring(
    "/store/user/pivarski/TracksAndMuons/2010B/events_1_1_7kI.root",
    ))
process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(100))

# load system includes for vertexing
process.load("Configuration/StandardSequences/FrontierConditions_GlobalTag_cff")
process.GlobalTag.globaltag = "GR_R_39X_V5::All"
process.load("TrackingTools.TransientTrack.TransientTrackBuilder_cfi")
process.load("Configuration.StandardSequences.Geometry_cff")
process.load("Configuration.StandardSequences.MagneticField_cff")

# load our module and put it in the path
process.MyAnalyzer = cms.EDAnalyzer("MyAnalyzer")
process.Path = cms.Path(process.MyAnalyzer)

# set up a ROOT file for ntuple output
process.TFileService = cms.Service("TFileService", fileName = cms.string("MyAnalyzer.root"))

(7) Finally, compile and start cmsRun:

scramv1 build
cmsRun MyAnalyzer_cfg.py

You should see it process 100 events (slowly). When it's done, draw the plots in ROOT (e.g. with the TBrowser). They should look like the following:

There's a huge concentration of vertices at (x, y) = (0.1 cm, 0.02 cm): these are particles emerging from the beamspot. However, there are plenty of good vertices at least as far out as 3 cm. (Also note that the LHC beams don't collide exactly at (0, 0)!)

In the z projection (along the beamline), the distribution is much broader: about 6 cm. We should be thinking of the interaction region as a long, narrow strip, rather than a circular blob.

When we look at the invariant mass distribution of all good vertices, there's no hint of anything special happening at the mass, 0.498 GeV/c².

For good verticies 2 cm or more from (0, 0), however, there's a 4.5-sigma excess at the mass (from 100 events). Try fitting it to a Gaussian-plus-constant.

Hunting for resonances

In this final exercise, we will scale up our example from 100 events to 32 pb^-1 (run 2010B). It would take MyAnalyzer a very long time to process the whole data sample because it attempts to calculate a vertex for every opposite-sign pair of tracks. In the MyAnalyzer output, we saw that the vast majority of intersections occur close to the beamline, so we can increase our sensitivity to the long-lived particles, at least, by considering only vertices that are far from the beamline. Vertex-calculation is an expensive operation, so we should filter out beamline-pointing tracks before calculating a vertex: that is, require |d_xy| > 3 cm.

An analyzer called TrackResonanceNtuple scales up what we did in MyAnalyzer and writes its results to ntuples (ROOT TTrees). In this exercise, we will be analyzing the ntuples in ROOT (or PyROOT, if you prefer) to search for resonance peaks. Look over the TrackResonanceNtuple source code to understand how the ntuple variables were defined. TrackResonanceNtuple uses the same tools as MyAnalyzer, but it produces a format that will help us search for resonances quickly.

/uscms/home/pivarski/public/TrackingTutorial/TrackResonanceNtuple/src/TrackResonanceNtuple.cc

The ROOT file that we will use can be found here:

/uscms/home/pivarski/public/TrackResonanceNtuple.root

The ROOT file contains two independent TTrees, "twoTrack" (displaced verticies) and "twoMuon" (muon pairs). The "twoTrack" ntuple has the following basic features:

• track cuts: |d_xy| > 5 cm, quality="tight", number of valid hits > 7, normalized < 5;
• vertex formed with probability > 1%;
• resonance mass calculated for each of six daughter mass hypotheses.

The "twoMuon" TTree has invariant masses for all pairs of clean muons.

In ROOT, you can load and start drawing the ntuple data with

TFile *_file0 = new TFile("/uscms/home/pivarski/public/TrackResonanceNtuple.root");
TTree *twoTrack = _file0->Get("TrackResonanceNtuple/twoTrack");
TTree *twoMuon = _file0->Get("TrackResonanceNtuple/twoMuon");
twoTrack->Draw("mass_pipi", "mass_pipi < 1");

You can do the same thing by starting PyROOT (within a Python session), like this:

>>> import ROOT
>>> _file0 = ROOT.TFile("/uscms/home/pivarski/public/TrackResonanceNtuple.root")
>>> twoTrack = _file0.Get("TrackResonanceNtuple/twoTrack")
>>> twoMuon = _file0.Get("TrackResonanceNtuple/twoMuon")
>>> twoTrack.Draw("mass_pipi", "mass_pipi < 1")

For reference, "twoTrack" contains the following variables:

Variables	Brief description
mass_pipi, mass_Ppi, mass_piP, mass_KK, mass_Kpi, masspiK	invariant mass of a pair of tracks, assuming different mass hypotheses for the tracks themselves (see the exercise below for more explanation)
px, py, pz	the momentum of the pair of tracks (vector sum of the two tracks' momenta)
vertexProb	the vertex probability (uniformly distributed between 0 and 1 for intersecting track pairs, peaks at 0 for tracks that do not cross within uncertainties)
vx, vy, vz	location of the vertex (intersection of the two tracks)

and "twoMuon" contains the following:

Variables	Brief description
mass_mumu	invariant mass of the pair of muons
px, py, pz	the momentum of the dimuon (vector sum of the two muons' momenta)

The projects below give you a chance to "rediscover" a Standard Model particle with real LHC data. Unless you're very fast, don't try to do all of them--- I encourage you to pick whatever sounds interesting to you. You don't need to do all of them in the session. I've labeled each question with my own estimate of "easy," "medium," and "hard."

Searches for long-lived hadrons

According to the PDG,

Particle	Mass	Decay mode	Branching fraction	Lifetime
	0.498 GeV/c2		69.2%	2.68 cm
	1.116 GeV/c2		63.9%	7.89 cm
	1.019 GeV/c2		48.9%	4.26 MeV/c2

Can you find these particles in the data? In the "twoTrack" ntuple, each mass calculation makes a different assumption about the identities of the two tracks.

• mass_pipi assumes that they are both pions;
• mass_Ppi assumes that the highest-p_T track is a proton and the lowest-p_T track is a pion;
• mass_piP assumes the reverse;
• mass_KK assumes that they are both kaons;
• mass_Kpi assumes that the highest-p_T track is a kaon and the lowest-p_T track is a pion;
• mass_piK assumes the reverse.

Strange Vees from the 60's

Easy: Can you find the ? How about the ? What happens if you assume the wrong mass hypotheses for the tracks (and why)? What happens if you assume a different highest-p_T/lowest-p_T ordering? Why should that make a difference?

Background suppression

Medium: Now that you have found a resonance or two, can you improve the signal to background ratio using one or more of the following variables? If so, why is that? (That is, explain the geometrical reason for the improvement.)

Variable	Meaning	Ntuple variables
vertex probability	uniform distribution for real resonances, peaks at zero for fakes	vertexProb
	distance of vertex from the beamline	sqrt(vx**2 + vy**2)
	cosine of the angle between the displacement vector and the resonance momentum vector	(vx*px + vy*py)/sqrt(vx**2 + vy**2)/sqrt(px**2 + py**2) or 3-D equivalent
polar angle	angle with respect to the beamline; equivalent to pseudorapidity	pz/sqrt(px**2 + py**2)

The very low mass excess

Medium: What is the peak with invariant mass below 0.3 GeV/c² in mass_pipi? Hint: try selecting these events by mass and look at the distribution of vertices vx, vy, vz, preferably in two dimensions.

You can make 2-D plots in ROOT by

>>> twoTrack.Draw("vy:vx", "mass_pipi < 10")     # vy-versus-vx (vx is on the horizontal axis)

Why do we see a ? (Or do we?)

Hard: When you plot

>>> twoTrack.Draw("mass_KK", "0.6 < mass_KK && mass_KK < 1.5")

can you explain the peak that you see? Is it the ? Note that the 's lifetime is so short that it is measured in MeV/c², rather than cm: it is a prompt resonance. How can particles be found so far from the collision? Hint: consider using mass_pipi and mass_KK together to explain what's going on.

And finally, what is that???

Very hard (because I haven't figured it out yet): What is the peak at a mass_pipi of about 1.28 GeV/c²? Can you use the tools introduced above to figure out what it is?

Search for dimuon resonances

According to the PDG,

Particle	Mass	Decay mode	Branching fraction
	91.2 GeV/c2		3.37%
	9.46, 10.02, 10.36 GeV/c2		2.48, 1.93, 2.18%
,	3.096, 3.686 GeV/c2		5.93, 0.77%
	1.019 GeV/c2
	0.783 GeV/c2
	0.775 with = 0.149 GeV/c2
	0.548 GeV/c2

Two lumps or three?

Easy: Can you find the family? The muon system is composed of two major parts: the barrel and the endcap (see CMS quarter-view). The mass resolution is better in the barrel than it is in the endcap. Can you improve the separation between the peaks by selecting barrel-only?

Intrinsic width and resolution

Medium (because it requires ROOT knowledge for fitting): Can you find the ? If you're familiar with ROOT, try fitting the peak. Is it a pure Gaussian? Try fitting it to a Breit-Wigner distribution, or better yet, a Voigt distribution.

Scraping the bottom of the dataset

Hard (because it requires ROOT knowledge for fitting and familiarity with statistical analysis): Can you find the or the resonances, and provide statistical evidence proving that you have?

-- JimPivarski - 18-Jan-2011