20 #include "TDatabasePDG.h" 23 #include "TPrincipal.h" 37 #include "art_root_io/TFileService.h" 66 const double* xyz1 = h1->
XYZ();
67 const double* xyz2 = h2->
XYZ();
68 return xyz1[2] < xyz2[2];
73 const double* xyz1 = h1->
XYZ();
74 const double* xyz2 = h2->
XYZ();
75 return xyz1[1] < xyz2[1];
80 const double* xyz1 = h1->
XYZ();
81 const double* xyz2 = h2->
XYZ();
82 return xyz1[0] < xyz2[0];
88 const unsigned int s1 = h1.
size();
89 const unsigned int s2 = h2.
size();
105 void rotationCov(TMatrixT<Double_t>& cov,
const TVector3& u,
const TVector3& v);
219 fPosErr = pset.get<std::vector<double>>(
"PosErr3");
220 fMomErr = pset.get<std::vector<double>>(
"MomErr3");
221 fMomStart = pset.get<std::vector<double>>(
"MomStart3");
222 fPerpLim = pset.get<
double>(
"PerpLimit", 1.e6);
223 fDoFit = pset.get<
bool>(
"DoFit",
true);
224 fNumIt = pset.get<
int>(
"NumIt", 5);
226 pset.get<
int>(
"MinNumSppts", 5);
227 fErrScaleS = pset.get<
double>(
"ErrScaleSim", 1.0);
228 fErrScaleM = pset.get<
double>(
"ErrScaleMeas", 1.0);
229 fDecimate = pset.get<
int>(
"DecimateC", 40);
230 fMaxUpdate = pset.get<
double>(
"MaxUpdateC", 0.1);
233 pset.get<
double>(
"DistanceU", 10.0);
234 fMaxUpdateU = pset.get<
double>(
"MaxUpdateU", 0.02);
235 fMomLow = pset.get<
double>(
"MomLow", 0.01);
236 fMomHigh = pset.get<
double>(
"MomHigh", 20.);
237 fPdg = pset.get<
int>(
"PdgCode", -13);
238 fChi2Thresh = pset.get<
double>(
"Chi2HitThresh", 12.0E12);
239 fSortDim = pset.get<std::string>(
"SortDirection",
"z");
240 fMaxPass = pset.get<
int>(
"MaxPass", 2);
242 if (pset.get_if_present(
"GenfPRINT", fGenfPRINT)) {
244 <<
"Parameter 'GenfPRINT' has been deprecated.\n" 245 "Please use the standard message facility to enable GenFit debug output.";
256 produces<std::vector<recob::Track>>();
257 produces<art::Assns<recob::Track, recob::Cluster>>();
258 produces<art::Assns<recob::Track, recob::SpacePoint>>();
259 produces<art::Assns<recob::Track, recob::Hit>>();
266 const double charge(1.0);
267 const double mEE(188.);
268 const double matZ(18.);
269 const double matA(40.);
270 const double matDensity(1.4);
271 const double me(0.000511);
273 double beta = p / std::sqrt(mass * mass + p * p);
274 double gammaSquare = 1. / (1.0 - beta * beta);
276 double dedx = 0.307075 * matDensity * matZ / matA / (beta * beta) * charge * charge;
277 double massRatio = me / mass;
279 double argument = gammaSquare * beta * beta * me * 1.E3 * 2. /
280 ((1.E-6 * mEE) * std::sqrt(1 + 2 * std::sqrt(gammaSquare) * massRatio +
281 massRatio * massRatio));
283 if (mass == 0.0)
return (0.0);
284 if (argument <= exp(beta * beta)) { dedx = 0.; }
286 dedx *= (log(argument) - beta * beta);
288 if (dedx < 0.) dedx = 0.;
295 TVector3 xhat(1.0, 0.0, 0.0);
296 TVector3 yhat(0.0, 1.0, 0.0);
297 TVector3 zhat(0.0, 0.0, 1.0);
298 TVector3
w(u.Cross(v));
300 TVector3 vprime(
w.Cross(xhat));
301 Double_t angle(v.Angle(vprime));
303 uprime.Rotate(angle,
w);
304 if (uprime * xhat < 0) {
305 uprime.Rotate(TMath::Pi(),
w);
306 vprime.Rotate(TMath::Pi(),
w);
307 angle += TMath::Pi();
310 double c = TMath::Cos(angle), s = TMath::Sin(angle);
311 TMatrixT<Double_t> rot(5, 5);
337 std::vector<double> v;
339 double mindist(100.0);
340 auto spptminIt(sppt);
341 while (sppt != s.
end()) {
342 if (((**sppt).XYZ() - loc).Mag() < mindist) {
343 double dist = ((**sppt).XYZ() - loc).Mag();
346 if (dist < mindist) {
349 if (mindist < 0.01)
break;
355 unsigned int ind(std::distance(s.
begin(), spptminIt));
357 std::vector<art::Ptr<recob::Hit>> hitlist = h.at(ind);
359 double wirePitch = 0.;
360 double angleToVert = 0;
366 ihit != hitlist.end();
372 wirePitch = geom->
WirePitch(hit1PlaneID);
373 angleToVert = geom->
Plane(hit1PlaneID).
Wire(0).
ThetaZ(
false) - 0.5 * TMath::Pi();
377 TMath::Abs(TMath::Sin(angleToVert) * dir.Y() + TMath::Cos(angleToVert) * dir.Z());
379 v.push_back(charge / wirePitch / cosgamma);
390 stMCT =
new TMatrixT<Double_t>(5, 1);
391 covMCT =
new TMatrixT<Double_t>(5, 5);
392 stREC =
new TMatrixT<Double_t>(5, 1);
393 covREC =
new TMatrixT<Double_t>(5, 5);
397 fpREC =
new Float_t[4];
400 fCov0 =
new Float_t[25];
421 fPC1 =
new Float_t[3];
422 fPC2 =
new Float_t[3];
423 fPC3 =
new Float_t[3];
430 tree = tfs->make<TTree>(
"GENFITttree",
"GENFITttree");
435 tree->Branch(
"covMCT",
"TMatrixD", &covMCT, 64000, 0);
438 tree->Branch(
"covREC",
fCov0,
"covREC[25]/F");
443 tree->Branch(
"chi2", &
chi2,
"chi2/F");
444 tree->Branch(
"nfail", &
nfail,
"nfail/I");
445 tree->Branch(
"ndf", &
ndf,
"ndf/I");
446 tree->Branch(
"evtNo", &
evtt,
"evtNo/I");
450 tree->Branch(
"trkNo", &
nTrks,
"trkNo/I");
453 tree->Branch(
"shx",
fshx,
"shx[ptsNo]/F");
454 tree->Branch(
"shy",
fshy,
"shy[ptsNo]/F");
455 tree->Branch(
"shz",
fshz,
"shz[ptsNo]/F");
456 tree->Branch(
"sep",
fsep,
"sep[ptsNo]/F");
457 tree->Branch(
"dQdx",
fdQdx,
"dQdx[ptsNo]/F");
458 tree->Branch(
"eshx",
feshx,
"eshx[ptsNo]/F");
459 tree->Branch(
"eshy",
feshy,
"eshy[ptsNo]/F");
460 tree->Branch(
"eshz",
feshz,
"eshz[ptsNo]/F");
461 tree->Branch(
"eshyz",
feshyz,
"eshyz[ptsNo]/F");
462 tree->Branch(
"update",
fupdate,
"update[ptsNo]/F");
464 tree->Branch(
"th",
fth,
"th[ptsNo]/F");
465 tree->Branch(
"eth",
feth,
"eth[ptsNo]/F");
466 tree->Branch(
"edudw",
fedudw,
"edudw[ptsNo]/F");
467 tree->Branch(
"edvdw",
fedvdw,
"edvdw[ptsNo]/F");
468 tree->Branch(
"eu",
feu,
"eu[ptsNo]/F");
469 tree->Branch(
"ev",
fev,
"ev[ptsNo]/F");
474 tree->Branch(
"pcEvec1",
fPC1,
"pcEvec1[3]/F");
475 tree->Branch(
"pcEvec2",
fPC2,
"pcEvec2[3]/F");
476 tree->Branch(
"pcEvec3",
fPC3,
"pcEvec3[3]/F");
478 tree->Branch(
"pMCMom", fpMCMom,
"pMCMom[4]/F");
480 tree->Branch(
"pRECKalF",
fpREC,
"pRECKalF[4]/F");
481 tree->Branch(
"pRECKalL",
fpRECL,
"pRECKalL[4]/F");
546 std::unique_ptr<std::vector<recob::Track>> tcol(
new std::vector<recob::Track>);
547 std::unique_ptr<art::Assns<recob::Track, recob::SpacePoint>> tspassn(
549 std::unique_ptr<art::Assns<recob::Track, recob::Hit>> thassn(
552 unsigned int tcnt = 0;
574 for (
unsigned int ii = 0; ii < mctruthListHandle->size(); ++ii) {
580 std::vector<art::PtrVector<recob::SpacePoint>> spptIn(spptListHandle->begin(),
581 spptListHandle->end());
604 for (
unsigned int ii = 0; ii < mclist.
size(); ++ii) {
606 for (
int jj = 0; jj < mc->
NParticles(); ++jj) {
611 <<
"FROM MC TRUTH, the particle's pdg code is: " <<
part.PdgCode()
612 <<
" with energy = " <<
part.E() <<
", with energy = " <<
part.E()
615 <<
"\n (both in Global (not volTPC) coords)";
628 <<
" repMC, covMC are ... \n" 633 TParticlePDG*
part = TDatabasePDG::Instance()->GetParticle(
fPdg);
634 Double_t mass = part->Mass();
637 while (sppt != spptIn.end()) {
646 unsigned int nTailPoints = 0;
647 if (spacepoints.
size() < 5) {
654 <<
"\n\t found " << spacepoints.
size()
655 <<
" 3D spacepoint(s) for this element of std::vector<art:PtrVector> spacepoints. \n";
661 TPrincipal* principal =
new TPrincipal(3,
"ND");
673 for (
unsigned int point = 0; point < spacepointss.
size(); ++point) {
674 if (point < (spacepointss.
size() - nTailPoints)) {
675 principal->AddRow(spacepointss[point]->XYZ());
678 principal->MakePrincipals();
680 const TVectorD* evals = principal->GetEigenValues();
681 const TMatrixD* evecs = principal->GetEigenVectors();
682 const TVectorD* means = principal->GetMeanValues();
683 const TVectorD* sigmas = principal->GetSigmas();
685 Double_t
tmp[3], tmp2[3];
686 principal->X2P((Double_t*)(means->GetMatrixArray()), tmp);
687 principal->X2P((Double_t*)(sigmas->GetMatrixArray()), tmp2);
688 for (
unsigned int ii = 0; ii < 3; ++ii) {
691 fPCevals[ii] = (Float_t)(evals->GetMatrixArray())[ii];
696 evecs->ExtractRow(ii, 0, w);
701 Double_t tmp3[3], tmp4[3], tmp5[3];
702 principal->X2P((Double_t*)
fPC1, tmp3);
703 principal->X2P((Double_t*)
fPC2, tmp4);
704 principal->X2P((Double_t*)
fPC3, tmp5);
708 fMomStart[0] = spacepointss[spacepointss.
size() - 1]->XYZ()[0] - spacepointss[0]->XYZ()[0];
709 fMomStart[1] = spacepointss[spacepointss.
size() - 1]->XYZ()[1] - spacepointss[0]->XYZ()[1];
710 fMomStart[2] = spacepointss[spacepointss.
size() - 1]->XYZ()[2] - spacepointss[0]->XYZ()[2];
714 TVector3 mom(dEdx *
fMomStart[0], dEdx * fMomStart[1], dEdx * fMomStart[2]);
715 double pmag2 = pow(mom.Mag() + mass, 2. - mass * mass);
716 mom.SetMag(std::sqrt(pmag2));
718 mom.SetMag(1.0 * mom.Mag());
722 bool uncontained(
false);
724 double epsMag(0.001);
729 spacepointss[spacepointss.
size() - 1]->XYZ()[0] < close ||
731 spacepointss[0]->XYZ()[0] < close ||
737 spacepointss[spacepointss.
size() - 1]->XYZ()[2] < close ||
739 spacepointss[0]->XYZ()[2] <
close)
748 mom.SetMag(2.0 * mom.Mag());
749 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
"Uncontained track ... ";
755 <<
"Contained track ... Run " << evt.
run() <<
" Event " << evt.
id().
event();
762 fcont = (int)(!uncontained);
765 unsigned short rePass = rePass0;
767 unsigned short tcnt1(0);
768 while (rePass <= maxPass) {
771 TVector3 momErrFit(momM[0] / 3.0, momM[1] / 3.0,
791 std::vector<unsigned int> spptSurvivedIndex;
792 std::vector<unsigned int> spptSkippedIndex;
793 unsigned int ppoint(0);
794 for (
unsigned int point = 0; point < spacepointss.
size(); ++point) {
800 principal->X2P((Double_t*)(spacepointss[point]->XYZ()), tmp);
801 sep = std::sqrt(tmp[1] * tmp[1] /
fPCevals[1] + tmp[2] * tmp[2] /
fPCevals[2]);
804 spptSkippedIndex.push_back(point);
809 TVector3 one(spacepointss[point]->XYZ());
810 TVector3 two(spacepointss[ppoint]->XYZ());
811 if (rePass == 2 && uncontained) {
814 fErrScaleMHere = 0.1;
820 else if (rePass == 2 && !uncontained) {}
822 ((one - two).Mag() < epsMag ||
823 ((one - two).Mag() > 8.0 && rePass == 1) ||
824 std::abs(spacepointss[point]->XYZ()[2] - spacepointss[ppoint]->XYZ()[2]) < epsZ ||
825 std::abs(spacepointss[point]->XYZ()[0] - spacepointss[ppoint]->XYZ()[0]) > epsX)) {
826 spptSkippedIndex.push_back(point);
831 point % fDecimateHere &&
840 TVector3 spt3 = (TVector3)(spacepointss[point]->XYZ());
841 std::vector<double> err3;
842 err3.push_back(spacepointss[point]->ErrXYZ()[0]);
843 err3.push_back(spacepointss[point]->ErrXYZ()[2]);
844 err3.push_back(spacepointss[point]->ErrXYZ()[4]);
845 err3.push_back(spacepointss[point]->ErrXYZ()[5]);
863 fth[
fptsNo] = (pointer.Unit()).Angle(pointerPrev.Unit());
874 <<
"ihit xyz..." << spt3[0] <<
"," << spt3[1] <<
"," << spt3[2];
879 spptSurvivedIndex.push_back(point);
887 <<
"Bailing cuz only " <<
fptsNo <<
" spacepoints.";
891 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
"Fitting on " <<
fptsNo <<
" spacepoints.";
904 bool skipFill =
false;
905 std::vector<TMatrixT<double>> hitMeasCov;
906 std::vector<TMatrixT<double>> hitUpdate;
907 std::vector<TMatrixT<double>> hitCov;
908 std::vector<TMatrixT<double>> hitCov7x7;
909 std::vector<TMatrixT<double>> hitState;
910 std::vector<double> hitChi2;
911 std::vector<TVector3> hitPlaneXYZ;
912 std::vector<TVector3> hitPlaneUxUyUz;
913 std::vector<TVector3> hitPlaneU;
914 std::vector<TVector3> hitPlaneV;
921 <<
"just caught a cet::exception: " << e.what()
922 <<
"\nExceptions won't be further handled; skip filling big chunks of the TTree.";
930 std::ostringstream dbgmsg;
931 dbgmsg <<
"Original plane:";
932 planeG.
Print(dbgmsg);
934 dbgmsg <<
"Current (fit) reference Plane:";
937 dbgmsg <<
"Last reference Plane:";
941 dbgmsg <<
" => original hit plane (not surprisingly) not current reference Plane!";
943 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit") << dbgmsg.str();
956 unsigned int totHits = hitState.size();
959 unsigned int jhit = 0;
960 for (
unsigned int ihit = totHits - 2 * totHits / (2 *
fNumIt);
961 ihit < (totHits - totHits / (2 *
fNumIt));
964 feth[jhit] = (Float_t)(hitMeasCov.at(ihit)[0][0]);
965 fedudw[jhit] = (Float_t)(hitMeasCov.at(ihit)[1][1]);
966 fedvdw[jhit] = (Float_t)(hitMeasCov.at(ihit)[2][2]);
967 feu[jhit] = (Float_t)(hitMeasCov.at(ihit)[3][3]);
968 fev[jhit] = (Float_t)(hitMeasCov.at(ihit)[4][4]);
969 fupdate[jhit] = (Float_t)(hitUpdate.at(ihit)[0][0]);
970 fchi2hit[jhit] = (Float_t)(hitChi2.at(ihit));
980 for (
unsigned int ii = 0; ii < 5; ii++) {
981 stREC->ExtractRow(ii, 0, dum);
983 covREC->ExtractRow(ii, 0, dum2);
984 for (
unsigned int jj = 0; jj < 5; jj++) {
985 fCov0[ii * 5 + jj] = dum2[jj];
995 ispptvec = 1 + std::distance(spptB, sppt);
1000 <<
"Track3DKalmanSPS about to do tree->Fill(). Chi2/ndf is " <<
chi2 /
ndf <<
".";
1001 fpMCMom[3] = MCMomentum.Mag();
1002 for (
int ii = 0; ii < 3; ++ii) {
1005 fpREC[ii] = hitPlaneUxUyUz.at(totHits - 2 * totHits / (2 *
fNumIt))[ii];
1006 fpRECL[ii] = hitPlaneUxUyUz.at(totHits - totHits / (2 *
fNumIt) - 1)[ii];
1010 nspptvec = (
unsigned int)spptListHandle->size();
1013 std::vector<std::vector<double>> dQdx;
1015 std::vector<TMatrixT<double>> hitCovLFP;
1016 std::vector<TVector3> hitPlaneXYZLFP;
1017 std::vector<TVector3> hitPlaneUxUyUzLFP;
1018 std::vector<TVector3> hitPlanePxPyPzLFP;
1019 std::vector<TVector3> hitPlaneULFP;
1020 std::vector<TVector3> hitPlaneVLFP;
1021 std::vector<double> pLFP;
1022 std::vector<TMatrixT<double>> c7x7LFP;
1025 for (
unsigned int ii = 0; ii < totHits / (2 *
fNumIt); ii++) {
1026 pLFP.push_back(1. / hitState.at(totHits - 2 * totHits / (2 *
fNumIt) + ii)[0][0]);
1028 c7x7LFP.push_back(hitCov7x7.at(totHits - 2 * totHits / (2 *
fNumIt) + ii));
1029 hitCovLFP.push_back(hitCov.at(totHits - 2 * totHits / (2 *
fNumIt) + ii));
1030 hitPlaneXYZLFP.push_back(hitPlaneXYZ.at(totHits - 2 * totHits / (2 *
fNumIt) + ii));
1031 hitPlaneUxUyUzLFP.push_back(
1032 hitPlaneUxUyUz.at(totHits - 2 * totHits / (2 *
fNumIt) + ii));
1033 hitPlanePxPyPzLFP.push_back(hitPlaneUxUyUzLFP.back() * pLFP.back());
1034 hitPlaneULFP.push_back(hitPlaneU.at(totHits - 2 * totHits / (2 *
fNumIt) + ii));
1035 hitPlaneVLFP.push_back(hitPlaneV.at(totHits - 2 * totHits / (2 *
fNumIt) + ii));
1041 rotationCov(hitCovLFP.back(), hitPlaneULFP.back(), hitPlaneVLFP.back());
1043 dQdxCalc(hitAssns, spacepointss, hitPlaneUxUyUzLFP.back(), hitPlaneXYZLFP.back()));
1044 fdQdx[ii] = dQdx.back().back();
1054 for (
unsigned int i = 0; i < 5; i++) {
1055 for (
unsigned int j = i; j < 5; j++) {
1056 covVtx(i, j) = hitCovLFP.front()(i, j);
1057 covEnd(i, j) = hitCovLFP.back()(i, j);
1071 if (rePass == 1) tcnt1++;
1072 if (rePass != 1 && tcnt1) tcol->pop_back();
1073 tcol->push_back(the3DTrack);
1076 for (
unsigned int ii = 0; ii < spacepointss.
size(); ++ii) {
1077 hits.
insert(hits.
end(), hitAssns.at(ii).begin(), hitAssns.at(ii).end());
1088 for (
unsigned int ind = 0; ind < spptSurvivedIndex.size(); ++ind) {
1097 spacepointss.
begin() + spptSurvivedIndex[ind];
1104 for (
unsigned int ind = 0; ind < spptSkippedIndex.size(); ++ind) {
1106 spacepointss.
begin() + spptSkippedIndex[ind];
1110 std::stable_sort(spacepointss.
begin(), spacepointss.
end());
1111 std::stable_sort(spacepointssExcise.
begin(), spacepointssExcise.
end());
1112 std::set_union(spacepointssExcise.
begin(),
1113 spacepointssExcise.
end(),
1114 spacepointssExcise.
begin(),
1115 spacepointssExcise.
end(),
1116 spacepointssExcise.
begin());
1119 std::set_difference(spacepointss.
begin(),
1121 spacepointssExcise.
begin(),
1122 spacepointssExcise.
end(),
1123 spacepointss.
begin());
1124 spacepointss.
erase(diffSpptIt, spacepointss.
end());
1135 if (uncontained) kick = 0.5;
1136 for (
int ii = 0; ii < 3; ++ii) {
1137 mom[ii] = momM[ii] * kick;
1140 else if (uncontained) {
1141 double unstick(1.0);
1143 for (
int ii = 0; ii < 3; ++ii) {
1144 mom[ii] = momM[ii] * unstick;
1148 for (
int ii = 0; ii < 3; ++ii) {
1149 mom[ii] = 1.1 * momM[ii];
1161 evt.
put(std::move(tcol));
1163 evt.
put(std::move(tspassn));
1166 evt.
put(std::move(thassn));
TMatrixT< Double_t > * stREC
unsigned int getNDF() const
std::vector< double > fMomErr
void setMomLow(Double_t f)
std::string GetLArTPCVolumeName(TPCID const &tpcid=tpc_zero) const
Return the name of specified LAr TPC volume.
std::vector< TVector3 > getHitPlaneUxUyUz()
void setMaxUpdate(Double_t f)
geo::SigType_t SignalType() const
Signal type for the plane of the hit.
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
std::string fG4ModuleLabel
GFDetPlane getLastPlane() const
void rotationCov(TMatrixT< Double_t > &cov, const TVector3 &u, const TVector3 &v)
WireGeo const & Wire(unsigned int iwire) const
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Length_t DetHalfWidth(TPCID const &tpcid=tpc_zero) const
Returns the half width of the active volume of the specified TPC.
The data type to uniquely identify a Plane.
genf::GFAbsTrackRep * rep
void setMomHigh(Double_t f)
iterator erase(iterator position)
TrackTrajectory::Flags_t Flags_t
std::string fGenieGenModuleLabel
std::vector< Double_t > getHitChi2()
constexpr auto abs(T v)
Returns the absolute value of the argument.
#define MF_LOG_ERROR(category)
std::vector< TMatrixT< Double_t > > getHitCov7x7()
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 >> SMatrixSym55
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
geo::WireID const & WireID() const
Initial tdc tick for hit.
std::vector< double > fPosErr
void Print(std::ostream &out=std::cout) const
const TMatrixT< Double_t > & getState() const
std::vector< TMatrixT< Double_t > > getHitState()
std::vector< TVector3 > getHitPlaneU()
std::string fSpptModuleLabel
static bool sp_sort_3dy(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
std::string fClusterModuleLabel
void produce(art::Event &evt)
Length_t DetLength(TPCID const &tpcid=tpc_zero) const
Returns the length of the active volume of the specified TPC.
double ThetaZ() const
Returns angle of wire with respect to z axis in the Y-Z plane in radians.
Base Class for genfit track representations. Defines interface for track parameterizations.
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
typename data_t::const_iterator const_iterator
std::vector< double > fMomStart
#define DEFINE_ART_MODULE(klass)
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.
enum geo::_plane_sigtype SigType_t
Enumerate the possible plane projections.
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.
int getFailedHits(int repId=-1)
return the number of failed Hits in track fit repId == -1 will use cardinal rep
TMatrixT< Double_t > * covREC
void processTrack(GFTrack *)
Performs fit on a GFTrack.
Provides recob::Track data product.
static GFFieldManager * getInstance()
TMatrixT< Double_t > * covMCT
const Double32_t * XYZ() const
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
Declaration of cluster object.
const TMatrixT< Double_t > & getCov() const
std::vector< TVector3 > getHitPlaneV()
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.
Encapsulate the geometry of a wire.
const simb::MCParticle & GetParticle(int i) const
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
Utility object to perform functions of association.
Encapsulate the construction of a single detector plane.
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
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
std::vector< TMatrixT< Double_t > > getHitUpdate()
Length_t DetHalfHeight(TPCID const &tpcid=tpc_zero) const
Returns the half height of the active volume of the specified TPC.
void setErrorScaleSTh(Double_t f)
2D representation of charge deposited in the TDC/wire plane
void addHit(GFAbsRecoHit *theHit)
deprecated!
std::vector< TMatrixT< Double_t > > getHitCov()
Track3DKalmanSPS(fhicl::ParameterSet const &pset)
#define MF_LOG_WARNING(category)
std::vector< TMatrixT< Double_t > > getHitMeasuredCov()
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
art framework interface to geometry description
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
Signal from collection planes.
std::vector< TVector3 > getHitPlaneXYZ()