Mattione GlueX Kinematic Fitting

From GlueXWiki
Revision as of 16:10, 9 February 2016 by Pmatt (Talk | contribs) (Manual Kinematic Fit Example)

Jump to: navigation, search

Quick Start

  • The kinematic fitter (DKinFitter) is in the KINFITTER library, and is fully integrated into the analysis library (via DKinFitUtils_GlueX and the DKinFitResults_factory) is fully integrated into the analysis library. So, to kinematically fit a DReaction, you only have to specify the desired kinematic fit type in the DReaction.
    • See GlueX Analysis Factories for details on the analysis libraries.
    • See sim-recon/src/programs/Analysis/plugins/b1pi_hists for an example of using the kinematic fitter in an analysis.
  • Details on the kinematic fitting equations and constraints are in GlueX-doc-2112.

Not Yet Supported

  • Spacetime constraints
  • Fits with massive neutral showers.

Constraints and Unknowns

  • Up to one particle may be missing: 3 unknowns (px, py, pz)
  • Four-momentum constraint: 4 constraints (Initial p4 - Final p4 = 0)
  • Mass constraints: 1 constraint per particle (E2 - p2 - m2 = 0)
    • Note that the ANALYSIS library does not constrain the masses of the ω or φ mesons since their widths are non-negligible.
    • Mass constraints are automatically applied when a "P4" fit is chosen.
  • Vertex constraints: 3 unknowns per vertex (x, y, z), 2 constraints per constraining particle
    • Constraining particles: Tracks, beam (when uncertainties are defined), decaying particles (when their positions are defined by a different vertex constraint)
    • E.g. for no magnetic field, p x (Δx) = 0: The path (Δx) and momentum vectors must be parallel
      • Δx = xMeasured - xUnknown
  • Spacetime constraints (NOT YET SUPPORTED): 4 unknowns per spacetime point (x, y, z, t), 3 constraints per constraining particle (vertex + time)
    • Also, one constraint per constraining neutral shower (time)

Notes for those manually performing their own fits


  • These utilities are used to setup the kinematic fit that is indicated by the DReaction::dKinFitType variable. This variable is an enum type called DKinFitType, which is defined in libraries/ANALYSIS/DKinFitResults.h
  • These utilities automatically determines what p4, vertex, and spacetime (not yet supported!) constraints to setup such that each the data in each DParticleCombo are maximally constrained.
    • It will also automatically create additional vertex constraints for decaying particles with long lifetimes, as long as the IsDetachedVertex() function in libraries/include/particleType.h is kept updated.
    • If a missing or decaying particle of type Unknown is present in the DReaction, the fit will not constrain the overall p4 of the reaction, but it will still apply the other constraints.
  • To exclude decaying particles from the mass constraints of a kinematic fit, call the DReactionStep::Set_KinFitConstrainInitMassFlag(bool) function with an argument of true on the step where the decaying particle is in the initial state.

Detected Neutral Particles

  • The DNeutralParticleHypothesis objects from the default factory sets the particle momentum assuming that the particle came from the position defined by the DVertex (kinematic fit of all tracks).
  • The best way to accurately perform kinematic fits with neutral particles is to perform a simultaneous vertex & p4 fit, with the neutral particles included in the vertex constraints.
    • Instead of constraining the vertex, the neutral particle vertex will be defined by the latest iteration of the kinematic fit, and it's momentum will be recalculated at each iteration accordingly.
      • Technically the neutral shower momentum isn't directly included in the p4 fit if it's also in a vertex fit: the shower hit position and energy (and their uncertainties) are used along with the unknown vertex instead.

Vertex Fits

  • The initial vertex guesses must be supplied to the kinematic fitter, and it is extremely important to have good initial values for them.
  • It is recommended that the following procedure be used to determine initial guesses for the vertex fits:
    • Use the DOCA routines in DAnalysisUtilities to find an initial vertex guess for a set of tracks.
    • Prior to performing a simultaneous vertex-p4 fit, first kinematically fit each vertex individually using the above initial vertex guesses. Then use the results from these kinematic fit as inputs to the vertex-p4 fit.
  • To disable linking of detached vertices (e.g. using the Λ to fit the production vertex in γp→K+Λ, Λ→pπ-), run with the command-line argument: -PKINFIT:LINKVERTICES=0

Code Documentation

Primary Classes

  • The DKinFitter C++ class is designed to be usable for any experiment, fixed-target or multi-beam, and with or without a magnetic field present. It's only dependency is ROOT, from which it uses the vector and matrix classes. The only assumption that the fitter makes is that the uncertainties in the input particle covariance matrices are in the following order:
    • Particles: px, py, pz, xx, xy, xz, t
    • Showers: E, xx, xy, xz, t
  • The DKinFitUtils C++ class is used to supplement the kinematic fitter with utility functions used to help manage the fit, and to perform some calculations.
    • It is used for creating the particles, constraints, and other objects used in the fit, and is responsible for managing their memory.
    • It is also used to provide access to information about the presence, strength, or direction of the magnetic field.
  • However, since the magnetic field can vary from experiment to experiment, the DKinFitUtils class itself knows nothing about it. Instead, it is an abstract class, and each experiment must define a class inheriting from it that defines the following member functions (including two that are useful for helping to setup constraints):
virtual bool Get_IncludeBeamlineInVertexFitFlag(void) const = 0;
virtual bool Get_IsDetachedVertex(int locPDG_PID) const = 0;
virtual bool Get_IsBFieldNearBeamline(void) const = 0;
virtual TVector3 Get_BField(const TVector3& locPosition) const = 0; //must return in units of Tesla!!
  • Thus DKinFitUtils is an abstract base class; it cannot be instantiated directly. Instead, for each experiment a class should inherit from DKinFitUtils and define these methods; for GlueX, this is the DKinFitUtils_GlueX class.

DKinFitter Data Classes

The following classes are used to hold and handle data for the kinematic fit. The memory for all of these classes is managed internally (with resource pools) by the DKinFitUtils class.

  • DKinFitParticle: One for each particle in the fit. There are several types of particles, which are differentiated by the DKinFitParticleType enum.
    • DKinFitParticleType: An enum used to distinguish different particle types. The different types of particles are: d_DetectedParticle, d_BeamParticle, d_TargetParticle, d_DecayingParticle, d_MissingParticle
  • DKinFitConstraint: This is an abstract base class that has the following derived classes:
  • DKinFitConstraint_P4: Constrains a set of initial- and final-state particles to conserve four-momentum.
  • DKinFitConstraint_Mass: Constrains the four-momentum of a decaying particle to have a fixed mass.
  • DKinFitConstraint_Vertex: Used for constraining a set of particles to a common vertex.
  • DKinFitConstraint_Spacetime: Used for constraining a set of particles to a common spacetime. This class inherits from DKinFitConstraint_Vertex.
  • DKinFitChain & DKinFitChainStep: Used for representing the reaction to be fit as a series of decay-steps in a decay chain, containing the created DKinFitParticle's.
    • This class is not necessary to perform the fit, however, it can be useful for using some of the (optional) DKinFitUtils functions to help setup the constraints for your fits.

In addition, two resource pools of TMatrixDSym objects are used to manage memory for the covariance matrices: One for the matrices for the individual particles, and one for the large matrices used for the fit parameters and their correlations.

Fit Setup - New Run

  • Required: First, create pointers to the DKinFitter & DKinFitUtils_GlueX classes in your header file:
DKinFitter dKinFitter;
DKinFitUtils_GlueX dKinFitUtils;
  • Required: Create the objects in your .cc file. This also sets the magnetic field.
dKinFitUtils = new DKinFitUtils_GlueX(locEventLoop);
dKinFitter = new DKinFitter(dKinFitUtils);
  • Optional: Set sizes for resource pools, set the debug flag, etc.

Fit Setup - New Event

  • Call DKinFitter::Reset_NewEvent() to reset the fitter (and the utilities) at the beginning of each new event. This recycles the memory (& output/input mapping) of all kinfit objects for reuse. These fit results objects are stored in DKinFitResults, so they must be valid for the entire event.

Fit Setup - New Particle Combination

  • Create DKinFitParticle objects for your particles. These can be created by calling one of the following methods of DKinUtils_GlueX:
DKinFitParticle* Make_BeamParticle(const DBeamPhoton* locBeamPhoton);
DKinFitParticle* Make_BeamParticle(const DBeamPhoton* locBeamPhoton, const DEventRFBunch* locEventRFBunch); //sets rf time for photon
DKinFitParticle* Make_TargetParticle(Particle_t locPID);
DKinFitParticle* Make_DetectedParticle(const DKinematicData* locKinematicData);
DKinFitParticle* Make_DetectedShower(const DNeutralShower* locNeutralShower, Particle_t locPID); //DO NOT call this unless the neutral is also in a vertex fit!
DKinFitParticle* Make_DecayingParticle(Particle_t locPID, const set<DKinFitParticle*>& locFromInitialState, const set<DKinFitParticle*>& locFromFinalState);
DKinFitParticle* Make_MissingParticle(Particle_t locPID);
  • Note that for decaying particles, you must pass in the pre-created DKinFitParticle objects that are used to define the momentum of the decaying particle.

Fit Setup - New Fit

  • Reset the kinematic fitter for the new fit. This clears the list of "current-fit" particles and constraints, gets new resources for the fit covariance matrices, etc.:
  • Create the constraints used for the fit. These can be created manually by calling the following functions:
DKinFitConstraint_Mass* DKinFitUtils::Make_MassConstraint(DKinFitParticle* locDecayingParticle);
DKinFitConstraint_P4* DKinFitUtils::Make_P4Constraint(const set<DKinFitParticle*>& locInitialParticles, const set<DKinFitParticle*>& locFinalParticles);
DKinFitConstraint_Vertex* DKinFitUtils::Make_VertexConstraint(const set<DKinFitParticle*>& locFullConstrainParticles, const set<DKinFitParticle*>& locNoConstrainParticles, TVector3 locVertexGuess = TVector3());
DKinFitConstraint_Spacetime* DKinFitUtils::Make_SpacetimeConstraint(const set<DKinFitParticle*>& locFullConstrainParticles, const set<DKinFitParticle*>& locOnlyConstrainTimeParticles, 
	const set<DKinFitParticle*>& locNoConstrainParticles, TLorentzVector locSpacetimeGuess = TLorentzVector());
  • Or, the constraints can be created automatically by calling the following function:
set<DKinFitConstraint*> DKinFitUtils_GlueX::Create_Constraints(const DParticleCombo* locParticleCombo, const DKinFitChain* locKinFitChain, DKinFitType locKinFitType, deque<DKinFitConstraint_Vertex*>& locSortedVertexConstraints);
  • Next, set initial guesses for any common vertices, common times, and the p3 of a missing or open-ended decaying particle, if present. Do this by calling:
void DKinFitConstraint_P4::Set_InitP3Guess(const TVector3& locInitP3Guess);
void DKinFitConstraint_Vertex::Set_InitVertexGuess(const TVector3& locInitVertexGuess);
void DKinFitConstraint_Spacetime::Set_InitVertexGuess(const TVector3& locInitVertexGuess);
void DKinFitConstraint_Spacetime::Set_InitTimeGuess(double locInitTimeGuess);
  • Or, the vertex and time guesses can be calculated and set automatically via (where locSortedVertexConstraints should be in the order that you want the guesses to be obtained):
void DKinFitUtils_GlueX::Set_SpacetimeGuesses(const deque<DKinFitConstraint_Vertex*>& locSortedVertexConstraints, bool locIsP4FitFlag);
  • Of the constraints created above, set the constraints you would like to apply for the current fit, by calling either:
void DKinFitter::Add_Constraint(DKinFitConstraint* locKinFitConstraint);
void DKinFitter::Add_Constraints(const set<DKinFitConstraint*>& locKinFitConstraints);

Performing the Kinematic Fit

  • Perform the kinematic fit by calling:
bool DKinFitter::Fit_Reaction(void);
  • This function returns true if the fit was successful, and false if it failed. The fit can fail for one of several reasons:
    • The input constraints were not setup correctly. For example, if a detected particle is constrained to multiple vertices, or is in multiple four-momentum constraints, or if not enough particles are available to constrain a vertex (2 are needed), etc.
    • The fit fails to converge to a minimum χ2.
    • If during the kinematic fit, one of matrices is not invertible (e.g. the determinant is zero).

Retrieving Fit Results

  • The following methods can be called to retrieve information on the quality of the most-recently-completed fit.
double DKinFitter::Get_ChiSq(void) const;
double DKinFitter::Get_ConfidenceLevel(void) const;
unsigned int DKinFitter::Get_NDF(void) const;
  • Many different fits with different constraints can be placed on the same set of particles in a given event or particle combination. This includes re-using the DKinFitParticle and DKinFitConstraint objects that were created earlier. However, each fit needs to update the "active" DKinFitParticle and DKinFitConstraint objects as each iteration improves the particle/vertex/etc. properties. Thus, in order to avoid over-writing the user-created DKinFitParticle and DKinFitConstraint objects, clones of these objects are made when DKinFitter::Fit_Reaction() is called. The cloned objects are used in the fit and contain the fit results. Thus, to obtain the fit results for a given particle / constraint, the cloned, output object needs to be retrieved via:
set<DKinFitConstraint*> DKinFitter::Get_KinFitConstraints(void) const;
set<DKinFitParticle*> Get_KinFitParticles(void) const;
  • The output DKinFitParticle's can be mapped to the input DKinFitParticle's via:
DKinFitParticle* DKinFitUtils::Get_InputKinFitParticle(DKinFitParticle* locKinFitParticle) const;
  • And, the input DKinFitParticle's can be mapped to the original JObject's via:
const JObject* DKinFitUtils_GlueX::Get_SourceJObject(DKinFitParticle* locInputKinFitParticle) const;
  • For a more detailed look at how the fit parameters have changed with respect to their input, measured values, the pulls can be retrieved:
//1st key is particle, 2nd key is param type (defined in DKinFitParticle.h)
void DKinFitter::Get_Pulls(map<DKinFitParticle*, map<DKinFitPullType, double> >& locPulls) const;
  • The full covariance matrices for the fit parameters can be retrieved via:
const TMatrixDSym* Get_VEta(void); //covariance matrix of the measured parameters used in the fit
const TMatrixDSym* Get_VXi(void); //covariance matrix of the unknown parameters used in the fit
const TMatrixDSym* Get_V(void); //combined covariance matrix of the measured and unknown parameters used in the fit.

Duplicate Fit Results

The DKinFitResults_factory class will auto-detect whether the requested fit for the given particle combination is identical to a fit that has been previously performed that event. For example, fits can be identical when:

  • Two DParticleCombo's from the same DReaction differ only by a final-state photon, and a vertex-fit is requested (the photon is not used in the vertex fit).
  • Multiple DReaction's are used for the same DReactionStep contents. For example, each DReaction may have different actions but the same particles.
  • Two DReaction's are identical except for the inclusion of an ω or φ meson. For example: γ p → π+ π- π0 p vs. γ p → ω p; ω → π+ π- π0. The ω and φ mesons don't have a fixed mass so they are not constrained, nor do they have a detached vertex. Therefore, the kinematic fit results are independent of their inclusion.
  • Etc.

In these cases, the DKinFitResults_factory will simply assign the fit results from the previously-fit DParticleCombo to the new DParticleCombo.

Auto-detection of these cases is performed by creating the DKinFitConstraint's for the new DParticleCombo and comparing them to the previously-fit DParticleCombo.

Manual Kinematic Fit Example

  • To see an example on how to manually perform a kinematic fit (and cut on it's confidence level) inside of a DAnalysisAction, see DCustomAction_KinFitExample.* at Link
    • This assumes that the DReaction is setup as γ, p → K+, Λ; Λ → p, π-


#ifndef _DKinFitAction_Example_
#define _DKinFitAction_Example_
#include "JANA/JEventLoop.h"
#include "DANA/DApplication.h"
#include "ANALYSIS/DAnalysisAction.h"
#include "ANALYSIS/DKinFitter_GlueX.h"
#include "ANALYSIS/DKinFitParticle.h"
#include "ANALYSIS/DParticleCombo.h"
using namespace std;
using namespace jana;
class DKinFitAction_Example : public DAnalysisAction
    DKinFitAction_Example(const DReaction* locReaction, string locActionUniqueString = "") : 
    DAnalysisAction(locReaction, "KinFit_Example", false, locActionUniqueString){}
    bool Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo);
    void Initialize(JEventLoop* locEventLoop);
    DKinFitter_GlueX dKinFitter;
    //define any histograms here
#endif // _DKinFitAction_Example_

#include "DKinFitAction_Example.h"
void DKinFitAction_Example::Initialize(JEventLoop* locEventLoop)
  DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
  const DMagneticFieldMap* locMagneticFieldMap = locApplication->GetBfield(locEventLoop->GetJEvent().GetRunNumber());
  double locTargetZCenter = 65.0;
  DGeometry* locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
  //Only set magnetic field if non-zero!
  double locBx, locBy, locBz;
  locMagneticFieldMap->GetField(0.0, 0.0, locTargetZCenter, locBx, locBy, locBz);
  TVector3 locBField(locBx, locBy, locBz);
  if(locBField.Mag() > 0.0)
  //create any histograms here
bool DKinFitAction_Example::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
  //THIS ASSUMES DParticleCombo (read: DReaction) is setup as:
    //g, p -> K+, Lambda
    //Lambda -> p, pi-
  dKinFitter.Reset_NewEvent(); //need to call prior to use in each event (cleans up memory allocated from last event)
  const DBeamPhoton* locBeamPhoton = static_cast<const DBeamPhoton*>(locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle_Measured());
  const DKinFitParticle* locKinFitParticle_BeamPhoton = dKinFitter.Make_BeamParticle(locBeamPhoton);
  const DKinFitParticle* locKinFitParticle_Target = dKinFitter.Make_TargetParticle(Proton);
  const DKinFitParticle* locKinFitParticle_Lambda = dKinFitter.Make_DecayingParticle(Lambda);
  const DChargedTrackHypothesis* locChargedTrackHypothesis_KPlus = static_cast<const DChargedTrackHypothesis*>(locParticleCombo->Get_ParticleComboStep(0)->Get_FinalParticle_Measured(0));
  const DKinFitParticle* locKinFitParticle_KPlus = dKinFitter.Make_DetectedParticle(locChargedTrackHypothesis_KPlus);
  const DChargedTrackHypothesis* locChargedTrackHypothesis_Proton = static_cast<const DChargedTrackHypothesis*>(locParticleCombo->Get_ParticleComboStep(1)->Get_FinalParticle_Measured(0));
  const DKinFitParticle* locKinFitParticle_Proton = dKinFitter.Make_DetectedParticle(locChargedTrackHypothesis_Proton);
  const DChargedTrackHypothesis* locChargedTrackHypothesis_PiMinus = static_cast<const DChargedTrackHypothesis*>(locParticleCombo->Get_ParticleComboStep(1)->Get_FinalParticle_Measured(1));
  const DKinFitParticle* locKinFitParticle_PiMinus = dKinFitter.Make_DetectedParticle(locChargedTrackHypothesis_PiMinus);
  dKinFitter.Reset_NewFit(); //disregards the constraints from the previous kinematic fit
  deque<const DKinFitParticle*> locInitialKinFitParticles, locFinalKinFitParticles;
  // first p4 constraint:
  // make constraint
  DKinFitConstraint_P4* locP4Constraint = dKinFitter.Make_P4Constraint(locInitialKinFitParticles, locFinalKinFitParticles, false); //false: don't constrain mass of initial particle (beam photon)
  // set constraint
  // second p4 constraint:
  // make constraint
  locP4Constraint = dKinFitter.Make_P4Constraint(locInitialKinFitParticles, locFinalKinFitParticles, true); //true: constrain invariant mass of initial particle (Lambda)
  // set constraint
  // vertex constraint
  TVector3 locVertexGuess(locKinFitParticle_Proton->Get_Position()); //try to get a better vertex guess than this!
  locInitialKinFitParticles.push_back(locKinFitParticle_Lambda); //here the lambda isn't constraining the vertex (since it's only in one vertex fit); instead it's vertex will be defined by the fit result
  // make constraint
  DKinFitConstraint_Vertex* locVertexConstraint = dKinFitter.Make_VertexConstraint(locInitialKinFitParticles, locFinalKinFitParticles, locVertexGuess);
  // set constraint
  double locConfidenceLevel = dKinFitter.Get_ConfidenceLevel();
  double locChiSq = dKinFitter.Get_ChiSq();
  unsigned int locNDF = dKinFitter.Get_NDF();
  const TMatrixDSym* locMeasuredParametersKinFitCovarianceMatrix = dKinFitter.Get_VEta();
  const TMatrixDSym* locUnknownParametersCovarianceMatrix = dKinFitter.Get_VXi();
  const DKinFitParticle* locKinFitParticle_Proton_Output = dKinFitter.Get_OutputKinFitParticle(locKinFitParticle_Proton);
  const DKinFitParticle* locKinFitParticle_PiMinus_Output = dKinFitter.Get_OutputKinFitParticle(locKinFitParticle_PiMinus);
  const DKinFitParticle* locKinFitParticle_KPlus_Output = dKinFitter.Get_OutputKinFitParticle(locKinFitParticle_KPlus);
  const DKinFitParticle* locKinFitParticle_Lambda_Output = dKinFitter.Get_OutputKinFitParticle(locKinFitParticle_Lambda);
  const DKinFitParticle* locKinFitParticle_BeamPhoton_Output = dKinFitter.Get_OutputKinFitParticle(locKinFitParticle_BeamPhoton);
  TVector3 locKinFitVertex = locKinFitParticle_PiMinus_Output->Get_CommonVertex(); //could also be grabbed from the proton or the lambda
  map<const DKinematicData*, map<DKinFitPullType, double> > locPulls; //map from input particle data to pull map //pull map: map from pull type to value (DKinFitPullType defined in DKinFitParticle.h)
  double locProtonPxPull = (locPulls[locChargedTrackHypothesis_Proton])[d_PxPull];
  //fill histograms, etc. here
  return (locConfidenceLevel >= 0.01); //cut if < 1% confidence level
  • Somewhere in the DReaction_factory in your plugin:
#include "DKinFitAction_Example.h"
  locReaction->Add_AnalysisAction(new DKinFitAction_Example(locReaction));