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>
26 #include <calobase/RawTowerGeom.h>
27 #include <calobase/RawTowerGeomContainer.h>
29 #include <phparameter/PHParameters.h>
31 #include <phfield/PHFieldUtility.h>
51 #include <GenFit/AbsMeasurement.h>
52 #include <GenFit/EventDisplay.h>
53 #include <GenFit/MeasuredStateOnPlane.h>
54 #include <GenFit/RKTrackRep.h>
56 #include <GenFit/FitStatus.h>
57 #include <GenFit/GFRaveTrackParameters.h>
58 #include <GenFit/GFRaveVertex.h>
59 #include <GenFit/GFRaveVertexFactory.h>
62 #include <TMatrixDSymfwd.h>
63 #include <TMatrixTSym.h>
64 #include <TMatrixTUtils.h>
67 #include <TVectorDfwd.h>
70 #include <gsl/gsl_randist.h>
71 #include <gsl/gsl_rng.h>
87 #define LogDebug(exp) \
88 if (Verbosity()) std::cout << "PHG4TrackFastSim (DEBUG): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
89 #define LogError(exp) \
90 std::cout << "PHG4TrackFastSim (ERROR): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
91 #define LogWarning(exp) \
92 std::cout << "PHG4TrackFastSim (WARNING): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
105 , m_RaveVertexFactory(nullptr)
106 , m_TruthContainer(nullptr)
107 , m_SvtxTrackMapOut(nullptr)
108 , m_SvtxVertexMap(nullptr)
109 , m_SubTopnodeName(
"SVTX")
110 , m_TrackmapOutNodeName(
"SvtxTrackMap")
111 , m_VertexingMethod(
"kalman-smoothing:1")
112 , m_FitAlgoName(
"DafRef")
113 , m_VertexMinNdf(10.)
114 , m_VertexXYResolution(50
E-4)
115 , m_VertexZResolution(50
E-4)
117 , m_PrimaryAssumptionPid(211)
118 , m_SmearingFlag(
true)
119 , m_DoEvtDisplayFlag(
false)
120 , m_UseVertexInFittingFlag(
true)
121 , m_PrimaryTrackingFlag(1)
122 , m_DoVertexingFlag(
false)
125 m_RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
173 if (isfinite(iter->second.second))
177 switch (iter->second.first)
179 case DETECTOR_TYPE::Cylinder:
181 string nodename =
"TOWERGEOM_" + iter->first;
189 case DETECTOR_TYPE::Vertical_Plane:
191 string towergeonodename =
"TOWERGEOM_" + iter->first;
192 RawTowerGeomContainer* towergeo = findNode::getClass<RawTowerGeomContainer>(topNode, towergeonodename);
195 cout <<
PHWHERE <<
" ERROR: Can't find node " << towergeonodename << endl;
210 cout <<
"invalid state reference: " << iter->second.first << endl;
221 cout <<
PHWHERE <<
" no Vertex Finder" << endl;
248 std::cout <<
"PHG4TrackFastSim::process_event: " <<
m_EventCnt <<
".\n";
259 LogError(
"m_SvtxTrackMapOut not found!");
263 vector<PHGenFit::Track*> rf_tracks;
266 TVector3 vtxPoint(truthVtx->
get_x(), truthVtx->
get_y(), truthVtx->
get_z());
288 itr != itr_range.second; ++itr)
292 TVector3 seed_pos(vtxPoint.x(), vtxPoint.y(), vtxPoint.z());
293 TVector3 seed_mom(0, 0, 0);
294 TMatrixDSym seed_cov(6);
297 std::vector<PHGenFit::Measurement*> measurements;
308 measurements.push_back(vtx_meas);
314 svtx_track_out.get(),
318 if (measurements.size() < 3)
323 std::cout <<
"event: " << m_EventCnt <<
" : measurements.size() < 3"
328 for (
unsigned int im = 0; im < measurements.size(); im++)
330 delete measurements[im]->getMeasurement();
331 delete measurements[im];
357 rf_tracks.push_back(track);
367 if (fitting_err != 0)
372 std::cout <<
"event: " << m_EventCnt
373 <<
" : fitting_err != 0, next track."
379 TVector3 vtx(vtxPoint.x(), vtxPoint.y(), vtxPoint.z());
382 measurements.size(), vtx);
392 const unsigned int track_id = m_SvtxTrackMapOut->insert(svtx_track_out.get())->get_id();
403 cout << __PRETTY_FUNCTION__ <<
"Failed to init vertex finder" << endl;
408 cout << __PRETTY_FUNCTION__ <<
"Failed to init vertex map" << endl;
416 vector<genfit::GFRaveVertex*> rave_vertices;
417 if (rf_tracks.size() >= 2)
421 vector<genfit::Track*> rf_gf_tracks;
422 for (std::vector<PHGenFit::Track*>::iterator
it = rf_tracks.begin();
it != rf_tracks.end(); ++
it)
439 rf_gf_tracks.push_back(track);
446 std::cout <<
PHWHERE <<
"GFRaveVertexFactory::findVertices failed!";
452 cout << __PRETTY_FUNCTION__ << __LINE__ <<
" rf_tracks = " << rf_tracks.size() <<
" rave_vertices = " << rave_vertices.size() << endl;
460 vector<genfit::Track*> rf_gf_tracks;
461 for (std::vector<PHGenFit::Track*>::iterator
it = rf_tracks.begin();
it != rf_tracks.end(); ++
it)
463 rf_gf_tracks.push_back((*it)->getGenFitTrack());
469 for (std::vector<PHGenFit::Track*>::iterator
it = rf_tracks.begin();
it != rf_tracks.end(); ++
it)
489 const std::vector<genfit::GFRaveVertex*>& rave_vertices,
502 svtx_vtx->set_chisq(rave_vtx->getChi2());
503 svtx_vtx->set_ndof(rave_vtx->getNdf());
504 svtx_vtx->set_position(0, rave_vtx->getPos().X());
505 svtx_vtx->set_position(1, rave_vtx->getPos().Y());
506 svtx_vtx->set_position(2, rave_vtx->getPos().Z());
508 for (
int i = 0;
i < 3;
i++)
510 for (
int j = 0;
j < 3;
j++)
512 svtx_vtx->set_error(
i,
j, rave_vtx->getCov()[
i][
j]);
516 for (
unsigned int i = 0;
i < rave_vtx->getNTracks();
i++)
520 rave_vtx->getParameters(
i)->getTrack();
525 auto iter = gf_track_map.find(rave_track);
526 if (iter != gf_track_map.end())
528 svtx_vtx->insert_track(iter->second);
554 cout <<
PHWHERE <<
" DST Node missing, doing nothing." << endl;
581 if (!m_SvtxTrackMapOut)
589 cout << m_TrackmapOutNodeName <<
" node added" << endl;
593 m_SvtxVertexMap = findNode::getClass<SvtxVertexMap_v1>(topNode,
"SvtxVertexMap");
598 tb_node->
addNode(vertexes_node);
601 cout <<
"Svtx/SvtxVertexMap node added" << endl;
615 m_TruthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode,
"G4TruthInfo");
618 cout <<
PHWHERE <<
" PHG4TruthInfoContainer node not found on node tree"
629 <<
" node not found on node tree" << endl;
635 cout <<
"PHG4TrackFastSim::GetNodes - node added: " <<
m_PHG4HitsNames[
i] << endl;
649 "PHCompositeNode",
"RUN"));
652 std::cout <<
Name() <<
"::"
653 <<
"::" << __PRETTY_FUNCTION__
654 <<
"Run Node missing!" << std::endl;
655 throw std::runtime_error(
656 "Failed to find Run node in PHG4TrackFastSim::CreateNodes");
660 cout << __PRETTY_FUNCTION__ <<
" : ";
683 if (!m_SvtxTrackMapOut && m_EventCnt < 2)
685 cout <<
PHWHERE << m_TrackmapOutNodeName
686 <<
" node not found on node tree" << endl;
694 std::vector<PHGenFit::Measurement*>& meas_out,
697 TVector3& seed_mom, TMatrixDSym& seed_cov,
const bool do_smearing)
701 seed_cov.ResizeTo(6, 6);
703 seed_pos.SetXYZ(0, 0, 0);
705 seed_cov[0][0] = .1 * .1;
706 seed_cov[1][1] = .1 * .1;
707 seed_cov[2][2] = 30 * 30;
713 seed_mom.SetXYZ(0, 0, 10);
714 for (
int i = 3;
i < 6;
i++)
721 TVector3 True_mom(particle->
get_px(), particle->
get_py(),
728 const double momSmear = 3. / 180. * M_PI;
729 const double momMagSmear = 0.1;
733 momMagSmear * True_mom.Mag()));
734 seed_mom.SetTheta(True_mom.Theta() + gsl_ran_gaussian(
m_RandomGenerator, momSmear));
741 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - DEBUG: "
742 <<
"searching for hits from " <<
m_PHG4HitContainer.size() <<
" PHG4Hit nodes" << endl;
746 multimap<double, PHGenFit::Measurement*> ordered_measurements;
752 LogError(
"No m_PHG4HitContainer[i] found!");
764 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - DEBUG: "
768 <<
", detradres = " << detradres
769 <<
", detphires = " << detphires
770 <<
", detlonres = " << detlonres
771 <<
", dethiteff = " << dethiteff
772 <<
", detnoise = " << detnoise
799 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition -adding vertical plane hit ilayer: "
800 << ilayer <<
"; detphires: " << detphires <<
"; detradres: " << detradres <<
" \n";
804 detphires, detradres);
810 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition -adding cylinder hit ilayer: "
811 << ilayer <<
"; detphires: " << detphires <<
"; detlonres : " << detlonres <<
" \n";
815 detphires, detlonres);
823 ordered_measurements.insert(make_pair(hit->
get_avg_t(), meas));
832 for (
auto& pair : ordered_measurements)
834 meas_out.push_back(pair.second);
838 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - measruement at t = " << pair.first <<
" ns: ";
839 pair.second->getMeasurement()->Print();
845 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - meas_out.size = " << meas_out.size() <<
" for "
850 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - seed_pos = "
851 << seed_pos.x() <<
", " << seed_pos.y() <<
". " << seed_pos.z() << endl;
852 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - seed_pos = "
853 << seed_mom.x() <<
", " << seed_mom.y() <<
". " << seed_mom.z() << endl;
854 std::cout <<
"PHG4TrackFastSim::PseudoPatternRecognition - seed_cov = "
855 << sqrt(seed_cov[0][0]) <<
", " << sqrt(seed_cov[1][1]) <<
". " << sqrt(seed_cov[2][2])
857 << sqrt(seed_cov[3][3]) <<
", " << sqrt(seed_cov[4][4]) <<
". " << sqrt(seed_cov[5][5]) << endl;
865 const unsigned int truth_track_id,
866 const unsigned int nmeas,
871 double chi2 = phgf_track->
get_chi2();
872 double ndf = phgf_track->
get_ndf();
874 double pathlenth_from_first_meas = -999999;
892 double pathlenth_orig_from_first_meas = phgf_track->
extrapolateToLine(*gf_state, vtx,
893 TVector3(0., 0., 1.));
897 cout << __PRETTY_FUNCTION__ << __LINE__ <<
" pathlenth_orig_from_first_meas = " << pathlenth_orig_from_first_meas << endl;
899 if (pathlenth_orig_from_first_meas < -999990)
905 TVector3 mom = gf_state->
getMom();
906 TVector3 pos = gf_state->
getPos();
907 TMatrixDSym cov = gf_state->
get6DCov();
916 double dca2d = gf_state->
getState()[3];
919 double dca3d = sqrt(dca2d * dca2d + gf_state->
getState()[4] * gf_state->
getState()[4]);
928 out_track->
set_px(mom.Px());
929 out_track->
set_py(mom.Py());
930 out_track->
set_pz(mom.Pz());
932 out_track->
set_x(pos.X());
933 out_track->
set_y(pos.Y());
934 out_track->
set_z(pos.Z());
935 for (
int i = 0;
i < 6;
i++)
937 for (
int j =
i;
j < 6;
j++)
948 switch (iter->second.first)
950 case DETECTOR_TYPE::Cylinder:
951 pathlenth_from_first_meas = phgf_track->
extrapolateToCylinder(*gf_state, iter->second.second, TVector3(0., 0., 0.), TVector3(0., 0., 1.), nmeas - 1);
953 case DETECTOR_TYPE::Vertical_Plane:
954 pathlenth_from_first_meas = phgf_track->
extrapolateToPlane(*gf_state, TVector3(0., 0., iter->second.second), TVector3(0, 0., 1.), nmeas - 1);
957 cout <<
"how in the world did you get here??????" << endl;
960 if (pathlenth_from_first_meas < -999990)
974 for (
int i = 0;
i < 6;
i++)
976 for (
int j =
i;
j < 6;
j++)
990 const PHG4Hit* g4hit,
const double phi_resolution,
991 const double r_resolution)
995 TVector3
v(pos.X(), pos.Y(), 0);
998 TVector3
u = v.Cross(TVector3(0, 0, 1));
1001 double u_smear = 0.;
1002 double v_smear = 0.;
1008 pos.SetX(g4hit->
get_avg_x() + u_smear * u.X() + v_smear * v.X());
1009 pos.SetY(g4hit->
get_avg_y() + u_smear * u.Y() + v_smear * v.Y());
1026 const PHG4Hit* g4hit,
const double phi_resolution,
1027 const double z_resolution)
1031 TVector3
v(0, 0, 1);
1033 TVector3
u = v.Cross(TVector3(pos.X(), pos.Y(), 0));
1034 u = 1 / u.Mag() *
u;
1036 double u_smear = 0.;
1037 double v_smear = 0.;
1043 pos.SetX(g4hit->
get_avg_x() + u_smear * u.X());
1044 pos.SetY(g4hit->
get_avg_y() + u_smear * u.Y());
1065 cov(0, 0) = dxy * dxy;
1066 cov(1, 1) = dxy * dxy;
1067 cov(2, 2) = dz *
dz;
1092 m_ProjectionsMap.insert(make_pair(stateName, make_pair(DETECTOR_TYPE::Vertical_Plane, NAN)));
1096 m_ProjectionsMap.insert(make_pair(stateName, make_pair(DETECTOR_TYPE::Cylinder, NAN)));
1100 cout <<
PHWHERE <<
" Invalid stateName " << stateName << endl;
1102 <<
"These are implemented for cylinders" << endl;
1105 cout << iter << endl;
1108 <<
"These are implemented are for zplanes" << endl;
1111 cout << iter << endl;
1123 cout <<
PHWHERE <<
": " << stateName <<
" is a reserved name, used a different name for your cylinder projection" << endl;
1128 cout <<
PHWHERE <<
": " << stateName <<
" is already a projection, please rename" << endl;
1131 m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Cylinder, radius)));
1140 cout <<
PHWHERE <<
": " << stateName <<
" is a reserved name, used different name for your zplane projection" << endl;
1145 cout <<
PHWHERE <<
": " << stateName <<
" is already a projection, please rename" << endl;
1148 m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Vertical_Plane, zplane)));