PHGenFitPkg/Example/minimumTestPHGenFit.cc

0001 /*!
0002  *  \file       minimumTestPHGenFit.cc
0003  *  \brief      Minimum program to demonstrate the usage of PHGenFit.
0004  *  \author     Haiwang Yu <yuhw@nmsu.edu>
0005  */
0006
0007 // STL
0008 #include <vector>
0009
0010 // ROOT
0011 #include <TMatrixDSym.h>
0012 #include <TVector3.h>
0013
0014 // GenFit
0015 #include <GenFit/AbsTrackRep.h>
0016 #include <GenFit/RKTrackRep.h>
0017 #include <GenFit/StateOnPlane.h>
0018
0019 // PHGenFit
0020 #include <phgenfit/Fitter.h>
0021 #include <phgenfit/Measurement.h>
0022 #include <phgenfit/PlanarMeasurement.h>
0023 #include <phgenfit/SpacepointMeasurement.h>
0024 #include <phgenfit/Track.h>
0025
0026 #include <phfield/PHFieldUtility.h>
0027
0028 #define LogDEBUG std::cout << "DEBUG: " << __LINE__ << "\n"
0029
0030 void get_seed(TVector3& seed_pos, TVector3& seed_mom, TMatrixDSym& seed_cov)
0031 {
0032   seed_pos.SetXYZ(0, 0, 0);
0033   seed_mom.SetXYZ(10, -5, 0);
0034   seed_cov.ResizeTo(6, 6);
0035 }
0036
0037 std::vector<TVector3> get_raw_measurements()
0038 {
0039   std::vector<TVector3> v_pos;
0040   v_pos.emplace_back(2.22459, -1.54767, -2.37792);
0041   v_pos.emplace_back(3.80050, -2.64444, -2.16561);
0042   v_pos.emplace_back(7.80344, -5.41815, -1.98928);
0043   v_pos.emplace_back(8.63214, -5.97797, -1.59626);
0044   return v_pos;
0045 }
0046
0047 int main(int /*argc*/, char** /*argv*/)
0048 {
0049   //! Initiallize Geometry, Field, Fitter
0050   PHGenFit::Fitter* fitter = new PHGenFit::Fitter("sPHENIX_Geo.root",
0051                                                   PHFieldUtility::BuildFieldMap(PHFieldUtility::DefaultFieldConfig(), 1),
0052                                                   "KalmanFitter", "RKTrackRep", false);
0053
0054   //! Build TrackRep from particle assumption
0055   int pid = -13;  // mu+
0056   genfit::AbsTrackRep* rep = new genfit::RKTrackRep(pid);
0057
0058   //! Initiallize track with seed from pattern recognition
0059   TVector3 seed_pos;
0060   TVector3 seed_mom;
0061   TMatrixDSym seed_cov;
0062   get_seed(seed_pos, seed_mom, seed_cov);
0063   PHGenFit::Track* track = new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov);
0064
0065   //! Create measurements
0066   std::vector<TVector3> v_pos = get_raw_measurements();
0067   double res_phi = 0.005;  // cm
0068   double res_z = 0.04;     // cm
0069   std::vector<PHGenFit::Measurement*> measurements;
0070   for (const auto& pos : v_pos)
0071   {
0072     TVector3 n(pos.x(), pos.Y(), 0);
0073     PHGenFit::Measurement* meas = new PHGenFit::PlanarMeasurement(pos, n, res_phi, res_z);
0074     // PHGenFit::Measurement* meas = new PHGenFit::SpacepointMeasurement(pos,res_phi);
0075     meas->getMeasurement()->Print();
0076     measurements.push_back(meas);
0077   }
0078
0079   //! Add measurements to track
0080   track->addMeasurements(measurements);
0081
0082   //! Fit the track
0083   fitter->processTrack(track, false);
0084
0085   //! Extrapolate to beam line
0086   // genfit::MeasuredStateOnPlane* state = track->extrapolateToLine(TVector3(0, 0, 0), TVector3(0, 0, 1));
0087   // genfit::MeasuredStateOnPlane* state = track->extrapolateToPoint(TVector3(0, 0, 0));
0088   genfit::MeasuredStateOnPlane* state = track->extrapolateToCylinder(1., TVector3(0, 0, 0), TVector3(0, 0, 1));
0089   state->Print();
0090   delete state;
0091
0092   //! Event display, uncomment to use
0093   // fitter->displayEvent();
0094
0095   //! Comment off if want to keep event display.
0096
0097   delete track;
0098
0099   delete fitter;
0100
0101   return 0;
0102 }