Difference between revisions of "HOWTO do a kinematic fit for etapi0p events"

From GlueXWiki
Jump to: navigation, search
m
Line 90: Line 90:
 
       //photon multiplicity
 
       //photon multiplicity
 
       nPhotons = photons.size();
 
       nPhotons = photons.size();
      nPhotons_ps = 0;  
+
        nPhotons_ps = 0;  
 
       for(unsigned int i = 0; i < photons.size() ; i++) {
 
       for(unsigned int i = 0; i < photons.size() ; i++) {
 
         if ( photons[i]->getTag() ==  DPhoton::kCharge  ) continue;
 
         if ( photons[i]->getTag() ==  DPhoton::kCharge  ) continue;
  nPhotons_ps++;
+
        nPhotons_ps++;
 
       } //this will give us the the number of DPhotons after suppressing clusters due to charged particles.
 
       } //this will give us the the number of DPhotons after suppressing clusters due to charged particles.
 
  //   
 
  //   
Line 173: Line 173:
 
     constraintParticles.push_back(2);
 
     constraintParticles.push_back(2);
 
     constraintParticles.push_back(3);
 
     constraintParticles.push_back(3);
     kfit->SetMassConstraint(0.54775,constraintParticles); //constrains the second two particles to the eta mass
+
     kfit->SetMassConstraint(0.54775,constraintParticles); //constrains the second two photons to the eta mass
 
now to get the output from the fitter
 
now to get the output from the fitter
 
     kfit->Fit();  // the fit is actually done here
 
     kfit->Fit();  // the fit is actually done here
Line 213: Line 213:
 
       RECPT = P4_RECOIL.Pt();
 
       RECPT = P4_RECOIL.Pt();
 
       t = -( P4_RECOIL - P4_target ).Mag2();
 
       t = -( P4_RECOIL - P4_target ).Mag2();
 +
 +
now to select count and select the best candidate
 +
if (FitChi2[nPairs]/FitNDF[nPairs] < 10) nFitCandChi2_50++; //counts the candidates with Chi2 < 10
 +
if (FitProb[nPairs] > 0.01) {
 +
    nFitCandProb++; //counts the candidates with probability > 1%
 +
    write_out = true;
 +
    }
 +
//if the combination has the maximum probability choose this combination as the candidate
 +
    if (FitProb[nPairs] > ProbMax){
 +
      ProbMax = FitProb[nPairs];
 +
      ProbMax_FIT = FitProb[nPairs];
 +
      FITCHI2 = FitChi2[nPairs];
 +
      FITNDF = FitNDF[nPairs];
 +
      for(unsigned int pullk = 0; pullk < 15; pullk++){
 +
          PULL[pullk] = FitPull[pullk][nPairs];
 +
      }
 +
//
 +
    PI0MASS = P4_PI0.M();
 +
    ETAMASS = P4_ETA.M();
 +
    P4_RES = P4_PI0 + P4_ETA; //resonance 4-momenta
 +
    RESMASS = P4_RES.M();
 +
    RESENERGY = P4_RES.T();
 +
    RESMOM =  P4_RES.Vect().Mag();
 +
    RESPT = P4_RES.Pt();
 +
 +
now to calculate the cosine of the Gottfried-Jackson angle
 +
    if ( (P4_RES.Beta() > 0 && P4_RES.Beta() < 1) && P4_ETA.T() < 10 && P4_PI0.T() < 10 ) { //sanity check
 +
    TLorentzRotation resRestBoost( -P4_RES.BoostVector() );
 +
    DLorentzVector beam_res  = resRestBoost * P4_beam;
 +
    DLorentzVector recoil_res = resRestBoost * P4_RECOIL;
 +
    DLorentzVector p1_res = resRestBoost * P4_PI0;
 +
//
 +
    DVector3 z = beam_res.Vect().Unit();
 +
    DVector3 y = recoil_res.Vect().Cross(z).Unit();
 +
    DVector3 x = y.Cross(z);
 +
    DVector3 angles( (p1_res.Vect()).Dot(x),(p1_res.Vect()).Dot(y),(p1_res.Vect()).Dot(z));
 +
    cosGJ = angles.CosTheta(); //cosine of the Gottfried-Jackson angle
 +
    phiGJ = angles.Phi(); //azimuthal Gottfried-Jackson angle
 +
    }
 +
        }
 +
  //
 +
    nPairs++; //useful iterator
 +
                    }  //end of "for (unsigned int l...){" 
 +
                }  //end of "for (unsigned int j...){" 
 +
            }  //end of "for (unsigned int j...){" 
 +
          } //end of "for (unsigned int i...){" 
 +
        } //end of "if (nPhotons_ps >= 4...){"
 +
    }  //end of "if( recoilproton.size() == 1 && ...){"
 +
 +
if the write_out condition is true:
 +
if( write_out ){
 +
    flush_s_HDDM(hddm, file);
 +
    Nevents_written++;
 +
  }
 +
hddm_source->flush_on_free = !write_out;
 +
        m_outTree->Fill();
 +
        m_tree->Fill();
 +
        return NOERROR;
 +
}
 +
//------------------
 +
// erun
 +
//------------------
 +
jerror_t FCALTreeProc::erun(void)
 +
{
 +
// Any final calculations on histograms (like dividing them)
 +
// should be done here. This may get called more than once.
 +
return NOERROR;
 +
}
 +
//------------------
 +
// fini
 +
//------------------
 +
jerror_t FCALTreeProc::fini(void)
 +
{
 +
        if(file){
 +
        close_s_HDDM(file);
 +
        cout<<endl<<"Closed HDDM file"<<endl;
 +
        }
 +
        cout<<" "<<Nevents_written<<" events written to "<<filename<<endl;
 +
        // If another JEventProcessor is in the list ahead of this one, then
 +
        // it will have finished before this is called. e.g. closed the
 +
        // ROOT file!
 +
        m_treeFile->Write();
 +
        m_treeFile->Close();
 +
        cout<<endl<<"Closed ROOT file"<<endl;
 +
        return NOERROR;
 +
}

Revision as of 17:01, 27 May 2010

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 ηπ0 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

All the headers files you might need:

#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;

so while we're looping over all the combinations of 4 photons, lets do a kinematic fit

      ////////////////////////////////
      ///Fit Initial and Final States
      ////////////////////////////////
//   
   //dumerr(0,0) = 1.0; //dummy error matrix, not sure it's needed anymore
//
   vector<DKinematicData> kd_initialState; //initial state kinematics
   DKinematicData kd_beam = DKinematicData();//beam kinematics
   kd_beam.setMass(0.0);
   kd_beam.setCharge(0);
   kd_beam.setMassFixed();
   kd_beam.setMomentum( TVector3(0.0, 0.0, 9.0));
   kd_beam.clearErrorMatrix(); // set error matrix with entries equal to zero
   kd_initialState.push_back(kd_beam);
   kd_initialState[0].smearMCThrownMomentum(0.001); //not sure if this actually smears the beams by 0.1%
//
   DKinematicData kd_target = DKinematicData();  //target kinematics
   kd_target.setMass(0.93827);
   kd_target.setCharge(1);
   kd_target.setMassFixed();
   kd_target.setMomentum( TVector3(0.0, 0.0, 0.0));
   kd_target.clearErrorMatrix();  // set error matrix with entries equal to zero
   kd_initialState.push_back(kd_target);
   kd_initialState[1].smearMCThrownMomentum(0.00001);
//
   vector<DKinematicData> kd_finalState; //final state kinematics
//
 //build the 5 particles in the fit
   kd_finalState.push_back(*(photons[i]));
   kd_finalState.push_back(*(photons[j]));
   kd_finalState.push_back(*(photons[k]));
   kd_finalState.push_back(*(photons[l]));
   if( recoilproton.size() == 1  ) {
       P4_RECOIL = recoilproton[0]->lorentzMomentum();
       kd_finalState.push_back(*(recoilproton[0]));
       // kd_finalState[4].clearErrorMatrix(); //if you don't have errors in HDParSim, you will have problems
   }
//	
// now put the kinematic data into the fitter		 
   DKinFit *kfit = new DKinFit(); 
   kfit->ResetForNewFit();
   kfit->SetVerbose(0);  // 0 = quiet,1 = some info, 2 = all info
   kfit->UseTimingInformation(); // 
   kfit->SetInitial(kd_initialState);
   kfit->SetFinal(kd_finalState);
   if(recoilproton.size() == 0) { //the acceptance of the tracking system misses protons with momentum less than 450 MeV
    //   cout << "NPart= = " <<recoilproton.size() << endl;
   kfit->SetMissingParticle(0.938);
   }
//

we want to constrain the invariant mass of two of the photons to the pion mass and two to the eta mass

   vector<int> constraintParticles;
   constraintParticles.push_back(0);
   constraintParticles.push_back(1);
   kfit->SetMassConstraint(0.13498,constraintParticles); //constrains first two photons to the pi0 mass
   constraintParticles.clear();
   constraintParticles.push_back(2);
   constraintParticles.push_back(3);
   kfit->SetMassConstraint(0.54775,constraintParticles); //constrains the second two photons to the eta mass

now to get the output from the fitter

   kfit->Fit();  // the fit is actually done here
   //get the fit chi2 and probabily and pulls
   FitProb[nPairs] = kfit->Prob();
   FitChi2[nPairs] = kfit->Chi2();
   FitNDF[nPairs] = kfit->Ndf();
//I'm only interested in the pulls of the photons
   for(unsigned int pullk = 0; pullk<15;pullk++){
   FitPull[pullk][nPairs] = kfit->GetPull(pullk+6);
   }
  //CovMat_FIT = kfit->GetCovMat();
   //fitter has adjusted the 4-momenta of the particles within error. get these new values.
   vector<DKinematicData> kinout = kfit->GetFinal_out();
   vector<DKinematicData> kinin = kfit->GetFinal_in();
   delete kfit;  //cleanup

we'll make pi0s and eta out of these photons now by adding the 4-momenta

   TLorentzVector P4_PI0,P4_ETA,P4_REC,P4_RES,P4_MISS4M;		
    P4_PI0 = kinout[0].lorentzMomentum() + kinout[1].lorentzMomentum(); 
    P4_ETA = kinout[2].lorentzMomentum() + kinout[3].lorentzMomentum();
// 
   if( recoilproton.size() == 1  ) {
         P4_RECOIL = kinout[4].lorentzMomentum();
         P4_MISS4M = P4_beam + P4_target - P4_ETA - P4_PI0 -P4_RECOIL;
        }
        else {
          P4_RECOIL = P4_beam + P4_target - P4_ETA - P4_PI0;
          P4_MISS4M = P4_RECOIL;
        }
      MISMASS = P4_MISS4M.M();
      MISMASS2 = P4_MISS4M.M2();
      MISENERGY = P4_MISS4M.T();
      MISMOM = P4_MISS4M.Vect().Mag();
      MISPT = P4_MISS4M.Pt();
      RECMASS = P4_RECOIL.M();
      RECENERGY = P4_RECOIL.T();
      RECMOM = P4_RECOIL.Vect().Mag();
      RECPT = P4_RECOIL.Pt();
      t = -( P4_RECOIL - P4_target ).Mag2();

now to select count and select the best candidate

if (FitChi2[nPairs]/FitNDF[nPairs] < 10) nFitCandChi2_50++; //counts the candidates with Chi2 < 10
if (FitProb[nPairs] > 0.01) {
    nFitCandProb++; //counts the candidates with probability > 1%
    write_out = true;
    }
//if the combination has the maximum probability choose this combination as the candidate
   if (FitProb[nPairs] > ProbMax){
     ProbMax = FitProb[nPairs];
     ProbMax_FIT = FitProb[nPairs];
     FITCHI2 = FitChi2[nPairs];
     FITNDF = FitNDF[nPairs];
     for(unsigned int pullk = 0; pullk < 15; pullk++){
          PULL[pullk] = FitPull[pullk][nPairs];
      }
//
    PI0MASS = P4_PI0.M();
    ETAMASS = P4_ETA.M();
    P4_RES = P4_PI0 + P4_ETA; //resonance 4-momenta
    RESMASS = P4_RES.M();
    RESENERGY = P4_RES.T();
    RESMOM =  P4_RES.Vect().Mag();
    RESPT = P4_RES.Pt();

now to calculate the cosine of the Gottfried-Jackson angle

   if ( (P4_RES.Beta() > 0 && P4_RES.Beta() < 1) && P4_ETA.T() < 10 && P4_PI0.T() < 10 ) { //sanity check
   TLorentzRotation resRestBoost( -P4_RES.BoostVector() );
   DLorentzVector beam_res   = resRestBoost * P4_beam;
   DLorentzVector recoil_res = resRestBoost * P4_RECOIL;
   DLorentzVector p1_res = resRestBoost * P4_PI0;
//
   DVector3 z = beam_res.Vect().Unit();
   DVector3 y = recoil_res.Vect().Cross(z).Unit();
   DVector3 x = y.Cross(z);
   DVector3 angles( (p1_res.Vect()).Dot(x),(p1_res.Vect()).Dot(y),(p1_res.Vect()).Dot(z));
   cosGJ = angles.CosTheta(); //cosine of the Gottfried-Jackson angle
   phiGJ = angles.Phi(); //azimuthal Gottfried-Jackson angle
   }
       }
 //
   nPairs++; //useful iterator
                   }   //end of "for (unsigned int l...){"  
                }  //end of "for (unsigned int j...){"  
            }  //end of "for (unsigned int j...){"  
         } //end of "for (unsigned int i...){"  
       } //end of "if (nPhotons_ps >= 4...){"
    }  //end of "if( recoilproton.size() == 1 && ...){"

if the write_out condition is true:

if( write_out ){
   flush_s_HDDM(hddm, file);
   Nevents_written++;
 }
hddm_source->flush_on_free = !write_out;
       m_outTree->Fill();
       m_tree->Fill();
       return NOERROR;
}
//------------------
// erun
//------------------
jerror_t FCALTreeProc::erun(void)
{ 
// Any final calculations on histograms (like dividing them)
// should be done here. This may get called more than once.
return NOERROR;
}
//------------------
// fini
//------------------
jerror_t FCALTreeProc::fini(void)
{
       if(file){
        close_s_HDDM(file);
        cout<<endl<<"Closed HDDM file"<<endl;
       }
       cout<<" "<<Nevents_written<<" events written to "<<filename<<endl;
       // If another JEventProcessor is in the list ahead of this one, then
       // it will have finished before this is called. e.g. closed the
       // ROOT file!
       m_treeFile->Write();
       m_treeFile->Close();
       cout<<endl<<"Closed ROOT file"<<endl;
       return NOERROR;
}