Difference between revisions of "Mattione GlueX Kinematic Fitting"
From GlueXWiki
(→Manual Kinematic Fit Example) |
|||
Line 62: | Line 62: | ||
private: | private: | ||
bool Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, const deque<pair<const DParticleCombo*, bool> >& locPreviousParticleCombos); | 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 | //define any histograms here | ||
Line 73: | Line 75: | ||
<syntaxhighlight> | <syntaxhighlight> | ||
#include "DKinFitAction_Sample.h" | #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) | bool DKinFitAction_Sample::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, const deque<pair<const DParticleCombo*, bool> >& locPreviousParticleCombos) | ||
{ | { | ||
Line 79: | Line 90: | ||
//Lambda -> p, pi- | //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 | //CREATE DKINFITPARTICLE OBJECTS FOR EACH PARTICLE | ||
const DBeamPhoton* locBeamPhoton = static_cast<const DBeamPhoton*>(locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle_Measured()); | const DBeamPhoton* locBeamPhoton = static_cast<const DBeamPhoton*>(locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle_Measured()); | ||
− | const DKinFitParticle* locKinFitParticle_BeamPhoton = | + | const DKinFitParticle* locKinFitParticle_BeamPhoton = dKinFitter.Make_BeamParticle(locBeamPhoton); |
− | const DKinFitParticle* locKinFitParticle_Target = | + | const DKinFitParticle* locKinFitParticle_Target = dKinFitter.Make_TargetParticle(Proton); |
− | const DKinFitParticle* locKinFitParticle_Lambda = | + | 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 DChargedTrackHypothesis* locChargedTrackHypothesis_KPlus = static_cast<const DChargedTrackHypothesis*>(locParticleCombo->Get_ParticleComboStep(0)->Get_FinalParticle_Measured(0)); | ||
− | const DKinFitParticle* locKinFitParticle_KPlus = | + | 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 DChargedTrackHypothesis* locChargedTrackHypothesis_Proton = static_cast<const DChargedTrackHypothesis*>(locParticleCombo->Get_ParticleComboStep(1)->Get_FinalParticle_Measured(0)); | ||
− | const DKinFitParticle* locKinFitParticle_Proton = | + | 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 DChargedTrackHypothesis* locChargedTrackHypothesis_PiMinus = static_cast<const DChargedTrackHypothesis*>(locParticleCombo->Get_ParticleComboStep(1)->Get_FinalParticle_Measured(1)); | ||
− | const DKinFitParticle* locKinFitParticle_PiMinus = | + | const DKinFitParticle* locKinFitParticle_PiMinus = dKinFitter.Make_DetectedParticle(locChargedTrackHypothesis_PiMinus); |
// SETUP THE CONSTRAINTS | // SETUP THE CONSTRAINTS | ||
− | + | dKinFitter.Reset_NewFit(); //disregards the constraints from the previous kinematic fit | |
deque<const DKinFitParticle*> locInitialKinFitParticles, locFinalKinFitParticles; | deque<const DKinFitParticle*> locInitialKinFitParticles, locFinalKinFitParticles; | ||
Line 109: | Line 117: | ||
locFinalKinFitParticles.push_back(locKinFitParticle_KPlus); | locFinalKinFitParticles.push_back(locKinFitParticle_KPlus); | ||
locFinalKinFitParticles.push_back(locKinFitParticle_Lambda); | locFinalKinFitParticles.push_back(locKinFitParticle_Lambda); | ||
− | + | dKinFitter.Add_P4Constraint(locInitialKinFitParticles, locFinalKinFitParticles); | |
// second p4 constraint: | // second p4 constraint: | ||
Line 117: | Line 125: | ||
locFinalKinFitParticles.push_back(locKinFitParticle_PiMinus); | locFinalKinFitParticles.push_back(locKinFitParticle_PiMinus); | ||
locFinalKinFitParticles.push_back(locKinFitParticle_Proton); | locFinalKinFitParticles.push_back(locKinFitParticle_Proton); | ||
− | + | dKinFitter.Add_P4Constraint(locInitialKinFitParticles, locFinalKinFitParticles); | |
// vertex constraint | // vertex constraint | ||
Line 126: | Line 134: | ||
locFinalKinFitParticles.push_back(locKinFitParticle_PiMinus); | locFinalKinFitParticles.push_back(locKinFitParticle_PiMinus); | ||
locFinalKinFitParticles.push_back(locKinFitParticle_Proton); | locFinalKinFitParticles.push_back(locKinFitParticle_Proton); | ||
− | + | dKinFitter.Add_VertexConstraint(locInitialKinFitParticles, locFinalKinFitParticles, locVertexGuess); | |
// PERFORM THE KINEMATIC FIT | // PERFORM THE KINEMATIC FIT | ||
− | + | dKinFitter.Fit_Reaction(); | |
// GET THE FIT RESULTS | // GET THE FIT RESULTS | ||
− | double locConfidenceLevel = | + | double locConfidenceLevel = dKinFitter.Get_ConfidenceLevel(); |
− | double locChiSq = | + | double locChiSq = dKinFitter.Get_ChiSq(); |
− | unsigned int locNDF = | + | unsigned int locNDF = dKinFitter.Get_NDF(); |
− | const TMatrixDSym* locMeasuredParametersKinFitCovarianceMatrix = | + | const TMatrixDSym* locMeasuredParametersKinFitCovarianceMatrix = dKinFitter.Get_VEta(); |
− | const TMatrixDSym* locUnknownParametersCovarianceMatrix = | + | const TMatrixDSym* locUnknownParametersCovarianceMatrix = dKinFitter.Get_VXi(); |
− | DKinFitParticle* locKinFitParticle_Proton_Output = | + | const DKinFitParticle* locKinFitParticle_Proton_Output = dKinFitter.Get_OutputKinFitParticle(locKinFitParticle_Proton); |
− | DKinFitParticle* locKinFitParticle_PiMinus_Output = | + | const DKinFitParticle* locKinFitParticle_PiMinus_Output = dKinFitter.Get_OutputKinFitParticle(locKinFitParticle_PiMinus); |
− | DKinFitParticle* locKinFitParticle_KPlus_Output = | + | const DKinFitParticle* locKinFitParticle_KPlus_Output = dKinFitter.Get_OutputKinFitParticle(locKinFitParticle_KPlus); |
− | DKinFitParticle* locKinFitParticle_Lambda_Output = | + | const DKinFitParticle* locKinFitParticle_Lambda_Output = dKinFitter.Get_OutputKinFitParticle(locKinFitParticle_Lambda); |
− | DKinFitParticle* locKinFitParticle_BeamPhoton_Output = | + | 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 | 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) | 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]; | double locProtonPxPull = (locPulls[locChargedTrackHypothesis_Proton])[d_PxPull]; | ||
Line 152: | Line 160: | ||
return (locConfidenceLevel >= 0.01); //cut if < 1% confidence level | return (locConfidenceLevel >= 0.01); //cut if < 1% confidence level | ||
} | } | ||
− | |||
</syntaxhighlight> | </syntaxhighlight> |
Revision as of 15:15, 14 October 2012
Contents
Quick Start
- The kinematic fitter 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.
Currently Supported
- Multiple, simultaneous 4-momentum constraints 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. γ, p → K+, Λ; Λ → p, π- works fine, but it doesn't work if any one of the charged particles is missing.
- Using decaying particles to constrain vertices together
- Performing any 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.
- 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.
- 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 }