HOWTO do a kinematic fit for etapi0p events

Introduction

This bit of code was taken from presentations by Matt Bellis which can be found on the DocDB and included in code I used to reconstruct ηπ⁰ events. I will try to make it as complete as possible where I can. Since there are previous examples of how to generate and simulate events I will not include this, but you can look in HOWTO simulate and analyze b1pi events for the steps. Everything that follows here assumes that the events have been passed through HDGEANT and mcsmear.

Reconstruction

The code that follows is built on the JANA framework and is compiled with a standard Makefile which I will include at the bottom of this page. I've cut and pasted this from a working program so I have not checked if it compiles.

FCALTreeProc.cc

#include "FCALTreeProc.h"
#include "JANA/JApplication.h"
#include "DANA/DApplication.h"
#include "PID/DKinematicData.h"
#include "TRACKING/DMCThrown.h"
#include "TRACKING/DMagneticFieldStepper.h"
#include "HDGEOMETRY/DMagneticFieldMap.h"
#include "FCAL/DFCALGeometry.h"
#include "FCAL/DFCALTruthShower.h"
#include "FCAL/DFCALHit.h"
#include "FCAL/DFCALCluster.h"
#include "FCAL/DFCALPhoton.h"
#include "BCAL/DBCALTruthShower.h"
#include "BCAL/DBCALShower.h"
#include "BCAL/DBCALPhoton.h" 
#include "PID/DPhoton.h"
#include "PID/DPi0.h"
#include "PID/DTwoGammaFit.h"
#include "PID/DParticle.h" 
#include "PID/DKinFit.h"
#include "PID/DTwoGammaFit_factory.h"
#include "TBranch.h"
#include "TH1F.h"
#include "TH2F.h"
#include <DLorentzVector.h>
#include <TLorentzRotation.h>
#ifndef PI
 #define PI 3.1415927
#endif

Now to we'll define a ROOT tree and a few branches to save some data to. You can add any branch you feel like, just make sure to add the appropriate lines in the header file.

//------------------
// init
//------------------
jerror_t FCALTreeProc::init(void)
{
      // Create histograms here
       m_treeFile = new TFile( "fcaltree4.root", "RECREATE" ); //the ROOT file everything is written out to
       m_tree = new TTree("Fcal","FCAL Tree");  // a tree in the ROOT file we'll add branches too
      //     m_outFile = new TFile( outFile.c_str(), "recreate" );
      m_outTree = new TTree( "kin", "Kinematics" ); // a second tree with the branches/leaves for IU PWA reconstruction
      m_tree->Branch( "ProbMax",&ProbMax,"ProbMax/F");
      m_tree->Branch( "ProbMax_PI0",&ProbMax_PI0,"ProbMax_PI0/F");
      m_tree->Branch( "ProbMax_ETA",&ProbMax_ETA,"ProbMax_ETA/F");
      m_tree->Branch( "ProbMax_FIT",&ProbMax_FIT,"ProbMax_FIT/F");
      m_tree->Branch( "PI0CHI2",&PI0CHI2,"PI0CHI2/F");
      m_tree->Branch( "ETACHI2",&ETACHI2,"ETACHI2/F");
      m_tree->Branch( "FITCHI2",&FITCHI2,"FITCHI2/F");
      m_tree->Branch( "FITNDF",&FITNDF,"FITNDF/F");
      m_tree->Branch( "ETAMASS",&ETAMASS,"ETAMASS/F");
      m_tree->Branch( "PI0MASS",&PI0MASS,"PI0MASS/F");
      m_tree->Branch( "RESMASS",&RESMASS,"RESMASS/F");
      m_tree->Branch( "RESENERGY",&RESENERGY,"RESENERGY/F");
      m_tree->Branch( "RESMOM",&RESMOM,"RESMOM/F");

etc... I have a few hundred or so more which I'm not going to include here, but you can easily add your own, depending on what information you want to record. Here's the branches in the IU PWA kin tree.

      m_outTree->Branch( "nPart", &nPart, "nPart/I" );
      m_outTree->Branch( "e", m_e, "e[nPart]/F" );
      m_outTree->Branch( "px", m_px, "px[nPart]/F" );
      m_outTree->Branch( "py", m_py, "py[nPart]/F" );
      m_outTree->Branch( "pz", m_pz, "pz[nPart]/F" );
      m_outTree->Branch( "eBeam", &m_eBeam, "eBeam/F" );
      m_outTree->Branch( "pxBeam", &m_pxBeam, "pxBeam/F" );
      m_outTree->Branch( "pyBeam", &m_pyBeam, "pyBeam/F" );
      m_outTree->Branch( "pzBeam", &m_pzBeam, "pzBeam/F" );
      m_outTree->Branch( "eRecoil", &m_eRecoil, "eRecoil/F" );
      m_outTree->Branch( "pxRecoil", &m_pxRecoil, "pxRecoil/F" );
      m_outTree->Branch( "pyRecoil", &m_pyRecoil, "pyRecoil/F" );
      m_outTree->Branch( "pzRecoil", &m_pzRecoil, "pzRecoil/F" );
      return NOERROR;
}

Here's the juicy bits.

      /////////////////////////////////////////
      //////begin reconstructing HDGEANT events
      /////////////////////////////////////////
//
      vector<const DPhoton*> photons;   //we get anything that looks like a DPhoton and put it in a vector
      loop->Get(photons);
      vector<const DParticle*> recoilproton; //we get anything that looks like a DParticle and put it in a vector
      loop->Get(recoilproton);
//
      //photon multiplicity
      nPhotons = photons.size();
      nPhotons_ps = 0; 
      for(unsigned int i = 0; i < photons.size() ; i++) {
       if ( photons[i]->getTag() ==  DPhoton::kCharge  ) continue;

nPhotons_ps++;

      } //this will give us the the number of DPhotons after suppressing clusters due to charged particles.
//  
//now we define some Lorentz Vectors we'll need for the kinematic fit
TLorentzVector P4_beam( 0.0, 0.0, 9.0, 9.0 ); // the 9 GeV photon beam in the z-direction
TLorentzVector P4_target( 0, 0, 0, 0.938 ); // the stationary proton target
TLorentzVector P4_RECOIL(0,0,0,0); //initializes the recoil TLorentzVector
//
 if( recoilproton.size() == 1 || recoilproton.size() == 0) { //we only want one charged particle, or we missed it
    if ( nPhotons_ps >= 4 &&  nPhotons_ps <= 8){ //a reasonable constraint on the number of showers we might see
// now loop over all possible unique combinations of 4 photons, excluding charged clusters
for(unsigned int i = 0; i < photons.size() ; i++) {
         if ( photons[i]->getTag() == DPhoton::kCharge  ) continue;
           for (unsigned int j = i+1; j < photons.size() ; j++) {
               if ( photons[j]->getTag() == DPhoton::kCharge  ) continue;
                  for (unsigned int k = 0; k < photons.size() ; k++) {
                      if ( photons[k]->getTag() == DPhoton::kCharge ||  k == i || k == j ) continue;

for (unsigned int l = k+1; l < photons.size() ; l++) {

                              if ( photons[l]->getTag() == DPhoton::kCharge || l == i || l == j )continue;