533 std::unique_ptr<std::vector<recob::Track>> tcol(
new std::vector<recob::Track>);
534 std::unique_ptr<art::Assns<recob::Track, recob::SpacePoint>> tspassn(
536 std::unique_ptr<art::Assns<recob::Track, recob::Hit>> thassn(
539 unsigned int tcnt = 0;
558 for (
unsigned int ii = 0; ii < mctruthListHandle->size(); ++ii) {
564 std::vector<art::PtrVector<recob::SpacePoint>> spptIn(spptListHandle->begin(),
565 spptListHandle->end());
588 for (
unsigned int ii = 0; ii < mclist.
size(); ++ii) {
590 for (
int jj = 0; jj < mc->NParticles(); ++jj) {
595 <<
"FROM MC TRUTH, the particle's pdg code is: " <<
part.PdgCode()
596 <<
" with energy = " <<
part.E() <<
", with energy = " <<
part.E()
599 <<
"\n (both in Global (not volTPC) coords)";
612 <<
" repMC, covMC are ... \n" 617 TParticlePDG*
part = TDatabasePDG::Instance()->GetParticle(
fPdg);
618 Double_t mass = part->Mass();
620 auto const& tpc = geom->
TPC({0, 0});
622 while (sppt != spptIn.end()) {
631 unsigned int nTailPoints = 0;
632 if (spacepoints.
size() < 5) {
639 <<
"\n\t found " << spacepoints.
size()
640 <<
" 3D spacepoint(s) for this element of std::vector<art:PtrVector> spacepoints. \n";
646 TPrincipal* principal =
new TPrincipal(3,
"ND");
654 std::sort(spacepointss.begin(), spacepointss.end(),
sp_sort_3dz);
655 if (!
fSortDim.compare(
"y")) std::sort(spacepointss.begin(), spacepointss.end(),
sp_sort_3dy);
656 if (!
fSortDim.compare(
"x")) std::sort(spacepointss.begin(), spacepointss.end(),
sp_sort_3dx);
658 for (
unsigned int point = 0; point < spacepointss.size(); ++point) {
659 if (point < (spacepointss.size() - nTailPoints)) {
660 principal->AddRow(spacepointss[point]->XYZ());
663 principal->MakePrincipals();
665 const TVectorD* evals = principal->GetEigenValues();
666 const TMatrixD* evecs = principal->GetEigenVectors();
667 const TVectorD* means = principal->GetMeanValues();
668 const TVectorD* sigmas = principal->GetSigmas();
670 Double_t
tmp[3], tmp2[3];
671 principal->X2P((Double_t*)(means->GetMatrixArray()), tmp);
672 principal->X2P((Double_t*)(sigmas->GetMatrixArray()), tmp2);
673 for (
unsigned int ii = 0; ii < 3; ++ii) {
676 fPCevals[ii] = (Float_t)(evals->GetMatrixArray())[ii];
681 evecs->ExtractRow(ii, 0, w);
686 Double_t tmp3[3], tmp4[3], tmp5[3];
687 principal->X2P((Double_t*)
fPC1, tmp3);
688 principal->X2P((Double_t*)
fPC2, tmp4);
689 principal->X2P((Double_t*)
fPC3, tmp5);
693 fMomStart[0] = spacepointss[spacepointss.size() - 1]->XYZ()[0] - spacepointss[0]->XYZ()[0];
694 fMomStart[1] = spacepointss[spacepointss.size() - 1]->XYZ()[1] - spacepointss[0]->XYZ()[1];
695 fMomStart[2] = spacepointss[spacepointss.size() - 1]->XYZ()[2] - spacepointss[0]->XYZ()[2];
699 TVector3 mom(dEdx *
fMomStart[0], dEdx * fMomStart[1], dEdx * fMomStart[2]);
700 double pmag2 = pow(mom.Mag() + mass, 2. - mass * mass);
701 mom.SetMag(std::sqrt(pmag2));
703 mom.SetMag(1.0 * mom.Mag());
707 bool uncontained(
false);
709 double epsMag(0.001);
713 if (spacepointss[spacepointss.size() - 1]->XYZ()[0] > (2. * tpc.HalfWidth() -
close) ||
714 spacepointss[spacepointss.size() - 1]->XYZ()[0] <
close ||
715 spacepointss[0]->XYZ()[0] > (2. * tpc.HalfWidth() -
close) ||
716 spacepointss[0]->XYZ()[0] <
close ||
717 spacepointss[spacepointss.size() - 1]->XYZ()[1] > (1. * tpc.HalfHeight() -
close) ||
718 (spacepointss[spacepointss.size() - 1]->XYZ()[1] < -1. * tpc.HalfHeight() +
close) ||
719 spacepointss[0]->XYZ()[1] > (1. * tpc.HalfHeight() -
close) ||
720 spacepointss[0]->XYZ()[1] < (-1. * tpc.HalfHeight() +
close) ||
721 spacepointss[spacepointss.size() - 1]->XYZ()[2] > (tpc.Length() -
close) ||
722 spacepointss[spacepointss.size() - 1]->XYZ()[2] <
close ||
723 spacepointss[0]->XYZ()[2] > (tpc.Length() -
close) || spacepointss[0]->XYZ()[2] <
close)
732 mom.SetMag(2.0 * mom.Mag());
733 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
"Uncontained track ... ";
739 <<
"Contained track ... Run " << evt.
run() <<
" Event " << evt.
id().
event();
746 fcont = (int)(!uncontained);
749 unsigned short rePass = rePass0;
751 unsigned short tcnt1(0);
752 while (rePass <= maxPass) {
755 TVector3 momErrFit(momM[0] / 3.0, momM[1] / 3.0,
769 fitTrack.setPDG(
fPdg);
775 std::vector<unsigned int> spptSurvivedIndex;
776 std::vector<unsigned int> spptSkippedIndex;
777 unsigned int ppoint(0);
778 for (
unsigned int point = 0; point < spacepointss.size(); ++point) {
784 principal->X2P((Double_t*)(spacepointss[point]->XYZ()), tmp);
785 sep = std::sqrt(tmp[1] * tmp[1] /
fPCevals[1] + tmp[2] * tmp[2] /
fPCevals[2]);
786 if ((
std::abs(sep) >
fPerpLim) && (point < (spacepointss.size() - nTailPoints)) &&
788 spptSkippedIndex.push_back(point);
793 TVector3 one(spacepointss[point]->XYZ());
794 TVector3 two(spacepointss[ppoint]->XYZ());
795 if (rePass == 2 && uncontained) {
798 fErrScaleMHere = 0.1;
804 else if (rePass == 2 && !uncontained) {}
806 ((one - two).Mag() < epsMag ||
807 ((one - two).Mag() > 8.0 && rePass == 1) ||
808 std::abs(spacepointss[point]->XYZ()[2] - spacepointss[ppoint]->XYZ()[2]) < epsZ ||
809 std::abs(spacepointss[point]->XYZ()[0] - spacepointss[ppoint]->XYZ()[0]) > epsX)) {
810 spptSkippedIndex.push_back(point);
815 point % fDecimateHere &&
824 TVector3 spt3 = (TVector3)(spacepointss[point]->XYZ());
825 std::vector<double> err3;
826 err3.push_back(spacepointss[point]->ErrXYZ()[0]);
827 err3.push_back(spacepointss[point]->ErrXYZ()[2]);
828 err3.push_back(spacepointss[point]->ErrXYZ()[4]);
829 err3.push_back(spacepointss[point]->ErrXYZ()[5]);
847 fth[
fptsNo] = (pointer.Unit()).Angle(pointerPrev.Unit());
858 <<
"ihit xyz..." << spt3[0] <<
"," << spt3[1] <<
"," << spt3[2];
863 spptSurvivedIndex.push_back(point);
871 <<
"Bailing cuz only " <<
fptsNo <<
" spacepoints.";
875 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
"Fitting on " <<
fptsNo <<
" spacepoints.";
888 bool skipFill =
false;
889 std::vector<TMatrixT<double>> hitMeasCov;
890 std::vector<TMatrixT<double>> hitUpdate;
891 std::vector<TMatrixT<double>> hitCov;
892 std::vector<TMatrixT<double>> hitCov7x7;
893 std::vector<TMatrixT<double>> hitState;
894 std::vector<double> hitChi2;
895 std::vector<TVector3> hitPlaneXYZ;
896 std::vector<TVector3> hitPlaneUxUyUz;
897 std::vector<TVector3> hitPlaneU;
898 std::vector<TVector3> hitPlaneV;
905 <<
"just caught a cet::exception: " << e.what()
906 <<
"\nExceptions won't be further handled; skip filling big chunks of the TTree.";
914 std::ostringstream dbgmsg;
915 dbgmsg <<
"Original plane:";
916 planeG.Print(dbgmsg);
918 dbgmsg <<
"Current (fit) reference Plane:";
921 dbgmsg <<
"Last reference Plane:";
925 dbgmsg <<
" => original hit plane (not surprisingly) not current reference Plane!";
927 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit") << dbgmsg.str();
930 hitMeasCov = fitTrack.getHitMeasuredCov();
931 hitUpdate = fitTrack.getHitUpdate();
932 hitCov = fitTrack.getHitCov();
933 hitCov7x7 = fitTrack.getHitCov7x7();
934 hitState = fitTrack.getHitState();
935 hitChi2 = fitTrack.getHitChi2();
936 hitPlaneXYZ = fitTrack.getHitPlaneXYZ();
937 hitPlaneUxUyUz = fitTrack.getHitPlaneUxUyUz();
938 hitPlaneU = fitTrack.getHitPlaneU();
939 hitPlaneV = fitTrack.getHitPlaneV();
940 unsigned int totHits = hitState.size();
943 unsigned int jhit = 0;
944 for (
unsigned int ihit = totHits - 2 * totHits / (2 *
fNumIt);
945 ihit < (totHits - totHits / (2 *
fNumIt));
948 feth[jhit] = (Float_t)(hitMeasCov.at(ihit)[0][0]);
949 fedudw[jhit] = (Float_t)(hitMeasCov.at(ihit)[1][1]);
950 fedvdw[jhit] = (Float_t)(hitMeasCov.at(ihit)[2][2]);
951 feu[jhit] = (Float_t)(hitMeasCov.at(ihit)[3][3]);
952 fev[jhit] = (Float_t)(hitMeasCov.at(ihit)[4][4]);
953 fupdate[jhit] = (Float_t)(hitUpdate.at(ihit)[0][0]);
954 fchi2hit[jhit] = (Float_t)(hitChi2.at(ihit));
964 for (
unsigned int ii = 0; ii < 5; ii++) {
965 stREC->ExtractRow(ii, 0, dum);
967 covREC->ExtractRow(ii, 0, dum2);
968 for (
unsigned int jj = 0; jj < 5; jj++) {
969 fCov0[ii * 5 + jj] = dum2[jj];
977 nfail = fitTrack.getFailedHits();
979 ispptvec = 1 + std::distance(spptB, sppt);
984 <<
"Track3DKalmanSPS about to do tree->Fill(). Chi2/ndf is " <<
chi2 /
ndf <<
".";
986 for (
int ii = 0; ii < 3; ++ii) {
989 fpREC[ii] = hitPlaneUxUyUz.at(totHits - 2 * totHits / (2 *
fNumIt))[ii];
990 fpRECL[ii] = hitPlaneUxUyUz.at(totHits - totHits / (2 *
fNumIt) - 1)[ii];
994 nspptvec = (
unsigned int)spptListHandle->size();
997 std::vector<std::vector<double>> dQdx;
999 std::vector<TMatrixT<double>> hitCovLFP;
1000 std::vector<TVector3> hitPlaneXYZLFP;
1001 std::vector<TVector3> hitPlaneUxUyUzLFP;
1002 std::vector<TVector3> hitPlanePxPyPzLFP;
1003 std::vector<TVector3> hitPlaneULFP;
1004 std::vector<TVector3> hitPlaneVLFP;
1005 std::vector<double> pLFP;
1006 std::vector<TMatrixT<double>> c7x7LFP;
1009 for (
unsigned int ii = 0; ii < totHits / (2 *
fNumIt); ii++) {
1010 pLFP.push_back(1. / hitState.at(totHits - 2 * totHits / (2 *
fNumIt) + ii)[0][0]);
1012 c7x7LFP.push_back(hitCov7x7.at(totHits - 2 * totHits / (2 *
fNumIt) + ii));
1013 hitCovLFP.push_back(hitCov.at(totHits - 2 * totHits / (2 *
fNumIt) + ii));
1014 hitPlaneXYZLFP.push_back(hitPlaneXYZ.at(totHits - 2 * totHits / (2 *
fNumIt) + ii));
1015 hitPlaneUxUyUzLFP.push_back(
1016 hitPlaneUxUyUz.at(totHits - 2 * totHits / (2 *
fNumIt) + ii));
1017 hitPlanePxPyPzLFP.push_back(hitPlaneUxUyUzLFP.back() * pLFP.back());
1018 hitPlaneULFP.push_back(hitPlaneU.at(totHits - 2 * totHits / (2 *
fNumIt) + ii));
1019 hitPlaneVLFP.push_back(hitPlaneV.at(totHits - 2 * totHits / (2 *
fNumIt) + ii));
1025 rotationCov(hitCovLFP.back(), hitPlaneULFP.back(), hitPlaneVLFP.back());
1027 dQdxCalc(hitAssns, spacepointss, hitPlaneUxUyUzLFP.back(), hitPlaneXYZLFP.back()));
1028 fdQdx[ii] = dQdx.back().back();
1038 for (
unsigned int i = 0; i < 5; i++) {
1039 for (
unsigned int j = i; j < 5; j++) {
1040 covVtx(i, j) = hitCovLFP.front()(i, j);
1041 covEnd(i, j) = hitCovLFP.back()(i, j);
1055 if (rePass == 1) tcnt1++;
1056 if (rePass != 1 && tcnt1) tcol->pop_back();
1057 tcol->push_back(the3DTrack);
1060 for (
unsigned int ii = 0; ii < spacepointss.size(); ++ii) {
1061 hits.
insert(hits.
end(), hitAssns.at(ii).begin(), hitAssns.at(ii).end());
1072 for (
unsigned int ind = 0; ind < spptSurvivedIndex.size(); ++ind) {
1081 spacepointss.
begin() + spptSurvivedIndex[ind];
1088 for (
unsigned int ind = 0; ind < spptSkippedIndex.size(); ++ind) {
1090 spacepointss.
begin() + spptSkippedIndex[ind];
1094 std::stable_sort(spacepointss.begin(), spacepointss.end());
1095 std::stable_sort(spacepointssExcise.
begin(), spacepointssExcise.
end());
1096 std::set_union(spacepointssExcise.
begin(),
1097 spacepointssExcise.
end(),
1098 spacepointssExcise.
begin(),
1099 spacepointssExcise.
end(),
1100 spacepointssExcise.
begin());
1103 std::set_difference(spacepointss.begin(),
1105 spacepointssExcise.
begin(),
1106 spacepointssExcise.
end(),
1107 spacepointss.begin());
1108 spacepointss.erase(diffSpptIt, spacepointss.end());
1119 if (uncontained) kick = 0.5;
1120 for (
int ii = 0; ii < 3; ++ii) {
1121 mom[ii] = momM[ii] * kick;
1124 else if (uncontained) {
1125 double unstick(1.0);
1127 for (
int ii = 0; ii < 3; ++ii) {
1128 mom[ii] = momM[ii] * unstick;
1132 for (
int ii = 0; ii < 3; ++ii) {
1133 mom[ii] = 1.1 * momM[ii];
1145 evt.
put(std::move(tcol));
1147 evt.
put(std::move(tspassn));
1150 evt.
put(std::move(thassn));
TMatrixT< Double_t > * stREC
unsigned int getNDF() const
std::vector< double > fMomErr
void setMomLow(Double_t f)
void setMaxUpdate(Double_t f)
TMatrixT< Double_t > * stMCT
static bool sp_sort_3dz(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
typename data_t::iterator iterator
const GFDetPlane & getReferencePlane() const
GFDetPlane getLastPlane() const
void rotationCov(TMatrixT< Double_t > &cov, const TVector3 &u, const TVector3 &v)
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
genf::GFAbsTrackRep * rep
void setMomHigh(Double_t f)
TrackTrajectory::Flags_t Flags_t
std::string fGenieGenModuleLabel
constexpr auto abs(T v)
Returns the absolute value of the argument.
#define MF_LOG_ERROR(category)
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 >> SMatrixSym55
std::vector< double > fPosErr
void Print(std::ostream &out=std::cout) const
const TMatrixT< Double_t > & getState() const
std::string fSpptModuleLabel
static bool sp_sort_3dy(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
std::string fClusterModuleLabel
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
std::vector< double > fMomStart
void push_back(Ptr< U > const &p)
virtual TVector3 getMom(const GFDetPlane &pl)=0
static bool sp_sort_3dx(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
A trajectory in space reconstructed from hits.
void setBlowUpFactor(double f)
Set the blowup factor (see blowUpCovs() )
static bool sp_sort_nsppts(const art::PtrVector< recob::SpacePoint > &h1, const art::PtrVector< recob::SpacePoint > &h2)
double energyLossBetheBloch(const double &mass, const double p)
void setNumIterations(Int_t i)
Set number of iterations for Kalman Filter.
TMatrixT< Double_t > * covREC
void processTrack(GFTrack *)
Performs fit on a GFTrack.
static GFFieldManager * getInstance()
TMatrixT< Double_t > * covMCT
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
const TMatrixT< Double_t > & getCov() const
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
std::string ROOTobjectToString(const ROOTOBJ &obj)
Shortcut to write one ROOT object into a string.
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
iterator insert(iterator position, Ptr< U > const &p)
std::vector< double > dQdxCalc(const art::FindManyP< recob::Hit > &h, const art::PtrVector< recob::SpacePoint > &s, const TVector3 &p, const TVector3 &d)
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
genf::GFAbsTrackRep * repMC
void setInitialDirection(int d)
Sets the inital direction of the track fit (1 for inner to outer, or -1 for outer to inner)...
void init(GFAbsBField *b)
set the magntic field here. Magnetic field classes must be derived from GFAbsBField ...
EventNumber_t event() const
void setErrorScaleSTh(Double_t f)
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
void setErrorScaleMTh(Double_t f)
cet::coded_exception< error, detail::translate > exception