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 << "B0TrackFastSim (DEBUG): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
91 #define LogError(exp) \
92 std::cout << "B0TrackFastSim (ERROR): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
93 #define LogWarning(exp) \
94 std::cout << "B0TrackFastSim (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)
163 std::cerr <<
PHWHERE << std::endl;
177 switch (iter->second.first)
179 case DETECTOR_TYPE::Cylinder:
181 std::string nodename =
"TOWERGEOM_" + iter->first;
189 case DETECTOR_TYPE::Vertical_Plane:
191 std::string towergeonodename =
"TOWERGEOM_" + iter->first;
192 RawTowerGeomContainer* towergeo = findNode::getClass<RawTowerGeomContainer>(topNode, towergeonodename);
195 std::cout <<
PHWHERE <<
" ERROR: Can't find node " << towergeonodename << std::endl;
210 std::cout <<
"invalid state reference: " << iter->second.first << std::endl;
226 std::cout <<
PHWHERE <<
" no Vertex Finder" << std::endl;
248 std::cout <<
"B0TrackFastSim::process_event: " <<
m_EventCnt <<
".\n";
259 LogError(
"m_SvtxTrackMapOut not found!");
263 std::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 std::cout << __PRETTY_FUNCTION__ <<
"Failed to init vertex finder" << std::endl;
408 std::cout << __PRETTY_FUNCTION__ <<
"Failed to init vertex map" << std::endl;
416 std::vector<genfit::GFRaveVertex*> rave_vertices;
417 if (rf_tracks.size() >= 2)
421 std::vector<genfit::Track*> rf_gf_tracks;
422 for (std::vector<PHGenFit::Track*>::iterator
it = rf_tracks.begin();
it != rf_tracks.end(); ++
it)
424 genfit::Track*
track = (*it)->getGenFitTrack();
431 track->getFittedState().getPosMomCov(pos, mom, cov);
433 std::cout <<
"Track getCharge = " << track->getFitStatus()->getCharge() <<
" getChi2 = " << track->getFitStatus()->getChi2() <<
" getNdf = " << track->getFitStatus()->getNdf() << std::endl;
439 rf_gf_tracks.push_back(track);
446 std::cout <<
PHWHERE <<
"GFRaveVertexFactory::findVertices failed!";
452 std::cout << __PRETTY_FUNCTION__ << __LINE__ <<
" rf_tracks = " << rf_tracks.size() <<
" rave_vertices = " << rave_vertices.size() << std::endl;
460 std::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,
492 for (genfit::GFRaveVertex* rave_vtx : rave_vertices)
496 std::cerr <<
PHWHERE << std::endl;
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++)
519 const genfit::Track* rave_track =
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 std::cout <<
PHWHERE <<
" DST Node missing, doing nothing." << std::endl;
589 std::cout << m_TrackmapOutNodeName <<
" node added" << std::endl;
593 m_SvtxVertexMap = findNode::getClass<SvtxVertexMap_v1>(topNode,
"SvtxVertexMap");
598 tb_node->
addNode(vertexes_node);
601 std::cout <<
"Svtx/SvtxVertexMap node added" << std::endl;
615 m_TruthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode,
"G4TruthInfo");
618 std::cout <<
PHWHERE <<
" PHG4TruthInfoContainer node not found on node tree"
629 <<
" node not found on node tree" << std::endl;
635 std::cout <<
"B0TrackFastSim::GetNodes - node added: " <<
m_PHG4HitsNames[i] << std::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 B0TrackFastSim::CreateNodes");
660 std::cout << __PRETTY_FUNCTION__ <<
" : ";
686 <<
" node not found on node tree" << std::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 <<
"B0TrackFastSim::PseudoPatternRecognition - DEBUG: "
742 <<
"searching for hits from " <<
m_PHG4HitContainer.size() <<
" PHG4Hit nodes" << std::endl;
746 std::multimap<double, PHGenFit::Measurement*> ordered_measurements;
752 LogError(
"No m_PHG4HitContainer[i] found!");
764 std::cout <<
"B0TrackFastSim::PseudoPatternRecognition - DEBUG: "
768 <<
", detradres = " << detradres
769 <<
", detphires = " << detphires
770 <<
", detlonres = " << detlonres
771 <<
", dethiteff = " << dethiteff
772 <<
", detnoise = " << detnoise
799 std::cout <<
"B0TrackFastSim::PseudoPatternRecognition -adding vertical plane hit ilayer: "
800 << ilayer <<
"; detphires: " << detphires <<
"; detradres: " << detradres <<
" \n";
804 detphires, detradres);
810 std::cout <<
"B0TrackFastSim::PseudoPatternRecognition -adding cylinder hit ilayer: "
811 << ilayer <<
"; detphires: " << detphires <<
"; detlonres : " << detlonres <<
" \n";
815 detphires, detlonres);
823 ordered_measurements.insert(std::make_pair(hit->
get_avg_t(), meas));
832 for (
auto& pair : ordered_measurements)
834 meas_out.push_back(pair.second);
838 std::cout <<
"B0TrackFastSim::PseudoPatternRecognition - measruement at t = " << pair.first <<
" ns: ";
839 pair.second->getMeasurement()->Print();
845 std::cout <<
"B0TrackFastSim::PseudoPatternRecognition - meas_out.size = " << meas_out.size() <<
" for "
850 std::cout <<
"B0TrackFastSim::PseudoPatternRecognition - seed_pos = "
851 << seed_pos.x() <<
", " << seed_pos.y() <<
". " << seed_pos.z() << std::endl;
852 std::cout <<
"B0TrackFastSim::PseudoPatternRecognition - seed_pos = "
853 << seed_mom.x() <<
", " << seed_mom.y() <<
". " << seed_mom.z() << std::endl;
854 std::cout <<
"B0TrackFastSim::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]) << std::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;
875 std::unique_ptr<genfit::MeasuredStateOnPlane> gf_state(
new genfit::MeasuredStateOnPlane());
892 double pathlenth_orig_from_first_meas = phgf_track->
extrapolateToLine(*gf_state, vtx,
893 TVector3(0., 0., 1.));
897 std::cout << __PRETTY_FUNCTION__ << __LINE__ <<
" pathlenth_orig_from_first_meas = " << pathlenth_orig_from_first_meas << std::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.), 0);
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.), 0);
957 std::cout <<
"how in the world did you get here??????" << std::endl;
960 if (pathlenth_from_first_meas < -999990)
965 state->
set_x(gf_state->getPos().x());
966 state->
set_y(gf_state->getPos().y());
967 state->
set_z(gf_state->getPos().z());
969 state->
set_px(gf_state->getMom().x());
970 state->
set_py(gf_state->getMom().y());
971 state->
set_pz(gf_state->getMom().z());
974 for (
int i = 0; i < 6; i++)
976 for (
int j = i; j < 6; j++)
978 state->
set_error(i, j, gf_state->get6DCov()[i][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.;
1065 cov(0, 0) = dxy * dxy;
1066 cov(1, 1) = dxy * dxy;
1067 cov(2, 2) = dz *
dz;
1090 std::cout <<
"B0TrackFastSim: Adding state name " << stateName << std::endl;
1093 m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Vertical_Plane, NAN)));
1097 m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Cylinder, NAN)));
1101 std::cout <<
PHWHERE <<
" Invalid stateName " << stateName << std::endl;
1102 std::cout << std::endl
1103 <<
"These are implemented for cylinders" << std::endl;
1106 std::cout << iter << std::endl;
1108 std::cout << std::endl
1109 <<
"These are implemented are for zplanes" << std::endl;
1112 std::cout << iter << std::endl;
1124 std::cout <<
PHWHERE <<
": " << stateName <<
" is a reserved name, used a different name for your cylinder projection" << std::endl;
1129 std::cout <<
PHWHERE <<
": " << stateName <<
" is already a projection, please rename" << std::endl;
1132 m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Cylinder, radius)));
1138 std::cout <<
"B0TrackFastSim: Adding z plane state " << stateName <<
" " << zplane << std::endl;
1142 std::cout <<
PHWHERE <<
": " << stateName <<
" is a reserved name, used different name for your zplane projection" << std::endl;
1147 std::cout <<
PHWHERE <<
": " << stateName <<
" is already a projection, please rename" << std::endl;
1150 m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Vertical_Plane, zplane)));