Difference between revisions of "Mattione GlueX Kinematic Fitting"

From GlueXWiki
Jump to: navigation, search
(DKinFitResuls_factory)
Line 3: Line 3:
 
** See [[Mattione_GlueX_Analysis_Factories | GlueX Analysis Factories]] for details on the analysis libraries.
 
** 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.
 
== Currently Supported ==
 
* Multiple, simultaneous 4-momentum, invariant mass, and vertex constraints, except for the cases outlined in the "Currently Unsupported" section.
 
 
== Currently Unsupported ==
 
* P4 or vertex constraints with "few" detected particles.
 
** e.g. &gamma;, p &rarr; K<sup>+</sup>, &Lambda;; &Lambda; &rarr; p, &pi;<sup>-</sup> works fine, but it doesn't work if any one of the charged particles is missing.
 
* Performing any spacetime fits.
 
  
 
== Constraints and Unknowns ==
 
== Constraints and Unknowns ==
Line 20: Line 12:
 
*** &Delta;'''x''' = '''x'''<sub>Measured</sub> - '''x'''<sub>Unknown</sub>
 
*** &Delta;'''x''' = '''x'''<sub>Measured</sub> - '''x'''<sub>Unknown</sub>
 
* Spacetime constraints: 4 unknowns per spacetime point (x, y, z, t), 3 constraints per track (vertex + time)
 
* Spacetime constraints: 4 unknowns per spacetime point (x, y, z, t), 3 constraints per track (vertex + time)
 +
 +
== Supported / Unsupported ==
 +
 +
=== Supported ===
 +
* Multiple, simultaneous 4-momentum, invariant mass, and vertex constraints, except for the cases outlined in the "Currently Unsupported" section.
 +
 +
=== Unsupported ===
 +
* P4 or vertex constraints with "few" detected particles.
 +
** e.g. &gamma;, p &rarr; K<sup>+</sup>, &Lambda;; &Lambda; &rarr; p, &pi;<sup>-</sup> works fine, but it doesn't work if any one of the charged particles is missing.
 +
* Spacetime fits.
  
 
== DKinFitResuls_factory ==
 
== DKinFitResuls_factory ==

Revision as of 23:45, 16 February 2013

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.

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)
  • 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: 4 unknowns per spacetime point (x, y, z, t), 3 constraints per track (vertex + time)

Supported / Unsupported

Supported

  • Multiple, simultaneous 4-momentum, invariant mass, and vertex constraints, except for the cases outlined in the "Currently Unsupported" section.

Unsupported

  • P4 or vertex constraints with "few" detected particles.
    • e.g. γ, p → K+, Λ; Λ → p, π- works fine, but it doesn't work if any one of the charged particles is missing.
  • Spacetime fits.

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.
  • To exclude additional, non-resonant particles from the p4 constraints of a kinematic fit, call the DReaction::Exclude_DecayingParticleFromP4KinFit() function with the index of DReactionStep where the particle you want to exclude from the fit is the initial particle (i.e. it's decay step).

Important Note Regarding 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.

Important Note Regarding 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.cc or DVertex_factory.cc 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.

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, const deque<pair<const DParticleCombo*, bool> >& locPreviousParticleCombos);
    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());
  dKinFitter.Set_BField(locApplication->GetBfield()); //need call only once
 
  //create any histograms here
}
 
bool DKinFitAction_Sample::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, const deque<pair<const DParticleCombo*, bool> >& locPreviousParticleCombos)
{
  //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);
  dKinFitter.Add_P4Constraint(locInitialKinFitParticles, locFinalKinFitParticles);
 
  // second p4 constraint:
  locInitialKinFitParticles.clear();
  locFinalKinFitParticles.clear();
  locInitialKinFitParticles.push_back(locKinFitParticle_Lambda);
  locFinalKinFitParticles.push_back(locKinFitParticle_PiMinus);
  locFinalKinFitParticles.push_back(locKinFitParticle_Proton);
  dKinFitter.Add_P4Constraint(locInitialKinFitParticles, locFinalKinFitParticles);
 
  // 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);
  dKinFitter.Add_VertexConstraint(locInitialKinFitParticles, locFinalKinFitParticles, locVertexGuess);
 
  // 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));