Difference between revisions of "Mattione GlueX Kinematic Fitting"

From GlueXWiki
Jump to: navigation, search
(DKinFitResults_factory.cc)
(Duplicate Fit Results)
 
(90 intermediate revisions by the same user not shown)
Line 1: Line 1:
 
== Quick Start ==
 
== Quick Start ==
* The kinematic fitter (<span style="color:#0000FF">DKinFitter_GlueX</span>) is fully integrated into the analysis libraries (interfaced to <span style="color:#0000FF">DKinFitResults_factory</span>), so to kinematically fit a <span style="color:#0000FF">DReaction</span> you only have to specify the desired kinematic fit type in the <span style="color:#0000FF">DReaction</span>.
+
* The kinematic fitter (<span style="color:#0000FF">DKinFitter</span>) is in the KINFITTER library, and is fully integrated into the analysis library (via <span style="color:#0000FF">DKinFitUtils_GlueX</span> and the <span style="color:#0000FF">DKinFitResults_factory</span>). So, to kinematically fit a <span style="color:#0000FF">DReaction</span>, you only have to specify the desired kinematic fit type in the <span style="color:#0000FF">DReaction</span>.  
** See [[Mattione_GlueX_Analysis_Factories | 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.
 
** 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.
 
* Details on the kinematic fitting equations and constraints are in GlueX-doc-2112.
  
== Support ==
+
== Not Yet Supported ==
* Note that spacetime constraints and fits including final-state massive neutral particles are not yet supported.
+
* Spacetime constraints
  
 
== Constraints and Unknowns ==
 
== Constraints and Unknowns ==
 
* Up to one particle may be missing: 3 unknowns (px, py, pz)
 
* Up to one particle may be missing: 3 unknowns (px, py, pz)
 
* Four-momentum constraint: 4 constraints (Initial p4 - Final p4 = '''0''')
 
* Four-momentum constraint: 4 constraints (Initial p4 - Final p4 = '''0''')
* Invariant mass constraints: 1 constraint (E<sup>2</sup> - '''p'''<sup>2</sup> - m<sup>2</sup> = 0)
+
* Mass constraints: 1 constraint per particle (E<sup>2</sup> - '''p'''<sup>2</sup> - m<sup>2</sup> = 0)
** Note that by default, the ANALYSIS library does not constrain the masses of the <span style="color:red">&omega;</span> or <span style="color:red">&phi;</span> resonances, even though their widths are very narrowThey can be constrained though if you setup the kinematic fit manually.
+
** Note that the ANALYSIS library does not constrain the masses of the <span style="color:red">&omega;</span> or <span style="color:red">&phi;</span> mesons since their widths are non-negligible.   
* Vertex constraints: 3 unknowns per vertex (x, y, z), 2 constraints per track
+
** 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 (&Delta;'''x''') = '''0''': The path (&Delta;'''x''') and momentum vectors must be parallel
 
** E.g. for no magnetic field, '''p''' x (&Delta;'''x''') = '''0''': The path (&Delta;'''x''') and momentum vectors must be parallel
 
*** &Delta;'''x''' = '''x'''<sub>Measured</sub> - '''x'''<sub>Unknown</sub>
 
*** &Delta;'''x''' = '''x'''<sub>Measured</sub> - '''x'''<sub>Unknown</sub>
* Spacetime constraints ('''NOT YET SUPPORTED'''): 4 unknowns per spacetime point (x, y, z, t), 3 constraints per track (vertex + time)
+
* 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)
  
== DKinFitResuls_factory ==
+
== Notes for those manually performing their own fits ==
* This factory sets up and performs the kinematic fit that is indicated by the <span style="color:#0000FF">DReaction</span>::<span style="color:#008000">dKinFitType</span> variable.  This variable is an enum type called <span style="color:#0000FF">DKinFitType</span>, which is defined in libraries/ANALYSIS/DKinFitResults.h
+
* This factory automatically determines what p4, vertex, and spacetime (not yet supported!) constraints to setup such that each the data in each <span style="color:#0000FF">DParticleCombo</span> are maximally constrained.
+
** The factory will automatically create additional vertex constraints for decaying particles with long lifetimes, as long as the <span style="color:#008000">IsDetachedVertex</span>() function in libraries/include/particleType.h is kept updated.
+
** The factory will also automatically exclude resonances from the p4 constraints of kinematic fits, as long as the <span style="color:#008000">IsFixedMass</span>() function in libraries/include/particleType.h is kept updated.
+
** If a missing or decaying particle of type <span style="color:red">Unknown</span> is present in the <span style="color:#0000FF">DReaction</span>, the fit will not constrain the overall p4 of the reaction, but it will still apply the other constraints.
+
* To exclude additional, non-resonant particles from the p4 constraints of a kinematic fit, call the <span style="color:#0000FF">DReactionStep</span>::<span style="color:#008000">Set_ApplyKinFitMassConstraintOnInitialParticleFlag</span>(<span style="color:#0000FF">bool</span>) function with an argument of <span style="color:red">true</span> on the step where the particle is in the initial state.
+
  
== Important Notes ==
+
=== DKinFitUtils_GlueX ===
 +
* These utilities are used to setup the kinematic fit that is indicated by the <span style="color:#0000FF">DReaction</span>::<span style="color:#008000">dKinFitType</span> variable.  This variable is an enum type called <span style="color:#0000FF">DKinFitType</span>, which is defined in libraries/ANALYSIS/DReaction.h
 +
* These utilities automatically determines what p4, vertex, and spacetime (not yet supported!) constraints to setup such that each the data in each <span style="color:#0000FF">DParticleCombo</span> are maximally constrained.
 +
** It will also automatically create additional vertex constraints for decaying particles with long lifetimes, as long as the <span style="color:#008000">IsDetachedVertex</span>() function in libraries/include/particleType.h is kept updated.
 +
** If an inclusive channel, or if a decaying particle of type <span style="color:red">Unknown</span> is present in the <span style="color:#0000FF">DReaction</span>, 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 <span style="color:#0000FF">DReactionStep</span>::<span style="color:#008000">Set_KinFitConstrainInitMassFlag</span>(<span style="color:#0000FF">bool</span>) function with an argument of <span style="color:red">true</span> on the step where the decaying particle is in the initial state.
  
 
=== Detected Neutral Particles ===
 
=== Detected Neutral Particles ===
* The <span style="color:#0000FF">DNeutralParticleHypothesis</span> objects from the default factory sets the particle momentum assuming that the particle came from the center of the target.
+
* The "measured" <span style="color:#0000FF">DNeutralParticleHypothesis</span> objects in the combos sets the particle momentum assuming that the particle came from the position roughly defined during comboing (see DSourceComboVertexer).  
 
* 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.   
 
* 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.   
 
** 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.
+
*** 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.
 
+
* Choosing to use the neutral shower instead of the default <span style="color:#0000FF">DNeutralParticleHypothesis</span> object is toggled by the boolean flag when calling the <span style="color:#0000FF">DKinFitter_GlueX</span>::<span style="color:#008000">Make_DetectedParticle</span>(const <span style="color:#0000FF">DNeutralParticleHypothesis</span>*, <span style="color:#0000FF">bool</span>) function.
+
  
 
=== Vertex Fits ===
 
=== 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.   
 
* The initial vertex guesses must be supplied to the kinematic fitter, and it is extremely important to have good initial values for them.   
** This is especially true for fits that also contain p4 constraints, to ensure that the fit converges properly.
 
 
* It is recommended that the following procedure be used to determine initial guesses for the vertex fits:
 
* It is recommended that the following procedure be used to determine initial guesses for the vertex fits:
 
** Use the DOCA routines in <span style="color:#0000FF">DAnalysisUtilities</span> to find an initial vertex guess for a set of tracks.
 
** Use the DOCA routines in <span style="color:#0000FF">DAnalysisUtilities</span> 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.
 
** 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 &Lambda; to fit the production vertex in &gamma;p&rarr;K<sup>+</sup>&Lambda;, &Lambda;&rarr;p&pi;<sub>-</sub>), run with the command-line argument: -PKINFIT:LINKVERTICES=0
 
** This is useful for when you do not want to use a mass constraint on the decaying particle.
 
  
== Manual Kinematic Fit Example ==
+
== Code Documentation ==
* This example shows how to manually perform a kinematic fit (and cut on it's confidence level) inside of a DAnalysisAction.
+
** This assumes that the DReaction is setup as &gamma;, p &rarr; K<sup>+</sup>, &Lambda;; &Lambda; &rarr; p, &pi;<sup>-</sup>
+
=== DKinFitAction_Sample.h ===
+
<syntaxhighlight>
+
#ifndef _DKinFitAction_Sample_
+
#define _DKinFitAction_Sample_
+
  
#include "JANA/JEventLoop.h"
+
=== Primary Classes ===
#include "DANA/DApplication.h"
+
  
#include "ANALYSIS/DAnalysisAction.h"
+
* The <span style="color:#0000FF">DKinFitter</span> 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 primary dependency is ROOT, from which it uses the vector and matrix classes. It also depends on  <span style="color:#0000FF">DResettable</span> and <span style="color:#0000FF">DResourcePool</span> in sim-recon/src/libraries/include/ to help with memory management. The only assumption that the fitter makes is that the uncertainties in the input particle covariance matrices are in the following order:
#include "ANALYSIS/DKinFitter_GlueX.h"
+
** Particles: p<sub>x</sub>, p<sub>y</sub>, p<sub>z</sub>, x<sub>x</sub>, x<sub>y</sub>, x<sub>z</sub>, t
#include "ANALYSIS/DKinFitParticle.h"
+
** Showers: E, x<sub>x</sub>, x<sub>y</sub>, x<sub>z</sub>, t
#include "ANALYSIS/DParticleCombo.h"
+
  
using namespace std;
+
* The <span style="color:#0000FF">DKinFitUtils</span> C++ class is used to supplement the kinematic fitter with utility functions used to help manage the fit, and to perform some calculations. 
using namespace jana;
+
** 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.
  
class DKinFitAction_Sample : public DAnalysisAction
+
* However, since the magnetic field can vary from experiment to experiment, the <span style="color:#0000FF">DKinFitUtils</span> 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):
{
+
  public:
+
    DKinFitAction_Sample(const DReaction* locReaction, string locActionUniqueString = "") :  
+
    DAnalysisAction(locReaction, "KinFit_Sample", false, locActionUniqueString){}
+
  
  private:
+
<syntaxhighlight>
    bool Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo);
+
//FOR SETUP:
    void Initialize(JEventLoop* locEventLoop);
+
virtual bool Get_IncludeBeamlineInVertexFitFlag(void) const = 0;
 
+
    DKinFitter_GlueX dKinFitter;
+
 
+
    //define any histograms here
+
};
+
  
#endif // _DKinFitAction_Sample_
+
//FOR FITS:
 +
virtual bool Get_IsBFieldNearBeamline(void) const = 0;
 +
virtual TVector3 Get_BField(const TVector3& locPosition) const = 0; //must return in units of Tesla!!
 
</syntaxhighlight>
 
</syntaxhighlight>
  
=== DKinFitAction_Sample.cc ===
+
* Thus <span style="color:#0000FF">DKinFitUtils</span> is an abstract base class; it cannot be instantiated directly. Instead, for each experiment a class should inherit from <span style="color:#0000FF">DKinFitUtils</span> and define these methods; for GlueX, this is the <span style="color:#0000FF">DKinFitUtils_GlueX</span> class.
<syntaxhighlight>
+
#include "DKinFitAction_Sample.h"
+
  
void DKinFitAction_Sample::Initialize(JEventLoop* locEventLoop)
+
=== DKinFitter Data Classes ===
{
+
  DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
+
  
  double locTargetZCenter = 65.0;
+
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 <span style="color:#0000FF">DKinFitUtils</span> class.
  DGeometry* locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
+
  locGeometry->GetTargetZ(locTargetZCenter);
+
  
  //Only set magnetic field if non-zero!
+
* <span style="color:#0000FF">DKinFitParticle</span>: One for each particle in the fit. There are several types of particles, which are differentiated by the <span style="color:#0000FF">DKinFitParticleType</span> enum.
  double locBx, locBy, locBz;
+
** <span style="color:#0000FF">DKinFitParticleType</span>: An enum used to distinguish different particle types. The different types of particles are: <span style="color:red">d_DetectedParticle</span>, <span style="color:red">d_BeamParticle</span>, <span style="color:red">d_TargetParticle</span>, <span style="color:red">d_DecayingParticle</span>, <span style="color:red">d_MissingParticle</span>
  locMagneticFieldMap->GetField(0.0, 0.0, locTargetZCenter, locBx, locBy, locBz);
+
* <span style="color:#0000FF">DKinFitConstraint</span>: This is an abstract base class that has the following derived classes:
  TVector3 locBField(locBx, locBy, locBz);
+
** <span style="color:#0000FF">DKinFitConstraint_P4</span>: Constrains a set of initial- and final-state particles to conserve four-momentum.  
  if(locBField.Mag() > 0.0)
+
** <span style="color:#0000FF">DKinFitConstraint_Mass</span>: Constrains the four-momentum of a decaying particle to have a fixed mass.  
    dKinFitter.Set_BField(locMagneticFieldMap);
+
** <span style="color:#0000FF">DKinFitConstraint_Vertex</span>: Used for constraining a set of particles to a common vertex.  
 +
** <span style="color:#0000FF">DKinFitConstraint_Spacetime</span>: Used for constraining a set of particles to a common spacetime. This class inherits from <span style="color:#0000FF">DKinFitConstraint_Vertex</span>.
 +
* <span style="color:#0000FF">DKinFitChain</span> & <span style="color:#0000FF">DKinFitChainStep</span>: Used for representing the reaction to be fit as a series of decay-steps in a decay chain, containing the created <span style="color:#0000FF">DKinFitParticle</span>'s.
 +
** These classes are not necessary to perform the fit, however, they can be useful for using some of the (optional) <span style="color:#0000FF">DKinFitUtils</span> functions to help setup the constraints for your fits.
  
  //create any histograms here
+
In addition, resource pools of <span style="color:#0000FF">TMatrixFSym</span> objects are used to manage memory for the covariance matrices.
}
+
  
bool DKinFitAction_Sample::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
+
=== Fit Setup - New Run ===
{
+
* Required: First, create pointers to the <span style="color:#0000FF">DKinFitter</span> & <span style="color:#0000FF">DKinFitUtils_GlueX</span> classes in your header file:
  //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)
+
<syntaxhighlight>
 +
#include "KINFITTER/DKinFitter.h"
 +
#include "ANALYSIS/DKinFitUtils_GlueX.h"
  
  //CREATE DKINFITPARTICLE OBJECTS FOR EACH PARTICLE
+
//in class definition:
  const DBeamPhoton* locBeamPhoton = static_cast<const DBeamPhoton*>(locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle_Measured());
+
DKinFitter* dKinFitter;
  const DKinFitParticle* locKinFitParticle_BeamPhoton = dKinFitter.Make_BeamParticle(locBeamPhoton);
+
DKinFitUtils_GlueX* dKinFitUtils;
 +
</syntaxhighlight>
  
  const DKinFitParticle* locKinFitParticle_Target = dKinFitter.Make_TargetParticle(Proton);
+
* Required: Create the objects in your .cc file.  This also sets the magnetic field.  
  const DKinFitParticle* locKinFitParticle_Lambda = dKinFitter.Make_DecayingParticle(Lambda);
+
  
  const DChargedTrackHypothesis* locChargedTrackHypothesis_KPlus = static_cast<const DChargedTrackHypothesis*>(locParticleCombo->Get_ParticleComboStep(0)->Get_FinalParticle_Measured(0));
+
<syntaxhighlight>
  const DKinFitParticle* locKinFitParticle_KPlus = dKinFitter.Make_DetectedParticle(locChargedTrackHypothesis_KPlus);
+
dKinFitUtils = new DKinFitUtils_GlueX(locEventLoop);
 +
dKinFitter = new DKinFitter(dKinFitUtils);
 +
</syntaxhighlight>
  
  const DChargedTrackHypothesis* locChargedTrackHypothesis_Proton = static_cast<const DChargedTrackHypothesis*>(locParticleCombo->Get_ParticleComboStep(1)->Get_FinalParticle_Measured(0));
+
* Optional: Set sizes for resource pools, set the debug flag, etc.
  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));
+
=== Fit Setup - New Event ===
  const DKinFitParticle* locKinFitParticle_PiMinus = dKinFitter.Make_DetectedParticle(locChargedTrackHypothesis_PiMinus);
+
* At the beginning of each event, the kinematic fitter (and the utilities) need to be reset This recycles the memory (& output/input object mapping) of all objects for reuse. After the fit, these fit results objects are stored in <span style="color:#0000FF">DKinFitResults</span>, so they must stay valid for the entire event. To reset everything, call:
  
  // SETUP THE CONSTRAINTS
 
  dKinFitter.Reset_NewFit(); //disregards the constraints from the previous kinematic fit
 
  deque<const DKinFitParticle*> locInitialKinFitParticles, locFinalKinFitParticles;
 
 
  // first p4 constraint:
 
  locInitialKinFitParticles.push_back(locKinFitParticle_BeamPhoton);
 
  locInitialKinFitParticles.push_back(locKinFitParticle_Target);
 
  locFinalKinFitParticles.push_back(locKinFitParticle_KPlus);
 
  locFinalKinFitParticles.push_back(locKinFitParticle_Lambda);
 
  // make constraint
 
  DKinFitConstraint_P4* locP4Constraint = dKinFitter.Make_P4Constraint(locInitialKinFitParticles, locFinalKinFitParticles, false); //false: don't constrain mass of initial particle (beam photon)
 
  // set constraint
 
  dKinFitter.Set_Constraint(locP4Constraint);
 
 
  // second p4 constraint:
 
  locInitialKinFitParticles.clear();
 
  locFinalKinFitParticles.clear();
 
  locInitialKinFitParticles.push_back(locKinFitParticle_Lambda);
 
  locFinalKinFitParticles.push_back(locKinFitParticle_PiMinus);
 
  locFinalKinFitParticles.push_back(locKinFitParticle_Proton);
 
  // make constraint
 
  locP4Constraint = dKinFitter.Make_P4Constraint(locInitialKinFitParticles, locFinalKinFitParticles, true); //true: constrain invariant mass of initial particle (Lambda)
 
  // set constraint
 
  dKinFitter.Set_Constraint(locP4Constraint);
 
 
  // vertex constraint
 
  TVector3 locVertexGuess(locKinFitParticle_Proton->Get_Position()); //try to get a better vertex guess than this!
 
  locInitialKinFitParticles.clear();
 
  locFinalKinFitParticles.clear();
 
  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
 
  locFinalKinFitParticles.push_back(locKinFitParticle_PiMinus);
 
  locFinalKinFitParticles.push_back(locKinFitParticle_Proton);
 
  // make constraint
 
  DKinFitConstraint_Vertex* locVertexConstraint = dKinFitter.Make_VertexConstraint(locInitialKinFitParticles, locFinalKinFitParticles, locVertexGuess);
 
  // set constraint
 
  dKinFitter.Set_Constraint(locVertexConstraint);
 
 
  // PERFORM THE KINEMATIC FIT
 
  dKinFitter.Fit_Reaction();
 
 
  // GET THE FIT RESULTS
 
  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)
 
  dKinFitter.Get_Pulls(locPulls);
 
  double locProtonPxPull = (locPulls[locChargedTrackHypothesis_Proton])[d_PxPull];
 
 
  //fill histograms, etc. here
 
  return (locConfidenceLevel >= 0.01); //cut if < 1% confidence level
 
}
 
</syntaxhighlight>
 
 
* Somewhere in the DReaction_factory in your plugin:
 
 
<syntaxhighlight>
 
<syntaxhighlight>
#include "DKinFitAction_Sample.h"
+
void DKinFitUtils_GlueX::Reset_NewEvent(void);
//...
+
void DKinFitter::Reset_NewEvent(void);
  locReaction->Add_AnalysisAction(new DKinFitAction_Sample(locReaction));
+
 
</syntaxhighlight>
 
</syntaxhighlight>
  
== Code Documentation - For End Users ==
+
=== Fit Setup - New Particle Combination ===
  
=== DKinFitter_GlueX and DKinFitter ===
+
* Create <span style="color:#0000FF">DKinFitParticle</span> objects for your particles. These can be created by calling one of the following methods of <span style="color:#0000FF">DKinUtils_GlueX</span>:
  
* The <span style="color:#0000FF">DKinFitter</span> C++ class is designed to be usable for any fixed-target experiment, 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:
+
<syntaxhighlight>
** Particles: p<sub>x</sub>, p<sub>y</sub>, p<sub>z</sub>, x<sub>x</sub>, x<sub>y</sub>, x<sub>z</sub>, t
+
shared_ptr<DKinFitParticle> Make_BeamParticle(const DBeamPhoton* locBeamPhoton);
** Showers: E, x<sub>x</sub>, x<sub>y</sub>, x<sub>z</sub>, t
+
shared_ptr<DKinFitParticle> Make_BeamParticle(const DBeamPhoton* locBeamPhoton, const DEventRFBunch* locEventRFBunch); //sets rf time for photon
 +
shared_ptr<DKinFitParticle> Make_TargetParticle(Particle_t locPID);
  
* However, the <span style="color:#0000FF">DKinFitter</span> class knows nothing about the presence, strength, or direction of the magnetic field.  The following <span style="color:#0000FF">DKinFitter</span> member functions are used during the fit to access the magnetic field properties, but they have been declared as pure-virtual:
+
shared_ptr<DKinFitParticle> Make_DetectedParticle(const DKinematicData* locKinematicData);
 +
shared_ptr<DKinFitParticle> Make_DetectedShower(const DNeutralShower* locNeutralShower, Particle_t locPID); //DO NOT call this unless the neutral is also in a vertex fit!
  
<syntaxhighlight>
+
shared_ptr<DKinFitParticle> Make_DecayingParticle(Particle_t locPID, const set<shared_ptr<DKinFitParticle>>& locFromInitialState, const set<shared_ptr<DKinFitParticle>>& locFromFinalState);
virtual bool DKinFitter::Get_IsBFieldNearBeamline(void) const = 0;
+
shared_ptr<DKinFitParticle> Make_MissingParticle(Particle_t locPID);
virtual TVector3 DKinFitter::Get_BField(const TVector3& locPosition) const = 0;
+
 
</syntaxhighlight>
 
</syntaxhighlight>
  
* Thus <span style="color:#0000FF">DKinFitter</span> is an abstract base class; it cannot be instantiated directly. Instead, for each experiment a class should inherit from <span style="color:#0000FF">DKinFitter</span> and define these methods; for GlueX, this is the <span style="color:#0000FF">DKinFitter_GlueX</span> class.  
+
* Note that for decaying particles, you must pass in the pre-created <span style="color:#0000FF">DKinFitParticle</span> objects that are used to define the momentum of the decaying particle.
  
=== DKinFitter Data Classes ===
+
=== Fit Setup - New Fit ===
  
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 <span style="color:#0000FF">DKinFitter</span> class. 
+
* Create the constraints used for the fit. These can be created manually by calling the following functions on DKinFitUtils_GlueX:  
 
+
* <span style="color:#0000FF">DKinFitParticle</span>: One for each particle in the fit. There are several types of particles, which are differentiated by the <span style="color:#0000FF">DKinFitParticleType</span> enum.
+
** <span style="color:#0000FF">DKinFitParticleType</span>: An enum used to distinguish different particle types. The different types of particles are: <span style="color:red">d_DetectedParticle</span>, <span style="color:red">d_BeamParticle</span>, <span style="color:red">d_TargetParticle</span>, <span style="color:red">d_DecayingParticle</span>, <span style="color:red">d_MissingParticle</span>
+
* <span style="color:#0000FF">DKinFitConstraint</span>: This is an abstract base class that has the following derived classes:
+
* <span style="color:#0000FF">DKinFitConstraint_P4</span>: Either constrains a set of initial- and final-state particles to conserve four-momentum, or constrains the final-state particles to have the invariant mass of the initial particle.
+
* <span style="color:#0000FF">DKinFitConstraint_VertexBase</span>: An abstract base class for the vertex and spacetime constraints. Particles at the vertex are contained in one of two different member variables:
+
** <span style="color:#008000">dFullConstrainParticles</span>: Particles that are used to constrain the vertex/spacetime position (e.g. charged, decaying, or beam).
+
** <span style="color:#008000">dNoConstrainParticles</span>: Particles whose vertex/spacetime are DEFINED BY this constraint, rather than used to constrain it (e.g. missing particles, decaying particles, neutral showers).
+
* <span style="color:#0000FF">DKinFitConstraint_Vertex</span>: Used for constraining a set of particles to a common vertex.
+
* <span style="color:#0000FF">DKinFitConstraint_Spacetime</span>: Used for constraining a set of particles to a common spacetime. In addition to the particles defined in <span style="color:#0000FF">DKinFitConstraint_VertexBase</span>, it contains a new class of particles:
+
** <span style="color:#008000">dOnlyConstrainTimeParticles</span>: Particles who are used to constrain the time, but not the vertex (e.g. neutral showers).
+
 
+
In addition, two resource pools of <span style="color:#0000FF">TMatrixDSym</span> 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 an instance of the <span style="color:#0000FF">DKinFitter_GlueX</span> class in your header file:  
+
  
 
<syntaxhighlight>
 
<syntaxhighlight>
DKinFitter_GlueX dKinFitter;
+
shared_ptr<DKinFitConstraint_Mass> Make_MassConstraint(const shared_ptr<DKinFitParticle>& locDecayingParticle);
 +
shared_ptr<DKinFitConstraint_P4> Make_P4Constraint(const set<shared_ptr<DKinFitParticle>>& locInitialParticles, const set<shared_ptr<DKinFitParticle>>& locFinalParticles);
 +
shared_ptr<DKinFitConstraint_Vertex> Make_VertexConstraint(const set<shared_ptr<DKinFitParticle>>& locFullConstrainParticles, const set<shared_ptr<DKinFitParticle>>& locNoConstrainParticles, TVector3 locVertexGuess = TVector3());
 +
shared_ptr<DKinFitConstraint_Spacetime> Make_SpacetimeConstraint(const set<shared_ptr<DKinFitParticle>>& locFullConstrainParticles, const set<shared_ptr<DKinFitParticle>>& locOnlyConstrainTimeParticles,
 +
const set<shared_ptr<DKinFitParticle>>& locNoConstrainParticles, TLorentzVector locSpacetimeGuess = TLorentzVector());
 
</syntaxhighlight>
 
</syntaxhighlight>
  
* Required: Set the magnetic field map for the run:
+
* Or, the constraints can be created automatically by calling the following function:
  
 
<syntaxhighlight>
 
<syntaxhighlight>
DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
+
set<shared_ptr<DKinFitConstraint>> Create_Constraints(const DReactionVertexInfo* locReactionVertexInfo, const DReaction* locReaction, const DParticleCombo* locParticleCombo, const shared_ptr<const DKinFitChain>& locKinFitChain, DKinFitType locKinFitType, vector<shared_ptr<DKinFitConstraint_Vertex>>& locSortedVertexConstraints);
const DMagneticFieldMap* locMagneticFieldMap = locApplication->GetBfield(locRunNumber);
+
dKinFitter.Set_BField(locMagneticFieldMap);
+
 
</syntaxhighlight>
 
</syntaxhighlight>
  
* Optional: Set sizes for resource pools, set the debug flag, etc.
+
* 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:  
 
+
=== Fit Setup - New Event ===
+
* Call <span style="color:#0000FF">DKinFitter</span>::<span style="color:#008000">Reset_NewEvent</span>() to reset the fitter 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 <span style="color:#0000FF">DKinFitResults</span>, so they must be valid for the entire event.
+
 
+
=== Fit Setup - New Particle Combination ===
+
 
+
* Create <span style="color:#0000FF">DKinFitParticle</span> objects for your particles. These can be created by calling one of the following methods of <span style="color:#0000FF">DKinFitter_GlueX</span>:
+
 
+
 
<syntaxhighlight>
 
<syntaxhighlight>
const DKinFitParticle* DKinFitter_GlueX::Make_BeamParticle(const DBeamPhoton* locBeamPhoton);
+
void DKinFitConstraint_P4::Set_InitP3Guess(const TVector3& locInitP3Guess);
const DKinFitParticle* DKinFitter_GlueX::Make_DetectedParticle(const DKinematicData* locKinematicData);
+
void DKinFitConstraint_Vertex::Set_InitVertexGuess(const TVector3& locInitVertexGuess);
const DKinFitParticle* DKinFitter_GlueX::Make_DetectedShower(const DNeutralParticleHypothesis* locNeutralParticleHypothesis); //DO NOT call this unless the neutral is also in a vertex fit!
+
void DKinFitConstraint_Spacetime::Set_InitVertexGuess(const TVector3& locInitVertexGuess);
const DKinFitParticle* DKinFitter_GlueX::Make_DecayingParticle(Particle_t locPID);
+
void DKinFitConstraint_Spacetime::Set_InitTimeGuess(double locInitTimeGuess);
const DKinFitParticle* DKinFitter_GlueX::Make_MissingParticle(Particle_t locPID);
+
const DKinFitParticle* DKinFitter_GlueX::Make_TargetParticle(Particle_t locPID);
+
 
</syntaxhighlight>
 
</syntaxhighlight>
  
=== Fit Setup - New Fit ===
+
=== Performing the Kinematic 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.:
 
* 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.:
  
 
<syntaxhighlight>
 
<syntaxhighlight>
DKinFitter::Reset_NewFit();
+
void DKinFitter::Reset_NewFit(void);
 
</syntaxhighlight>
 
</syntaxhighlight>
  
* Create the constraints used for the fit. These can be created by calling one of the following functions:  
+
* Of the constraints created above, set the constraints you would like to apply for the current fit, by calling either of:  
  
 
<syntaxhighlight>
 
<syntaxhighlight>
DKinFitConstraint_Vertex* DKinFitter::Make_VertexConstraint(const deque<const DKinFitParticle*>& locFinalParticles, TVector3 locVertexGuess);
+
void DKinFitter::Add_Constraint(const shared_ptr<DKinFitConstraint>& locKinFitConstraint){dKinFitConstraints.insert(locKinFitConstraint);}
DKinFitConstraint_Vertex* DKinFitter::Make_VertexConstraint(const deque<const DKinFitParticle*>& locInitialParticles, const deque<const DKinFitParticle*>& locFinalParticles, TVector3 locVertexGuess);
+
void DKinFitter::Add_Constraints(const set<shared_ptr<DKinFitConstraint>>& locKinFitConstraints);
DKinFitConstraint_Spacetime* DKinFitter::Make_SpacetimeConstraint(const deque<const DKinFitParticle*>& locInitialParticles, const deque<const DKinFitParticle*>& locFinalParticles, bool locUseRFTimeFlag, TVector3 locVertexGuess, double locCommonTimeGuess);
+
//Note: below, locConstrainInitialParticleMassFlag is ignored if the parent particle is beam/detected/open-ended-decaying (enforce p4 constraint instead)
+
DKinFitConstraint_P4* DKinFitter::Make_P4Constraint(const deque<const DKinFitParticle*>& locInitialParticles, const deque<const DKinFitParticle*>& locFinalParticles, bool locConstrainInitialParticleMassFlag = true);
+
 
</syntaxhighlight>
 
</syntaxhighlight>
 
* Of the constraints created above, set the constraints you would like to apply for the current fit, by calling the appropriate one of the following functions:
 
 
<syntaxhighlight>
 
void DKinFitter::Set_Constraint(const DKinFitConstraint_P4* locConstraint);
 
void DKinFitter::Set_Constraint(const DKinFitConstraint_Vertex* locConstraint);
 
void DKinFitter::Set_Constraint(const DKinFitConstraint_Spacetime* locConstraint);
 
</syntaxhighlight>
 
 
=== Performing the Kinematic Fit ===
 
  
 
* Perform the kinematic fit by calling:
 
* Perform the kinematic fit by calling:
 
<syntaxhighlight>
 
<syntaxhighlight>
bool DKinFitter_GlueX::Fit_Reaction(void);
+
bool DKinFitter::Fit_Reaction(void);
 
</syntaxhighlight>
 
</syntaxhighlight>
  
Line 312: Line 189:
 
</syntaxhighlight>
 
</syntaxhighlight>
  
* 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 <span style="color:#0000FF">DKinFitParticle</span> and <span style="color:#0000FF">DKinFitConstraint</span> objects that were created earlier. However, each fit needs to update the "active" <span style="color:#0000FF">DKinFitParticle</span> and <span style="color:#0000FF">DKinFitConstraint</span> objects as each iteration improves the particle/vertex/etc. properties. Thus, in order to avoid over-writing the user-created <span style="color:#0000FF">DKinFitParticle</span> and <span style="color:#0000FF">DKinFitConstraint</span> objects, clones of these objects are made when <span style="color:#0000FF">DKinFitter</span>::<span style="color:#008000">Fit_Reaction</span>() 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 either:  
+
* 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 <span style="color:#0000FF">DKinFitParticle</span> and <span style="color:#0000FF">DKinFitConstraint</span> objects that were created earlier. However, each fit needs to update the "active" <span style="color:#0000FF">DKinFitParticle</span> and <span style="color:#0000FF">DKinFitConstraint</span> objects as each iteration improves the particle/vertex/etc. properties. Thus, in order to avoid over-writing the user-created <span style="color:#0000FF">DKinFitParticle</span> and <span style="color:#0000FF">DKinFitConstraint</span> objects, clones of these objects are made when <span style="color:#0000FF">DKinFitter</span>::<span style="color:#008000">Fit_Reaction</span>() 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:
  
 
<syntaxhighlight>
 
<syntaxhighlight>
const DKinFitParticle* DKinFitter::Get_OutputKinFitParticle(const DKinFitParticle* locInputKinFitParticle) const;
+
set<shared_ptr<DKinFitConstraint>> DKinFitter::Get_KinFitConstraints(void) const;
const DKinFitConstraint* DKinFitter::Get_OutputKinFitConstraint(const DKinFitConstraint* locInputKinFitConstraint) const;
+
set<shared_ptr<DKinFitParticle>> DKinFitter::Get_KinFitParticles(void) const;
 
</syntaxhighlight>
 
</syntaxhighlight>
  
* If you have lost track, the <span style="color:#0000FF">DKinFitParticle</span>'s can then be mapped to the original <span style="color:#0000FF">DKinematicData</span> objects via either:
+
* The output <span style="color:#0000FF">DKinFitParticle</span>'s can be mapped to the input <span style="color:#0000FF">DKinFitParticle</span>'s via:  
 
<syntaxhighlight>
 
<syntaxhighlight>
const DKinematicData* DKinFitter_GlueX::Get_Source_FromInput(const DKinFitParticle* locKinFitParticle) const;
+
shared_ptr<DKinFitParticle> DKinFitUtils::Get_InputKinFitParticle(shared_ptr<DKinFitParticle> locKinFitParticle) const;
const DKinematicData* DKinFitter_GlueX::Get_Source_FromOutput(const DKinFitParticle* locKinFitParticle) const;
+
</syntaxhighlight>
 +
 
 +
* And, the input <span style="color:#0000FF">DKinFitParticle</span>'s can be mapped to the original <span style="color:#0000FF">JObject</span>'s via:
 +
<syntaxhighlight>
 +
const JObject* DKinFitUtils_GlueX::Get_SourceJObject(shared_ptr<DKinFitParticle> locInputKinFitParticle) const;
 
</syntaxhighlight>
 
</syntaxhighlight>
  
 
* For a more detailed look at how the fit parameters have changed with respect to their input, measured values, the pulls can be retrieved:
 
* For a more detailed look at how the fit parameters have changed with respect to their input, measured values, the pulls can be retrieved:
 
<syntaxhighlight>
 
<syntaxhighlight>
//1st key is particle (NULL for rf-bunch), 2nd key is param type (defined in DKinFitParticle.h)
+
//1st key is particle, 2nd key is param type (defined in DKinFitParticle.h)
void DKinFitter::Get_Pulls(map<const DKinFitParticle*, map<DKinFitPullType, double> >& locPulls) const;  
+
void DKinFitter::Get_Pulls(map<shared_ptr<DKinFitParticle>, map<DKinFitPullType, double> >& locPulls) const;  
 
</syntaxhighlight>
 
</syntaxhighlight>
  
 
* The full covariance matrices for the fit parameters can be retrieved via:  
 
* The full covariance matrices for the fit parameters can be retrieved via:  
 
<syntaxhighlight>
 
<syntaxhighlight>
const TMatrixDSym* Get_VEta(void); //covariance matrix of the measured parameters used in the fit
+
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_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.  
+
const TMatrixDSym& Get_V(void); //combined covariance matrix of the measured and unknown parameters used in the fit.  
 
</syntaxhighlight>
 
</syntaxhighlight>
  
== Code Documentation (Under the Hood) - Resolving Constraints ==
+
== Notes on Manual Setup ==
  
When <span style="color:#0000FF">DKinFitter</span>::<span style="color:#008000">Fit_Reaction</span>() is called, there is a significant amount of work going on under the hood to resolve the constraints and prepare them for the kinematic fit. 
+
=== Magnetic Field ===
 +
* If the particle is charged, the magnetic field is taken into account if a common vertex is simultaneously kinematically fit
 +
** Magnetic field assumed to be constant over the propagation range
 +
* If a common vertex is not defined, the momentum is assumed to be constant
 +
* Energy loss between the common vertex and the particle position is ignored
  
=== Simple Example ===
+
=== Neutral Showers ===
 +
* For input neutral showers, the momentum is defined-by and updated-with the common vertex (see "Common Vertex Notes")
 +
* If a common vertex is not simultaneously kinematically fit for a neutral shower, the fit will fail
 +
* To include a neutral shower in a p4 fit without a vertex fit, create a neutral particle from it
 +
** With a full 7x7 covaraince matrix and input it instead of the neutral shower
 +
* Massive neutral showers (e.g. neutron) cannot be used in vertex constraints: only spacetime constraints. However, photons can.
 +
** This is because their momentum is defined by the vertex time
  
To illustrate this, consider the following reaction:
+
=== P4 Constraints ===
 +
* Don't do if studying an inclusive channel.
 +
* Pick the step of the decay chain that contains the missing or open-ended-decaying particle, to make sure it is constrained
 +
* If no missing or open-ended decaying particle, pick the initial step, to make sure the beam information is used
  
* &gamma; p &rarr; K<sup>+</sup> &Lambda;
+
=== Mass Constraints ===
** &Lambda; &rarr; p &pi;<sup>-</sup>
+
* Cannot use to constrain the mass of a missing particle: use a p4 constraint instead
  
When setting up a vertex-p4 fit, the user would create two <span style="color:#0000FF">DKinFitConstraint_P4</span> objects (one for each step) and two <span style="color:#0000FF">DKinFitConstraint_Vertex</span> objects (one for each vertex).  The kinematic fitter then needs to determine:  
+
* Mass of decaying particle constrained by either missing mass or invariant mass:
 +
** If no missing particles, can constrain either way
 +
*** By default, prefer using the invariant mass. This is to ensure consistency when combining with a p4 constraint.
 +
*** Failure to do so could lead to uninvertible matrices, and a failed kinematic fit.
 +
** If this particle has a missing decay product (inclusive or otherwise): Constrain by missing mass
 +
** If computation by missing mass would have a missing particle (inclusive or otherwise): Constrain by invariant mass
  
# Which <span style="color:#0000FF">DKinFitConstraint_Vertex</span> will define the <span style="color:red">&Lambda;</span> position, and which vertex will be CONSTRAINED-BY the <span style="color:red">&Lambda;</span>. Since two tracks are needed to define a vertex, it will choose the <span style="color:red">&Lambda;</span> decay and production vertex constraints for these, respectively.
+
* However, the way it is calculated is not defined in these constraints; it is defined when creating the decaying DKinFitParticle
# Which <span style="color:#0000FF">DKinFitConstraint_P4</span> object will be used to define the <span style="color:red">&Lambda;</span> four-momentum. Either step can work, but the momentum must be quoted at the defined <span style="color:red">&Lambda;</span>-position to stay self-consistent.
+
** So be sure to create them correctly!
## NOTE: This doesn't matter as much for the <span style="color:red">&Lambda;</span>, but charged decaying particles (such as the <span style="color:red">&Xi;<sup>-</sup></span>) bend in a magnetic field, and thus must have their momentum defined at the same position.
+
# Which <span style="color:#0000FF">DKinFitConstraint_P4</span> object will be used to use to constrain the four-momentum, and which to use to constrain the &Lambda; mass. It will choose the first and second steps for these, respectively.
+
  
=== Complex Example ===
+
=== Vertex Constraints ===
 +
* You can include the beam particle in the common vertex fit, but only do so if it has a non-zero error matrix.
 +
* Neutral and missing particles included in the constraint will not be used to constrain the vertex (no-constrain), but will be set with the fit common vertex
 +
** This is necessary if you want the neutral particle momentum (from an input shower) to change with the reconstructed vertex
 +
* Decaying particles should only be used to constrain a fit if the position is defined in another vertex constraint
 +
* Massive neutral showers (e.g. neutron) cannot be used in vertex constraints: only spacetime constraints. However, photons can.
 +
** This is because their momentum is defined by the vertex time
  
The above example was relatively simple and straight-forward. Now consider a more complicated reaction:
+
=== Spacetime Constraints ===
 +
* THESE ARE CURRENTLY DISABLED
 +
* It is not possible to fit a common time without simultaneously fitting a common vertex.
 +
** Requiring a common time at a non-common vertex has no meaning, and the fitter is not setup to input a pre-fit common vertex with uncertainties.
 +
* You can include the RF time in the fit by calling the Make_BeamParticle function with the RF information
 +
* Missing particles included in the constraint will not be used to constrain the spacetime (no-constrain), but will be set with the fit common spacetime
 +
* Decaying particles should only be used to constrain a fit if the spacetime is defined in another vertex constraint
  
* &gamma; p &rarr; K<sup>+</sup> K<sup>+</sup> &Xi;<sup>-</sup>
+
=== Objects & Memory ===
** &Xi;<sup>-</sup> &rarr; &pi;<sup>-</sup> &Lambda;
+
* User creates DKinFitParticles, adds them to the constraints
*** &Lambda; &rarr; p &pi;<sup>-</sup>
+
* When the fit starts, INTERNALLY ONLY, the particles are cloned and the constraints are cloned to contain the clones
 +
* User can reuse particles & constraints objects (since they are cloned, they are not overwritten)
  
Now, this example isn't particularly more difficult than the previous one; it just contains an extra step.  However, now consider the fact that the user could choose to treat any of the particles in this reaction as missing!  The kinematic fitter now needs to detect which direction to "traverse" the constraints to make sure that everything is properly defined in EVERY POSSIBLE CASE.  Detecting and resolving these cases is why the code is so complex. The possible cases are:
+
== 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 [https://halldsvn.jlab.org/repos/trunk/home/pmatt/plugins/kinfit_test/ Link]
* '''Missing <span style="color:red">K<sup>+</sup></span>''':
+
** Vertices/Momentum: The <span style="color:red">&Lambda;</span> decay step defines the <span style="color:red">&Lambda;</span> position and momentum, the <span style="color:red">&Lambda;</span> is used to constrain the <span style="color:red">&Xi;<sup>-</sup></span> decay vertex and momentum, and the <span style="color:red">&Xi;<sup>-</sup></span> is used to constrain the production vertex and missing-<span style="color:red">K<sup>+</sup></span> momentum.
+
** Masses/P4: The first step is used as the four-momentum constraint, and the following steps are used as invariant mass constraints.
+
 
+
* '''Missing <span style="color:red">&Lambda;</span>-decay product''': This is essentially the reverse of the above:
+
** Vertices/Momentum: The production step defines the <span style="color:red">&Xi;<sup>-</sup></span> momentum and position, the <span style="color:red">&Xi;<sup>-</sup></span> decay step is used to define the <span style="color:red">&Lambda;</span> position and momentum, and the final step is used to define the missing-decay-product position and momentum.
+
** Masses/P4: The first step is used as the four-momentum constraint, and the other steps are used as invariant-mass constraints.
+
 
+
* '''Missing <span style="color:red">&pi;<sup>-</sup></span> from <span style="color:red">&Xi;<sup>-</sup></span> decay''':
+
** Vertices/Momentum: The <span style="color:red">&Lambda;</span> decay step defines the <span style="color:red">&Lambda;</span> position and momentum, the production step defines the <span style="color:red">&Xi;<sup>-</sup></span> momentum and position, and the <span style="color:red">&Lambda;</span> and the <span style="color:red">&Xi;<sup>-</sup></span> define the missing <span style="color:red">&pi<sup>-</sup></span> position and momentum.
+
** Masses/P4: The first step is used as the four-momentum constraint, and the following steps are used as invariant mass constraints.
+
 
+
=== Inclusive Fits ===
+
 
+
Consider the following example of an inclusive kinematic fit:
+
 
+
* &gamma; p &rarr; K<sup>+</sup> &Lambda; X
+
** &Lambda; &rarr; p &pi;<sup>-</sup>
+
 
+
Where <span style="color:red">X</span> is an unknown particle.  In this case, two <span style="color:#0000FF">DKinFitConstraint_P4</span> objects would still be created, but neither would be used to apply a four-momentum-conservation constraint.  A 4-momentum constraint would have 4 constraints and 4 unknowns  (since the <span style="color:red">X</span> mass is unknown), which does not help constrain the system.
+
 
+
Finally, if the unknown particle <span style="color:red">X</span> is a decay product, then the mass of the decaying particle would be constrained via missing-mass rather than invariant mass.  At the moment I can't think of an example where this functionality would be desired, but it is still possible.  Otherwise, all mass constraints are by invariant mass rather than missing mass.
+
 
+
== Initial Vertex/Spacetime Guesses ==
+
Complex, multi-decay fits contain a lot of unknown fit parameters.  Choosing appropriate initial guesses for these parameters is essential for converging on the global &chi;<sup>2</sup> minimum.  To obtain initial guesses for the vertices/times, complex systems are broken down and each vertex is fit individually. For example, for a vertex-p4 fit of the following reaction:
+
 
+
* &gamma; p &rarr; K<sup>+</sup> K<sup>+</sup> &Xi;<sup>-</sup>
+
** &Xi;<sup>-</sup> &rarr; &pi;<sup>-</sup> &Lambda;
+
*** &Lambda; &rarr; p &pi;<sup>-</sup>
+
 
+
The initial guesses for the vertices are determined by fitting the reaction in several stages:
+
 
+
# &Lambda; decay vertex initial guess: Midpoint of the line between the decay-products (proton, &pi;<sup>-</sup>) at their POCA (position of closest approach) to each other.
+
# &Lambda; decay vertex better guess: &Lambda; &rarr; p &pi;<sup>-</sup> vertex fit, using the initial vertex guess.
+
# &Lambda; decay vertex best guess: &Lambda; &rarr; p &pi;<sup>-</sup> vertex-p4 fit, using the vertex from the previous fit as the initial guess.  
+
# &Xi;<sup>-</sup> decay vertex initial guess: Midpoint of the line between the decay-products (&Lambda; &pi;<sup>-</sup>) at their POCA to each other. The reconstructed &Lambda; information is derived from the earlier &Lambda; fits.
+
# &Xi;<sup>-</sup> decay vertex better guess: &Xi;<sup>-</sup> &rarr; &pi;<sup>-</sup> &Lambda; vertex fit, using the initial vertex guess.  The reconstructed &Lambda; information is derived from the earlier &Lambda; fits.
+
# &Xi;<sup>-</sup> decay vertex best guess: &Xi;<sup>-</sup> &rarr; &pi;<sup>-</sup> &Lambda; vertex fit, using the initial vertex guess.  The reconstructed &Lambda; information is derived from the earlier &Lambda; fits.
+
 
+
# &Xi;<sup>-</sup> &rarr; &pi;<sup>-</sup> &Lambda; vertex fit, using the results for the &Lambda; from the previous fit.
+
# &Xi;<sup>-</sup> &rarr; &pi;<sup>-</sup> &Lambda; vertex-p4 fit, using the results for the &Lambda; from the previous fit, and the vertex from the previous fit as the initial guess.
+
# &gamma; p &rarr; K<sup>+</sup> K<sup>+</sup> &Xi;<sup>-</sup> vertex fit, using the results for the &Xi;<sup>-</sup> from the previous fit.
+

Latest revision as of 08:32, 3 November 2017

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). So, to kinematically fit a DReaction, you only have to specify the desired kinematic fit type in the DReaction.
    • 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

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

DKinFitUtils_GlueX

  • 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/DReaction.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 an inclusive channel, or if a 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 "measured" DNeutralParticleHypothesis objects in the combos sets the particle momentum assuming that the particle came from the position roughly defined during comboing (see DSourceComboVertexer).
  • 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.

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 primary dependency is ROOT, from which it uses the vector and matrix classes. It also depends on DResettable and DResourcePool in sim-recon/src/libraries/include/ to help with memory management. 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):
//FOR SETUP:
virtual bool Get_IncludeBeamlineInVertexFitFlag(void) const = 0;
 
//FOR FITS:
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.
    • These classes are not necessary to perform the fit, however, they can be useful for using some of the (optional) DKinFitUtils functions to help setup the constraints for your fits.

In addition, resource pools of TMatrixFSym objects are used to manage memory for the covariance matrices.

Fit Setup - New Run

  • Required: First, create pointers to the DKinFitter & DKinFitUtils_GlueX classes in your header file:
#include "KINFITTER/DKinFitter.h"
#include "ANALYSIS/DKinFitUtils_GlueX.h"
 
//in class definition:
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

  • At the beginning of each event, the kinematic fitter (and the utilities) need to be reset This recycles the memory (& output/input object mapping) of all objects for reuse. After the fit, these fit results objects are stored in DKinFitResults, so they must stay valid for the entire event. To reset everything, call:
void DKinFitUtils_GlueX::Reset_NewEvent(void);
void DKinFitter::Reset_NewEvent(void);

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:
shared_ptr<DKinFitParticle> Make_BeamParticle(const DBeamPhoton* locBeamPhoton);
shared_ptr<DKinFitParticle> Make_BeamParticle(const DBeamPhoton* locBeamPhoton, const DEventRFBunch* locEventRFBunch); //sets rf time for photon
shared_ptr<DKinFitParticle> Make_TargetParticle(Particle_t locPID);
 
shared_ptr<DKinFitParticle> Make_DetectedParticle(const DKinematicData* locKinematicData);
shared_ptr<DKinFitParticle> Make_DetectedShower(const DNeutralShower* locNeutralShower, Particle_t locPID); //DO NOT call this unless the neutral is also in a vertex fit!
 
shared_ptr<DKinFitParticle> Make_DecayingParticle(Particle_t locPID, const set<shared_ptr<DKinFitParticle>>& locFromInitialState, const set<shared_ptr<DKinFitParticle>>& locFromFinalState);
shared_ptr<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

  • Create the constraints used for the fit. These can be created manually by calling the following functions on DKinFitUtils_GlueX:
shared_ptr<DKinFitConstraint_Mass> Make_MassConstraint(const shared_ptr<DKinFitParticle>& locDecayingParticle);
shared_ptr<DKinFitConstraint_P4> Make_P4Constraint(const set<shared_ptr<DKinFitParticle>>& locInitialParticles, const set<shared_ptr<DKinFitParticle>>& locFinalParticles);
shared_ptr<DKinFitConstraint_Vertex> Make_VertexConstraint(const set<shared_ptr<DKinFitParticle>>& locFullConstrainParticles, const set<shared_ptr<DKinFitParticle>>& locNoConstrainParticles, TVector3 locVertexGuess = TVector3());
shared_ptr<DKinFitConstraint_Spacetime> Make_SpacetimeConstraint(const set<shared_ptr<DKinFitParticle>>& locFullConstrainParticles, const set<shared_ptr<DKinFitParticle>>& locOnlyConstrainTimeParticles,
	const set<shared_ptr<DKinFitParticle>>& locNoConstrainParticles, TLorentzVector locSpacetimeGuess = TLorentzVector());
  • Or, the constraints can be created automatically by calling the following function:
set<shared_ptr<DKinFitConstraint>> Create_Constraints(const DReactionVertexInfo* locReactionVertexInfo, const DReaction* locReaction, const DParticleCombo* locParticleCombo, const shared_ptr<const DKinFitChain>& locKinFitChain, DKinFitType locKinFitType, vector<shared_ptr<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);

Performing the Kinematic 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.:
void DKinFitter::Reset_NewFit(void);
  • Of the constraints created above, set the constraints you would like to apply for the current fit, by calling either of:
void DKinFitter::Add_Constraint(const shared_ptr<DKinFitConstraint>& locKinFitConstraint){dKinFitConstraints.insert(locKinFitConstraint);}
void DKinFitter::Add_Constraints(const set<shared_ptr<DKinFitConstraint>>& locKinFitConstraints);
  • 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<shared_ptr<DKinFitConstraint>> DKinFitter::Get_KinFitConstraints(void) const;
set<shared_ptr<DKinFitParticle>> DKinFitter::Get_KinFitParticles(void) const;
  • The output DKinFitParticle's can be mapped to the input DKinFitParticle's via:
shared_ptr<DKinFitParticle> DKinFitUtils::Get_InputKinFitParticle(shared_ptr<DKinFitParticle> locKinFitParticle) const;
  • And, the input DKinFitParticle's can be mapped to the original JObject's via:
const JObject* DKinFitUtils_GlueX::Get_SourceJObject(shared_ptr<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<shared_ptr<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.

Notes on Manual Setup

Magnetic Field

  • If the particle is charged, the magnetic field is taken into account if a common vertex is simultaneously kinematically fit
    • Magnetic field assumed to be constant over the propagation range
  • If a common vertex is not defined, the momentum is assumed to be constant
  • Energy loss between the common vertex and the particle position is ignored

Neutral Showers

  • For input neutral showers, the momentum is defined-by and updated-with the common vertex (see "Common Vertex Notes")
  • If a common vertex is not simultaneously kinematically fit for a neutral shower, the fit will fail
  • To include a neutral shower in a p4 fit without a vertex fit, create a neutral particle from it
    • With a full 7x7 covaraince matrix and input it instead of the neutral shower
  • Massive neutral showers (e.g. neutron) cannot be used in vertex constraints: only spacetime constraints. However, photons can.
    • This is because their momentum is defined by the vertex time

P4 Constraints

  • Don't do if studying an inclusive channel.
  • Pick the step of the decay chain that contains the missing or open-ended-decaying particle, to make sure it is constrained
  • If no missing or open-ended decaying particle, pick the initial step, to make sure the beam information is used

Mass Constraints

  • Cannot use to constrain the mass of a missing particle: use a p4 constraint instead
  • Mass of decaying particle constrained by either missing mass or invariant mass:
    • If no missing particles, can constrain either way
      • By default, prefer using the invariant mass. This is to ensure consistency when combining with a p4 constraint.
      • Failure to do so could lead to uninvertible matrices, and a failed kinematic fit.
    • If this particle has a missing decay product (inclusive or otherwise): Constrain by missing mass
    • If computation by missing mass would have a missing particle (inclusive or otherwise): Constrain by invariant mass
  • However, the way it is calculated is not defined in these constraints; it is defined when creating the decaying DKinFitParticle
    • So be sure to create them correctly!

Vertex Constraints

  • You can include the beam particle in the common vertex fit, but only do so if it has a non-zero error matrix.
  • Neutral and missing particles included in the constraint will not be used to constrain the vertex (no-constrain), but will be set with the fit common vertex
    • This is necessary if you want the neutral particle momentum (from an input shower) to change with the reconstructed vertex
  • Decaying particles should only be used to constrain a fit if the position is defined in another vertex constraint
  • Massive neutral showers (e.g. neutron) cannot be used in vertex constraints: only spacetime constraints. However, photons can.
    • This is because their momentum is defined by the vertex time

Spacetime Constraints

  • THESE ARE CURRENTLY DISABLED
  • It is not possible to fit a common time without simultaneously fitting a common vertex.
    • Requiring a common time at a non-common vertex has no meaning, and the fitter is not setup to input a pre-fit common vertex with uncertainties.
  • You can include the RF time in the fit by calling the Make_BeamParticle function with the RF information
  • Missing particles included in the constraint will not be used to constrain the spacetime (no-constrain), but will be set with the fit common spacetime
  • Decaying particles should only be used to constrain a fit if the spacetime is defined in another vertex constraint

Objects & Memory

  • User creates DKinFitParticles, adds them to the constraints
  • When the fit starts, INTERNALLY ONLY, the particles are cloned and the constraints are cloned to contain the clones
  • User can reuse particles & constraints objects (since they are cloned, they are not overwritten)

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