Difference between revisions of "Mattione GlueX Kinematic Fitting"

From GlueXWiki
Jump to: navigation, search
(DKinFitter Data Classes)
(Fit Setup - New Event)
Line 248: Line 248:
  
 
=== Fit Setup - New Event ===
 
=== 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.  
+
* 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.  
# Creating <span style="color:#0000FF">DKinFitParticle</span> objects for your particles.
+
 
 +
=== 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>
 +
const DKinFitParticle* DKinFitter_GlueX::Make_BeamParticle(const DBeamPhoton* locBeamPhoton);
 +
const DKinFitParticle* DKinFitter_GlueX::Make_DetectedParticle(const DKinematicData* locKinematicData);
 +
const DKinFitParticle* DKinFitter_GlueX::Make_DetectedShower(const DNeutralParticleHypothesis* locNeutralParticleHypothesis); //DO NOT call this unless the neutral is also in a vertex fit!
 +
const DKinFitParticle* DKinFitter_GlueX::Make_DecayingParticle(Particle_t locPID);
 +
const DKinFitParticle* DKinFitter_GlueX::Make_MissingParticle(Particle_t locPID);
 +
const DKinFitParticle* DKinFitter_GlueX::Make_TargetParticle(Particle_t locPID);
 +
</syntaxhighlight>
  
 
=== Fit Setup - New Fit - User Input ===
 
=== Fit Setup - New Fit - User Input ===

Revision as of 11:37, 19 April 2015

Quick Start

  • The kinematic fitter (DKinFitter_GlueX) is fully integrated into the analysis libraries (interfaced to DKinFitResults_factory), 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.

Support

  • Note that spacetime constraints and fits including final-state massive neutral particles are not yet supported.

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)
  • Invariant mass constraints: 1 constraint (E2 - p2 - m2 = 0)
    • Note that by default, the ANALYSIS library does not constrain the masses of the ω or φ resonances, even though their widths are very narrow. They can be constrained though if you setup the kinematic fit manually.
  • Vertex constraints: 3 unknowns per vertex (x, y, z), 2 constraints per track
    • 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 track (vertex + time)

DKinFitResuls_factory

  • This factory sets up and performs 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
  • This factory automatically determines what p4, vertex, and spacetime (not yet supported!) constraints to setup such that each the data in each DParticleCombo are maximally constrained.
    • The factory will 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.
    • The factory will also automatically exclude resonances from the p4 constraints of kinematic fits, as long as the IsFixedMass() 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 additional, non-resonant particles from the p4 constraints of a kinematic fit, call the DReactionStep::Set_ApplyKinFitMassConstraintOnInitialParticleFlag(bool) function with an argument of true on the step where the particle is in the initial state.

Important Notes

Detected Neutral Particles

  • The DNeutralParticleHypothesis objects from the default factory sets the particle momentum assuming that the particle came from the center of the target.
  • 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.
  • Choosing to use the neutral shower instead of the default DNeutralParticleHypothesis object is toggled by the boolean flag when calling the DKinFitter_GlueX::Make_DetectedParticle(const DNeutralParticleHypothesis*, bool) function.

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.
    • 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:
    • 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
    • This is useful for when you do not want to use a mass constraint on the decaying particle.

Manual Kinematic Fit Example

  • 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 γ, p → K+, Λ; Λ → p, π-

DKinFitAction_Sample.h

#ifndef _DKinFitAction_Sample_
#define _DKinFitAction_Sample_
 
#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_Sample : public DAnalysisAction
{
  public:
    DKinFitAction_Sample(const DReaction* locReaction, string locActionUniqueString = "") : 
    DAnalysisAction(locReaction, "KinFit_Sample", false, locActionUniqueString){}
 
  private:
    bool Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo);
    void Initialize(JEventLoop* locEventLoop);
 
    DKinFitter_GlueX dKinFitter;
 
    //define any histograms here
};
 
#endif // _DKinFitAction_Sample_

DKinFitAction_Sample.cc

#include "DKinFitAction_Sample.h"
 
void DKinFitAction_Sample::Initialize(JEventLoop* locEventLoop)
{
  DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
 
  double locTargetZCenter = 65.0;
  DGeometry* locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
  locGeometry->GetTargetZ(locTargetZCenter);
 
  //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)
    dKinFitter.Set_BField(locMagneticFieldMap);
 
  //create any histograms here
}
 
bool DKinFitAction_Sample::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)
 
  //CREATE DKINFITPARTICLE OBJECTS FOR EACH PARTICLE
  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);
 
  // 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
}
  • Somewhere in the DReaction_factory in your plugin:
#include "DKinFitAction_Sample.h"
//...
  locReaction->Add_AnalysisAction(new DKinFitAction_Sample(locReaction));

Code Documentation

DKinFitter_GlueX and DKinFitter

  • The DKinFitter 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:
    • Particles: px, py, pz, xx, xy, xz, t
    • Showers: E, xx, xy, xz, t
  • However, the DKinFitter class knows nothing about the presence, strength, or direction of the magnetic field. The following DKinFitter member functions are used during the fit to access the magnetic field properties, but they have been declared as pure-virtual:
virtual bool DKinFitter::Get_IsBFieldNearBeamline(void) const = 0;
virtual TVector3 DKinFitter::Get_BField(const TVector3& locPosition) const = 0;
  • Thus DKinFitter is an abstract base class; it cannot be instantiated directly. Instead, for each experiment a class should inherit from DKinFitter and define these methods; for GlueX, this is the DKinFitter_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 DKinFitter 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: 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.
  • DKinFitConstraint_VertexBase: An abstract base class for the vertex and spacetime constraints. Particles at the vertex are contained in one of two different member variables:
    • dFullConstrainParticles: Particles that are used to constrain the vertex/spacetime position (e.g. charged, decaying, or beam).
    • dNoConstrainParticles: Particles whose vertex/spacetime are DEFINED BY this constraint, rather than used to constrain it (e.g. missing particles, decaying particles, neutral showers).
  • 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. In addition to the particles defined in DKinFitConstraint_VertexBase, it contains a new class of particles:
    • dOnlyConstrainTimeParticles: Particles who are used to constrain the time, but not the vertex (e.g. neutral showers).

Fit Setup - New Run

  • Required: First, create an instance of the DKinFitter_GlueX class in your header file:
DKinFitter_GlueX dKinFitter;
  • Required: Set the magnetic field map for the run:
DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
const DMagneticFieldMap* locMagneticFieldMap = locApplication->GetBfield(locRunNumber);
dKinFitter.Set_BField(locMagneticFieldMap);
  • Optional: Set sizes for resource pools, set the debug flag, etc.

Fit Setup - New Event

  • Call DKinFitter::Reset_NewEvent() 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 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 DKinFitter_GlueX:
const DKinFitParticle* DKinFitter_GlueX::Make_BeamParticle(const DBeamPhoton* locBeamPhoton);
const DKinFitParticle* DKinFitter_GlueX::Make_DetectedParticle(const DKinematicData* locKinematicData);
const DKinFitParticle* DKinFitter_GlueX::Make_DetectedShower(const DNeutralParticleHypothesis* locNeutralParticleHypothesis); //DO NOT call this unless the neutral is also in a vertex fit!
const DKinFitParticle* DKinFitter_GlueX::Make_DecayingParticle(Particle_t locPID);
const DKinFitParticle* DKinFitter_GlueX::Make_MissingParticle(Particle_t locPID);
const DKinFitParticle* DKinFitter_GlueX::Make_TargetParticle(Particle_t locPID);

Fit Setup - New Fit - User Input

The kinematic fitter supports conducting many fits with many different constraints in a given event. This includes re-using the DKinFitParticle and DKinFitConstraint objects that you've created. 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 either:

const DKinFitParticle* DKinFitter::Get_OutputKinFitParticle(const DKinFitParticle* locInputKinFitParticle) const;
const DKinFitConstraint* DKinFitter::Get_OutputKinFitConstraint(const DKinFitConstraint* locInputKinFitConstraint) const;

If you have lost track, the DKinFitParticle's can then be mapped to the original DKinematicData objects via either:

const DKinematicData* DKinFitter_GlueX::Get_Source_FromInput(const DKinFitParticle* locKinFitParticle) const;
const DKinematicData* DKinFitter_GlueX::Get_Source_FromOutput(const DKinFitParticle* locKinFitParticle) const;


  1. Creating DKinFitConstraint objects for your constraints, with your DKinFitParticle's as input.
  2. Resetting the fitter for each new fit. This clears the "current-fit" constraints / particles / matrices variables. Previous-fit objects are still valid.
  3. Setting the constraints for your current fit.
  4. Call DKinFitter::Fit_Reaction();

Fit Setup - New Fit - Internal Setup

At this stage, the fit has all of the user-input it needs: the constraints to apply, and the particles to constrain them with (members of the constraints).

  1. UBER-EXPERTS ONLY: Sort a constraints.

DKinFitResults_factory.cc