minimumTestPHGenFit.cc

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file minimumTestPHGenFit.cc

//STL
#include <vector>

//ROOT
#include <TVector3.h>
#include <TMatrixDSym.h>

//GenFit
#include <GenFit/AbsTrackRep.h>
#include <GenFit/RKTrackRep.h>
#include <GenFit/StateOnPlane.h>

//PHGenFit
#include <phgenfit/Fitter.h>
#include <phgenfit/Track.h>
#include <phgenfit/Measurement.h>
#include <phgenfit/PlanarMeasurement.h>
#include <phgenfit/SpacepointMeasurement.h>

#include <phfield/PHFieldUtility.h>

#define LogDEBUG    std::cout<<"DEBUG: "<<__LINE__<<"\n"

void get_seed(TVector3& seed_pos, TVector3& seed_mom, TMatrixDSym& seed_cov)
{
  seed_pos.SetXYZ(0,0,0);
  seed_mom.SetXYZ(10,-5,0);
  seed_cov.ResizeTo(6,6);
}

std::vector<TVector3> get_raw_measurements()
{
  std::vector<TVector3> v_pos;
  v_pos.push_back(TVector3(2.22459,-1.54767,-2.37792));
  v_pos.push_back(TVector3(3.80050,-2.64444,-2.16561));
  v_pos.push_back(TVector3(7.80344,-5.41815,-1.98928));
  v_pos.push_back(TVector3(8.63214,-5.97797,-1.59626));
  return v_pos;
}

int main(int argc, char**argv) {

  PHGenFit::Fitter* fitter = new PHGenFit::Fitter("sPHENIX_Geo.root",
      PHFieldUtility::BuildFieldMap(PHFieldUtility::DefaultFieldConfig(), 1),
      "KalmanFitter","RKTrackRep",false);

  int pid = -13; //mu+
  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(pid);

  TVector3 seed_pos;
  TVector3 seed_mom;
  TMatrixDSym seed_cov;
  get_seed(seed_pos,seed_mom, seed_cov);
  PHGenFit::Track* track = new PHGenFit::Track(rep, seed_pos,seed_mom, seed_cov);

  std::vector<TVector3> v_pos = get_raw_measurements();
  double res_phi = 0.005; //cm
  double res_z = 0.04; //cm
  std::vector<PHGenFit::Measurement*> measurements;
  for (unsigned int imeasurement = 0; imeasurement < v_pos.size(); imeasurement++) {
    TVector3 pos = v_pos[imeasurement];
    TVector3 n(pos.x(),pos.Y(),0);
    PHGenFit::Measurement* meas = new PHGenFit::PlanarMeasurement(pos,n,res_phi, res_z);
    //PHGenFit::Measurement* meas = new PHGenFit::SpacepointMeasurement(pos,res_phi);
    meas->getMeasurement()->Print();
    measurements.push_back(meas);
  }

  track->addMeasurements(measurements);

  fitter->processTrack(track, false);

  //genfit::MeasuredStateOnPlane* state = track->extrapolateToLine(TVector3(0, 0, 0), TVector3(0, 0, 1));
  //genfit::MeasuredStateOnPlane* state = track->extrapolateToPoint(TVector3(0, 0, 0));
  genfit::MeasuredStateOnPlane* state = track->extrapolateToCylinder(1.,TVector3(0, 0, 0), TVector3(0, 0, 1));
  state->Print();
  delete state;

  //fitter->displayEvent();


  delete track;

  delete fitter;

  return 0;
}