10 #include <phgenfit/Fitter.h>
11 #include <phgenfit/Measurement.h>
12 #include <phgenfit/PlanarMeasurement.h>
13 #include <phgenfit/SpacepointMeasurement.h>
14 #include <phgenfit/Track.h>
15 #include <phgeom/PHGeomUtility.h>
28 #include <calobase/RawTowerGeom.h>
29 #include <calobase/RawTowerGeomContainer.h>
31 #include <phparameter/PHParameters.h>
33 #include <phfield/PHFieldUtility.h>
53 #include <GenFit/AbsMeasurement.h>
54 #include <GenFit/EventDisplay.h>
55 #include <GenFit/MeasuredStateOnPlane.h>
56 #include <GenFit/RKTrackRep.h>
58 #include <GenFit/FitStatus.h>
59 #include <GenFit/GFRaveTrackParameters.h>
60 #include <GenFit/GFRaveVertex.h>
61 #include <GenFit/GFRaveVertexFactory.h>
64 #include <TMatrixDSymfwd.h>
65 #include <TMatrixTSym.h>
66 #include <TMatrixTUtils.h>
69 #include <TVectorDfwd.h>
72 #include <gsl/gsl_randist.h>
73 #include <gsl/gsl_rng.h>
89 #define LogDebug(exp) \
90 if (Verbosity()) std::cout << "PHG4TrackFastSim (DEBUG): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
91 #define LogError(exp) \
92 std::cout << "PHG4TrackFastSim (ERROR): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
93 #define LogWarning(exp) \
94 std::cout << "PHG4TrackFastSim (WARNING): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
107 , m_RaveVertexFactory(nullptr)
108 , m_TruthContainer(nullptr)
109 , m_SvtxTrackMapOut(nullptr)
110 , m_SvtxVertexMap(nullptr)
111 , m_SubTopnodeName(
"SVTX")
112 , m_TrackmapOutNodeName(
"SvtxTrackMap")
113 , m_VertexingMethod(
"kalman-smoothing:1")
114 , m_FitAlgoName(
"DafRef")
115 , m_VertexMinNdf(10.)
116 , m_VertexXYResolution(50
E-4)
117 , m_VertexZResolution(50
E-4)
119 , m_PrimaryAssumptionPid(211)
120 , m_SmearingFlag(
true)
121 , m_DoEvtDisplayFlag(
false)
122 , m_UseVertexInFittingFlag(
true)
123 , m_PrimaryTrackingFlag(1)
124 , m_DoVertexingFlag(
false)
179 switch (iter->second.first)
181 case DETECTOR_TYPE::Cylinder:
183 string nodename =
"TOWERGEOM_" + iter->first;
191 case DETECTOR_TYPE::Vertical_Plane:
193 string towergeonodename =
"TOWERGEOM_" + iter->first;
194 RawTowerGeomContainer* towergeo = findNode::getClass<RawTowerGeomContainer>(topNode, towergeonodename);
197 cout <<
PHWHERE <<
" ERROR: Can't find node " << towergeonodename << endl;
212 cout <<
"invalid state reference: " << iter->second.first << endl;
228 cout <<
PHWHERE <<
" no Vertex Finder" << endl;
250 std::cout <<
"PHG4TrackFastSim::process_event: " <<
m_EventCnt <<
".\n";
261 LogError(
"m_SvtxTrackMapOut not found!");
265 vector<PHGenFit::Track*> rf_tracks;
268 TVector3 vtxPoint(truthVtx->
get_x(), truthVtx->
get_y(), truthVtx->
get_z());
290 itr != itr_range.second; ++itr)
294 TVector3 seed_pos(vtxPoint.x(), vtxPoint.y(), vtxPoint.z());
295 TVector3 seed_mom(0, 0, 0);
296 TMatrixDSym seed_cov(6);
299 std::vector<PHGenFit::Measurement*> measurements;
310 measurements.push_back(vtx_meas);
316 svtx_track_out.get(),
320 if (measurements.size() < 3)
325 std::cout <<
"event: " << m_EventCnt <<
" : measurements.size() < 3"
330 for (
unsigned int im = 0; im < measurements.size(); im++)
332 delete measurements[im]->getMeasurement();
333 delete measurements[im];
359 rf_tracks.push_back(track);
369 if (fitting_err != 0)
374 std::cout <<
"event: " << m_EventCnt
375 <<
" : fitting_err != 0, next track."
381 TVector3 vtx(vtxPoint.x(), vtxPoint.y(), vtxPoint.z());
384 measurements.size(), vtx);
394 const unsigned int track_id = m_SvtxTrackMapOut->insert(svtx_track_out.get())->get_id();
405 cout << __PRETTY_FUNCTION__ <<
"Failed to init vertex finder" << endl;
410 cout << __PRETTY_FUNCTION__ <<
"Failed to init vertex map" << endl;
418 vector<genfit::GFRaveVertex*> rave_vertices;
419 if (rf_tracks.size() >= 2)
423 vector<genfit::Track*> rf_gf_tracks;
424 for (std::vector<PHGenFit::Track*>::iterator
it = rf_tracks.begin();
it != rf_tracks.end(); ++
it)
426 genfit::Track*
track = (*it)->getGenFitTrack();
433 track->getFittedState().getPosMomCov(pos, mom, cov);
435 cout <<
"Track getCharge = " << track->getFitStatus()->getCharge() <<
" getChi2 = " << track->getFitStatus()->getChi2() <<
" getNdf = " << track->getFitStatus()->getNdf() << endl;
441 rf_gf_tracks.push_back(track);
448 std::cout <<
PHWHERE <<
"GFRaveVertexFactory::findVertices failed!";
454 cout << __PRETTY_FUNCTION__ << __LINE__ <<
" rf_tracks = " << rf_tracks.size() <<
" rave_vertices = " << rave_vertices.size() << endl;
462 vector<genfit::Track*> rf_gf_tracks;
463 for (std::vector<PHGenFit::Track*>::iterator
it = rf_tracks.begin();
it != rf_tracks.end(); ++
it)
465 rf_gf_tracks.push_back((*it)->getGenFitTrack());
471 for (std::vector<PHGenFit::Track*>::iterator
it = rf_tracks.begin();
it != rf_tracks.end(); ++
it)
491 const std::vector<genfit::GFRaveVertex*>& rave_vertices,
494 for (genfit::GFRaveVertex* rave_vtx : rave_vertices)
504 svtx_vtx->set_chisq(rave_vtx->getChi2());
505 svtx_vtx->set_ndof(rave_vtx->getNdf());
506 svtx_vtx->set_position(0, rave_vtx->getPos().X());
507 svtx_vtx->set_position(1, rave_vtx->getPos().Y());
508 svtx_vtx->set_position(2, rave_vtx->getPos().Z());
510 for (
int i = 0; i < 3; i++)
512 for (
int j = 0; j < 3; j++)
514 svtx_vtx->set_error(i, j, rave_vtx->getCov()[i][j]);
518 for (
unsigned int i = 0; i < rave_vtx->getNTracks(); i++)
521 const genfit::Track* rave_track =
522 rave_vtx->getParameters(i)->getTrack();
527 auto iter = gf_track_map.find(rave_track);
528 if (iter != gf_track_map.end())
530 svtx_vtx->insert_track(iter->second);
556 cout <<
PHWHERE <<
" DST Node missing, doing nothing." << endl;
591 cout << m_TrackmapOutNodeName <<
" node added" << endl;
595 m_SvtxVertexMap = findNode::getClass<SvtxVertexMap_v1>(topNode,
"SvtxVertexMap");
600 tb_node->
addNode(vertexes_node);
603 cout <<
"Svtx/SvtxVertexMap node added" << endl;
617 m_TruthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode,
"G4TruthInfo");
620 cout <<
PHWHERE <<
" PHG4TruthInfoContainer node not found on node tree"
631 <<
" node not found on node tree" << endl;
637 cout <<
"PHG4TrackFastSim::GetNodes - node added: " <<
m_PHG4HitsNames[i] << endl;
651 "PHCompositeNode",
"RUN"));
654 std::cout <<
Name() <<
"::"
655 <<
"::" << __PRETTY_FUNCTION__
656 <<
"Run Node missing!" << std::endl;
657 throw std::runtime_error(
658 "Failed to find Run node in PHG4TrackFastSim::CreateNodes");
662 cout << __PRETTY_FUNCTION__ <<
" : ";
688 <<
" node not found on node tree" << endl;
696 std::vector<PHGenFit::Measurement*>& meas_out,
699 TVector3& seed_mom, TMatrixDSym& seed_cov,
const bool do_smearing)
703 seed_cov.ResizeTo(6, 6);
705 seed_pos.SetXYZ(0, 0, 0);
707 seed_cov[0][0] = .1 * .1;
708 seed_cov[1][1] = .1 * .1;
709 seed_cov[2][2] = 30 * 30;
715 seed_mom.SetXYZ(0, 0, 10);
716 for (
int i = 3; i < 6; i++)
723 TVector3 True_mom(particle->
get_px(), particle->
get_py(),
730 const double momSmear = 3. / 180. *
M_PI;
731 const double momMagSmear = 0.1;
735 momMagSmear * True_mom.Mag()));
736 seed_mom.SetTheta(True_mom.Theta() + gsl_ran_gaussian(
m_RandomGenerator, momSmear));
743 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - DEBUG: "
744 <<
"searching for hits from " <<
m_PHG4HitContainer.size() <<
" PHG4Hit nodes" << endl;
748 multimap<double, PHGenFit::Measurement*> ordered_measurements;
754 LogError(
"No m_PHG4HitContainer[i] found!");
766 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - DEBUG: "
770 <<
", detradres = " << detradres
771 <<
", detphires = " << detphires
772 <<
", detlonres = " << detlonres
773 <<
", dethiteff = " << dethiteff
774 <<
", detnoise = " << detnoise
801 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition -adding vertical plane hit ilayer: "
802 << ilayer <<
"; detphires: " << detphires <<
"; detradres: " << detradres <<
" \n";
806 detphires, detradres);
812 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition -adding cylinder hit ilayer: "
813 << ilayer <<
"; detphires: " << detphires <<
"; detlonres : " << detlonres <<
" \n";
817 detphires, detlonres);
825 ordered_measurements.insert(make_pair(hit->
get_avg_t(), meas));
834 for (
auto& pair : ordered_measurements)
836 meas_out.push_back(pair.second);
840 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - measruement at t = " << pair.first <<
" ns: ";
841 pair.second->getMeasurement()->Print();
847 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - meas_out.size = " << meas_out.size() <<
" for "
852 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - seed_pos = "
853 << seed_pos.x() <<
", " << seed_pos.y() <<
". " << seed_pos.z() << endl;
854 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - seed_pos = "
855 << seed_mom.x() <<
", " << seed_mom.y() <<
". " << seed_mom.z() << endl;
856 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - seed_cov = "
857 << sqrt(seed_cov[0][0]) <<
", " << sqrt(seed_cov[1][1]) <<
". " << sqrt(seed_cov[2][2])
859 << sqrt(seed_cov[3][3]) <<
", " << sqrt(seed_cov[4][4]) <<
". " << sqrt(seed_cov[5][5]) << endl;
867 const unsigned int truth_track_id,
868 const unsigned int nmeas,
873 double chi2 = phgf_track->
get_chi2();
874 double ndf = phgf_track->
get_ndf();
876 double pathlenth_from_first_meas = -999999;
877 unique_ptr<genfit::MeasuredStateOnPlane> gf_state(
new genfit::MeasuredStateOnPlane());
894 double pathlenth_orig_from_first_meas = phgf_track->
extrapolateToLine(*gf_state, vtx,
895 TVector3(0., 0., 1.));
899 cout << __PRETTY_FUNCTION__ << __LINE__ <<
" pathlenth_orig_from_first_meas = " << pathlenth_orig_from_first_meas << endl;
901 if (pathlenth_orig_from_first_meas < -999990)
907 TVector3
mom = gf_state->getMom();
908 TVector3
pos = gf_state->getPos();
909 TMatrixDSym cov = gf_state->get6DCov();
918 double dca2d = gf_state->getState()[3];
921 double dca3d = sqrt(dca2d * dca2d + gf_state->getState()[4] * gf_state->getState()[4]);
930 out_track->
set_px(mom.Px());
931 out_track->
set_py(mom.Py());
932 out_track->
set_pz(mom.Pz());
934 out_track->
set_x(pos.X());
935 out_track->
set_y(pos.Y());
936 out_track->
set_z(pos.Z());
937 for (
int i = 0; i < 6; i++)
939 for (
int j = i; j < 6; j++)
950 switch (iter->second.first)
952 case DETECTOR_TYPE::Cylinder:
953 pathlenth_from_first_meas = phgf_track->
extrapolateToCylinder(*gf_state, iter->second.second, TVector3(0., 0., 0.), TVector3(0., 0., 1.), nmeas - 1);
955 case DETECTOR_TYPE::Vertical_Plane:
956 pathlenth_from_first_meas = phgf_track->
extrapolateToPlane(*gf_state, TVector3(0., 0., iter->second.second), TVector3(0, 0., 1.), nmeas - 1);
959 cout <<
"how in the world did you get here??????" << endl;
962 if (pathlenth_from_first_meas < -999990)
967 state->
set_x(gf_state->getPos().x());
968 state->
set_y(gf_state->getPos().y());
969 state->
set_z(gf_state->getPos().z());
971 state->
set_px(gf_state->getMom().x());
972 state->
set_py(gf_state->getMom().y());
973 state->
set_pz(gf_state->getMom().z());
976 for (
int i = 0; i < 6; i++)
978 for (
int j = i; j < 6; j++)
980 state->
set_error(i, j, gf_state->get6DCov()[i][j]);
992 const PHG4Hit* g4hit,
const double phi_resolution,
993 const double r_resolution)
997 TVector3
v(
pos.X(),
pos.Y(), 0);
1000 TVector3
u = v.Cross(TVector3(0, 0, 1));
1001 u = 1 / u.Mag() *
u;
1003 double u_smear = 0.;
1004 double v_smear = 0.;
1010 pos.SetX(g4hit->
get_avg_x() + u_smear * u.X() + v_smear * v.X());
1011 pos.SetY(g4hit->
get_avg_y() + u_smear * u.Y() + v_smear * v.Y());
1028 const PHG4Hit* g4hit,
const double phi_resolution,
1029 const double z_resolution)
1033 TVector3
v(0, 0, 1);
1035 TVector3
u = v.Cross(TVector3(
pos.X(),
pos.Y(), 0));
1036 u = 1 / u.Mag() *
u;
1038 double u_smear = 0.;
1039 double v_smear = 0.;
1067 cov(0, 0) = dxy * dxy;
1068 cov(1, 1) = dxy * dxy;
1069 cov(2, 2) = dz *
dz;
1094 m_ProjectionsMap.insert(make_pair(stateName, make_pair(DETECTOR_TYPE::Vertical_Plane, NAN)));
1098 m_ProjectionsMap.insert(make_pair(stateName, make_pair(DETECTOR_TYPE::Cylinder, NAN)));
1102 cout <<
PHWHERE <<
" Invalid stateName " << stateName << endl;
1104 <<
"These are implemented for cylinders" << endl;
1107 cout << iter << endl;
1110 <<
"These are implemented are for zplanes" << endl;
1113 cout << iter << endl;
1125 cout <<
PHWHERE <<
": " << stateName <<
" is a reserved name, used a different name for your cylinder projection" << endl;
1130 cout <<
PHWHERE <<
": " << stateName <<
" is already a projection, please rename" << endl;
1133 m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Cylinder, radius)));
1142 cout <<
PHWHERE <<
": " << stateName <<
" is a reserved name, used different name for your zplane projection" << endl;
1147 cout <<
PHWHERE <<
": " << stateName <<
" is already a projection, please rename" << endl;
1150 m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Vertical_Plane, zplane)));