Brief Kinematic Fitting Example

From GlueXWiki
Revision as of 18:54, 19 March 2012 by Wilevine (Talk | contribs) (Brief summary)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

Brief summary

Kinematic fitting in the HallD sim-recon code is provided by the DKinFit class.

DKinFit kfit;

In the simplest mode, the fitter takes a fully specified reaction (4-momenta and covariance matrices for all initial and final state particles) and fits using the constraint of overall conservation of 4-momentum. Missing final state particles can be specified, and additional mass constraints can be applied. The results of the fit can be used to filter out background that does not satisfy the fit hypothesis and to improve the 4-momenta and covariance matrices.

The kinematic fitter takes as its input two vectors of DKinematicData objects, containing the initial and final state particles, specified by SetInitial() and SetFinal(). Each DKinematicData contains both a 4-momentum and a kinematic covariance matrix.

vector<DKinematicData> initial;
vector<DKinematicData> final;

DKinematicData beam_photon,target,reconstructed_proton,reconstructed_pi1,...;
//fill these objects appropriately using objects (DChargedTrackHypothesis, DNeutralParticleHypothesis, DBeamPhoton) from factories
//...

initial.push_back(beam_photon);
initial.push_back(target);

final.push_back(reconstructed_proton);
final.push_back(reconstructed_pi1);
final.push_back(...);

kfit.SetFinal(final);
kfit.SetInitial(initial);

If there is a missing particle in the final state, use SetMissingParticle():

const double neutron_mass=.940;
kfit.SetMissingParticle(neutron_mass);

If there is a mass constraint (e.g. two photons from a pi0 decay), use SetMassConstraint():

vector<DKinematicData> final(3);
final[0] = reconstructed_proton;
final[1] = reconstructed_photon1;
final[2] = reconstructed_photon2;

vector<int> constrained_particles_indices(2);
//1 and 2 are the indices of the photons in final
constrained_particles_indices[0]=1;
constrained_particles_indices[1]=2;

const double pi0_mass=0.135;
kfit.SetMassConstraint(pi0_mass,constrained_particles_indices);

kfit.SetFinal(final);

After specifying initial and final state particles, missing particles, and mass constraints, simply use the Fit() function:

kfit.Fit();

Useful functions for accessing the results of the fit:

  • Prob(): returns the confidence level of fit, based on \chi ^{2}
  • GetPull(int n): returns the n'th pull. There is one pull for each measured quantity, so three pulls for each measured particle (p_x, p_y, and p_z). Pulls 0-2 are the pulls of the first initial particle (0 for p_x, 1 for p_y, 2 for p_z), pulls 3-5 are the pulls of the second initial particle, pulls 6-8 are the pulls of the first final state particle, etc.
  • GetInitial_out(), GetFinal_out(): these two functions return vectors of DKinematicData with post-fit "corrected" 4-momentum and covariance matrices.
  • GetMissingError(): returns the error on the missing particle momentum. Unfortunately, it uses the (p,lambda,phi) parameterization at the moment for some reason, so it's somewhat useless.

b1pi example

The code below is a plugin designed to analyze a sample of b1pi events using the kinematic fitter ( see here for how to generate such events).

It loops thru all detected charged tracks and photons to produce combinations that correspond to a p pi+ pi+ pi- pi- gamma gamma final state. For each combination, the kinematic fit is applied. For each event, only the fit with the highest confidence is examined. If the confidence level is about 1%, some results of the fit are saved to histograms.

DEventProcessor_b1pi_kinfit_simple.cc

#include <iostream>
using namespace std;

#include <TDirectoryFile.h>

#include <JANA/JApplication.h>
using namespace jana;

#include <PID/DNeutralParticleHypothesis.h>
#include <PID/DBeamPhoton.h>
#include <PID/DKinFit.h>
#include <PID/DChargedTrackHypothesis.h>
#include <units.h>

#include "DEventProcessor_b1pi_kinfit_simple.h"
//==========================================================================
// This file contains code for a plugin that will create a few histograms
// based on the reconstructed data. It can mainly serve as an example
// of how one could access the reconstructed data in their own analysis
// program.
//
// Histograms are created in the init() method and they are filled in
// the evnt() method.
//==========================================================================

// Routine used to create our JEventProcessor
extern "C"{
  void InitPlugin(JApplication *app){
    InitJANAPlugin(app);
    app->AddProcessor(new DEventProcessor_b1pi_kinfit_simple());
  }
}

//------------------
// init
//------------------
jerror_t DEventProcessor_b1pi_kinfit_simple::init(void)
{

  // Create KIN_FIT directory
  TDirectory *dir = new TDirectoryFile("KIN_FIT","KIN_FIT");
  dir->cd();

  // Create histograms
  for (int j=0;j<3*n_final_particles;j++){
    stringstream name;
    name << "h_pulls_" << j;
    h_pulls[j] = new TH1F(name.str().c_str(),name.str().c_str(),50,-5.0,5.0);
  }

  h_cl = new TH1F("h_cl","Confidence Level",100,0.0,1.0);

  h_mass_2pip_2pim_2gamma = new TH1F("h_mass_2pip_2pim_2gamma",
         "#pi^{+}#pi^{+}#pi^{-}#pi^{-}#gamma#gamma Invariant Mass",50,1.8,2.3);

  // Go back up to the (ROOT) parent directory
  dir->cd("../");

  pthread_mutex_init(&mutex, NULL);

  return NOERROR;
}

//------------------
// evnt 
//------------------
jerror_t DEventProcessor_b1pi_kinfit_simple::evnt(JEventLoop *loop, int eventnumber)
{

  //Get detected particles from factories
  vector<const DBeamPhoton*> beam_photons;
  vector<const DChargedTrackHypothesis*> track_hypotheses;
  vector<const DNeutralParticleHypothesis*> neutrals;
  loop->Get(beam_photons);
  loop->Get(track_hypotheses);
  loop->Get(neutrals);

  //create containers for each particle type we are interested in
  vector<const DNeutralParticleHypothesis*> photons;
  vector<const DChargedTrackHypothesis*> pi_plus;
  vector<const DChargedTrackHypothesis*> pi_minus;
  vector<const DChargedTrackHypothesis*> protons;

  //only keep neutrals if they are photons
  for (int i=0; i<neutrals.size(); i++) {
    if (neutrals[i]->dPID == Gamma) {
      photons.push_back(neutrals[i]);
    }
  }

  //sort track hypotheses by particle type
  for (int i=0; i<track_hypotheses.size(); i++) {
    Particle_t PID=track_hypotheses[i]->dPID;
    if (PID == PiPlus) {
      pi_plus.push_back(track_hypotheses[i]);
    } else if (PID == PiMinus) {
      pi_minus.push_back(track_hypotheses[i]);
    } else if (PID == Proton) {
      protons.push_back(track_hypotheses[i]);
    }
  }

  // Some generators don't supply information on the beam photon. If there
  // are no DBeamPhoton objects, then create a 9GeV one here
  DBeamPhoton beam;
  if(beam_photons.size()==0) {
    DVector3 mom(0.0, 0.0, 9.0);
    DVector3 pos(0.0, 0.0, 65.0); // center of target
    beam.setMomentum(mom);
    beam.setPosition(pos);
    beam.setMass(0.0);
    beam.setCharge(0.0);
    beam.setMassFixed();
    beam.clearErrorMatrix();
  } else {
    beam = *beam_photons[0]; //assume for now that we only have one beam photon
  }

  int n_photons = photons.size();
  int n_pi_plus = pi_plus.size();
  int n_pi_minus = pi_minus.size();
  int n_protons = protons.size();

  //check if we have sufficient particles to form a p pi+pi+pi-pi-gamma gamma final state
  if (n_photons >= 2 && n_pi_plus >= 2 && n_pi_minus>=2 && n_protons >=1) {
    const double CLCut = .01; //only consider fits with confidence level > 1%
    double maxCL = 0; //highest CL recorded from all combinations
    //save pulls and X mass from only highest CL fit. use these for filling histograms and stuff
    double pulls_highCL[n_final_particles*3];
    double mass_2pip_2pim_2gamma_highCL = 0;

    //loop over all possible of combinations of particles that produce the desired final state
    for (int k=0; k<n_protons; k++) {
      for (int l=0; l<(n_pi_plus-1); l++) {
        //For each positive track, fits are attempted for both proton and pi+
        //mass hypotheses. If the proton and pi+ are from the same track, this
        //is not a legitimate final state, so skip this track
        if (protons[k]->candidateid == pi_plus[l]->candidateid) continue;
        for (int m=l+1; m<n_pi_plus; m++) {
          if (protons[k]->candidateid == pi_plus[m]->candidateid) continue;
          for (int n=0; n<(n_pi_minus-1); n++) {
            for (int o=n+1; o<n_pi_minus; o++) {
              for (int p=0; p<(n_photons-1); p++) {
                for (int q=p+1; q<n_photons; q++) {

                  //vectors for holding the initial and final state particles
                  vector<DKinematicData> initial;
                  vector<DKinematicData> final;

                  //Add the beam particle to the initial state
                  initial.push_back(beam);

                  //Add the target to the initial state
                  DKinematicData kd_targ;
                  kd_targ.setMass(0.93827);
                  kd_targ.setCharge(1.0);
                  kd_targ.setMassFixed();
                  kd_targ.setMomentum( DVector3(0.0, 0.0, 0.0) );
                  kd_targ.clearErrorMatrix();
                  initial.push_back(kd_targ);
	
                  //Add all our particles to the final state
                  final.push_back(*protons[k]);
                  final.push_back(*pi_plus[l]);
                  final.push_back(*pi_plus[m]);
                  final.push_back(*pi_minus[n]);
                  final.push_back(*pi_minus[o]);
                  final.push_back(*photons[p]);
                  final.push_back(*photons[q]);
	
                  //set up the fit
                  //the fit uses only constraints of conservation of momentum and the pi0 mass constraint
                  DKinFit kfit;
                  vector<int> photonNum; //which particles are the photons constrained to the pi0 mass
                  const double pi0_mass=.135;
                  //5 and 6 are the indices (in final) of the two photons to be constrained
                  photonNum.push_back(5);
                  photonNum.push_back(6);
                  kfit.SetMassConstraint(pi0_mass,photonNum);
                  kfit.SetFinal(final);
                  kfit.SetInitial(initial);

                  //do the fit
                  kfit.Fit();

                  //get the corrected, post-fit kinematic data
                  vector<DKinematicData> initial_post;
                  vector<DKinematicData> final_post;
                  final_post = kfit.GetFinal_out();
                  initial_post = kfit.GetInitial_out();

                  //get the confidence level of fit
                  double cl=kfit.Prob();

                  //if this is the highest c.l. fit of all combos in this event
                  //save the pulls and other info
                  if (cl>CLCut && cl > maxCL) {
                    maxCL=cl;                    

                    for (int r=0;r<n_final_particles*3;r++) {
                      //The first 6 pulls correspond to the 2 initial particles
                      //(there are three pulls--px, py, pz--per measured
                      //particle). Since in this example the momenta of
                      //initial state particles is fixed (i.e. covariance
                      //matrix = 0), we only consider the pulls of the final
                      //state particles, hence the 6 in the line below.
                      pulls_highCL[r]=kfit.GetPull(r+6);
                    }

                    //determine the invariant mass of the 5pi system
                    //the loop below starts at 1 to skip the proton
                    DLorentzVector vect_2pip_2pim_2gamma;
                    for (int r=1;r<7;r++) {
                      vect_2pip_2pim_2gamma += final_post[r].lorentzMomentum();
                    }
                    mass_2pip_2pim_2gamma_highCL = vect_2pip_2pim_2gamma.M();
                  }
                }
              }
            }
          }
        }
      }
    }

    //fill histograms with values from highest c.l. fit
    if (maxCL > CLCut) {
      // Fill histograms below here
      pthread_mutex_lock(&mutex);

      h_cl->Fill(maxCL);
      h_mass_2pip_2pim_2gamma->Fill(mass_2pip_2pim_2gamma_highCL);
      for (int q=0;q<n_final_particles*3;q++) {
        h_pulls[q]->Fill(pulls_highCL[q]);
      }

      pthread_mutex_unlock(&mutex);
    }
  }
    
  return NOERROR;

}

DEventProcessor_b1pi_kinfit_simple.h

#ifndef _DEventProcessor_b1pi_kinfit_simple_
#define _DEventProcessor_b1pi_kinfit_simple_

#include <pthread.h>

#include <TH1.h>
#include <TH2.h>

#include <JANA/JEventProcessor.h>
#include <JANA/JEventLoop.h>
using namespace jana;

class DEventProcessor_b1pi_kinfit_simple:public JEventProcessor{

 public:
  DEventProcessor_b1pi_kinfit_simple(){};
  ~DEventProcessor_b1pi_kinfit_simple(){};

  //Total number of final state particles: 1 proton, 2 pi+, 2 pi-, 2 photons
  static const int n_final_particles=7;

  //histograms

  //invariant mass of 5pi system = mass of X resonance. calculated only when
  //the confidence is above a specified threshold.
  TH1F *h_mass_2pip_2pim_2gamma;

  TH1F *h_cl; //highest confidence level of all fits in an event
  TH1F *h_pulls[n_final_particles*3]; //final particles*3 components of momentum

 private:
  jerror_t init(void);	///< Invoked via DEventProcessor virtual method
  jerror_t evnt(JEventLoop *loop, int eventnumber);	///< Invoked via DEventProcessor virtual method

  pthread_mutex_t mutex;
};

#endif // _DEventProcessor_b1pi_kinfit_simple_