HOWTO do a kinematic fit for etapi0p events
by Blake Leverington (UofR)
Introduction
This bit of code was taken from presentations by Matt Bellis which can be found on the DocDB , and code by Matt Shepherd, Mihajlo Kornicer, David Lawrence and probably others, which I included in code that I used to reconstruct ηπ0 events. You can find some of this under Kinematic Fitting. 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. HDParSim was also used to handle the charged particles.
Note: this will only work properly if HDParSim has error matrices for the particles. The error matrices were not part of the repository version at the time I was working on this and was provided by Will Levine (CMU?).
In the end, what you should get out is reconstructed 4-momenta for a resonance with the kinematic fit probabilities and the adjusted photons to go along with it, though you also retain the unadjusted photon momenta information if desired.
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. Hopefully it should make more clear than when I started what is required to reconstruct eta-pi0 events. I've cut and pasted this from pieces of a working program (link at the bottom of the page) so it will not compile as it is, but should contain all the important pieces.
FCALTreeProc.cc: the good parts
All the headers files you might need to include (maybe some extra):
#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 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 too.
//------------------ // 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 couple 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; }
some extra stuff to open the data file that is to be analyzed
//------------------ // brun //------------------ jerror_t FCALTreeProc::brun(JEventLoop *eventLoop, int runnumber) { return NOERROR; } //------------------ // evnt //------------------ jerror_t FCALTreeProc::evnt(JEventLoop *loop, int eventnumber) { // Fill histograms here int nHits = 0; float Etot=0; write_out = true; // This is a little complicated. We need to get a hold of the s_HDDM_t // structure pointer for this event so we can pass it to flush_s_HDDM() // along with our ouput stream pointer. The flush routine frees up the // memory in the s_HDDM_t structure. When the framework tries "flush"ing // a second time, we get a seg fault. To prevent the framework from // flushing, we have to clear the free_on_flush flag (by default set // to true). This means we need to get the DEventSource pointer and // downcast to a DEventSourceHDDM structure. It a little strange setting // this for every event, but we have no way of knowing when the event // source changes and this at least guarantees it for all event sources. JEvent& event = loop->GetJEvent(); JEventSource *source = event.GetJEventSource(); DEventSourceHDDM *hddm_source = dynamic_cast<DEventSourceHDDM*>(source); if(!hddm_source){ cerr<<" This program MUST be used with an HDDM file as input!"<<endl; exit(-1); } s_HDDM_t *hddm = (s_HDDM_t*)event.GetRef(); if(!hddm)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. // TLorentzVector P4_RECOIL(0,0,0,0); //initializes the recoil TLorentzVector //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 // if( recoilproton.size() == 1 || recoilproton.size() == 0) { //we only want one charged particle, or we maybe missed it if ( nPhotons_ps >= 4 && nPhotons_ps <= 8){ //a reasonable constraint on the number of showers we might have // // 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 to check if they make a good ηπ0:
//////////////////////////////// ///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();
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; //maybe you only want to write out events with some good probability } //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(); t = -( P4_beam - P4_RES ).Mag2();
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, write the events to file:
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; }
A tarball of the full code with header and Makefile can be found here: fcalTree4.tgz It is not as well documented as this wiki page, and includes a lot of redundant or useless pieces of code. It was a work in progress for my thesis work.
To run this code on a simulation data set you would do something like:
home/s4/leverin/gluex/svn/bin/Linux_CentOS4-athlon-gcc3.4.4/fcalTree4 --plugin=hdparsim -PDEFTAG:DParticle=HDParSim \ -PDEFTAG:DTrackTimeBased=HDParSim --plugin=mcthrown_hists -PTHREAD_TIMEOUT= 60 -PPRINT_PLUGIN_PATHS=1 \ hdgeant_mcsmeared.hddm
This uses HDParSim plugin to handle resolution and acceptance for the charged particles and for suppressing clusters due to charged particles. This speeds up the reconstruction by a factor of 40 or so if you are not interested in fully reconstructing charged particles.