207 CLHEP::RandGaussQ gauss(
fEngine);
213 std::unique_ptr<std::vector<recob::Track>> tcol(
new std::vector<recob::Track>);
214 std::unique_ptr<std::vector<recob::SpacePoint>> spcol(
new std::vector<recob::SpacePoint>);
215 std::unique_ptr<art::Assns<recob::Track, recob::SpacePoint>> tspassn(
217 std::unique_ptr<art::Assns<recob::Track, recob::Cluster>> tcassn(
219 std::unique_ptr<art::Assns<recob::Track, recob::Hit>> thassn(
223 TString tpcName = geom->
TPC({0, 0}).ActiveVolume()->GetName();
235 for (
unsigned int ii = 0; ii < mctruthListHandle->size(); ++ii) {
242 std::vector<art::Ptr<recob::SpacePoint>> spacepoints;
246 <<
"There are " << trackListHandle->size()
247 <<
" Track3Dreco/SpacePt tracks/groups (whichever) in this event.";
253 for (
unsigned int ii = 0; ii < trackListHandle->size(); ++ii) {
271 for (
unsigned int ii = 0; ii < mclist.
size(); ++ii) {
274 for (
int jj = 0; jj < mc->NParticles(); ++jj) {
277 <<
"FROM MC TRUTH, the particle's pdg code is: " <<
part.PdgCode()
278 <<
" with energy = " <<
part.E() <<
", with energy = " <<
part.E()
279 <<
" and vtx and momentum in Global (not volTPC) coords are ";
294 mf::LogInfo(
"Track3DKalman: ") <<
" repMC, covMC are ... ";
303 while (trackIter != trackIn.
end()) {
305 spacepoints = fmsp.at(
nTrks);
307 mf::LogInfo(
"Track3DKalman: ") <<
"found " << spacepoints.size() <<
" 3D spacepoint(s).";
310 if (spacepoints.size() > 0) {
312 const double resolution = 0.5;
319 momM.SetX(gauss.fire(momM.X(), momErr.X() ));
320 momM.SetY(gauss.fire(momM.Y(), momErr.Y() ));
321 momM.SetZ(gauss.fire(momM.Z(), momErr.Z() ));
324 std::sort(spacepoints.begin(), spacepoints.end(),
sp_sort_3dz);
334 (TVector3)(spacepoints[0]->XYZ()),
347 for (
unsigned int point = 0; point < spacepoints.size(); ++point) {
349 TVector3 spt3(spacepoints[point]->XYZ());
352 TVector3 magNew(spt3[0], spt3[1], spt3[2]);
353 TVector3 magLast(spacepoints.back()->XYZ()[0],
354 spacepoints.back()->XYZ()[1],
355 spacepoints.back()->XYZ()[2]);
356 if (!(magNew.Mag() >= magLast.Mag() + eps || magNew.Mag() <= magLast.Mag() - eps))
375 <<
"ihit xyz..." << spt3[0] <<
"," << spt3[1] <<
"," << spt3[2];
393 mf::LogError(
"Track3DKalman: ") <<
"just caught a GFException." << std::endl;
396 <<
"Exceptions won't be further handled ->exit(1) " << __LINE__;
403 MF_LOG_DEBUG(
"Track3DKalman: ") << __FILE__ <<
" " << __LINE__;
404 MF_LOG_DEBUG(
"Track3DKalman: ") <<
"Track3DKalman.cxx: Original plane:";
407 MF_LOG_DEBUG(
"Track3DKalman: ") <<
"Current (fit) reference Plane:";
410 MF_LOG_DEBUG(
"Track3DKalman: ") <<
"Track3DKalman.cxx: Last reference Plane:";
415 MF_LOG_DEBUG(
"Track3DKalman: ") <<
"Track3DKalman: Original hit plane (not " 416 "surprisingly) not current reference Plane!" 425 MF_LOG_DEBUG(
"Track3DKalman: ") <<
" Final State and Cov:";
431 nfail = fitTrack.getFailedHits();
433 TVector3 dircoss = (*trackIter)->VertexDirection<TVector3>();
435 for (
int ii = 0; ii < 3; ++ii) {
436 fpMCT[ii] = MCMomentum[ii] / MCMomentum.Mag();
441 fpMCT[3] = MCMomentum.Mag();
442 fpREC[3] = -1.0 / (*stREC)[0][0];
446 <<
"Track3DKalman about to do tree->Fill(). Chi2/ndf is " <<
chi2 /
ndf 447 <<
". All in volTPC coords .... pMCT[0-3] is " <<
fpMCT[0] <<
", " <<
fpMCT[1] <<
", " 449 <<
", " <<
fpREC[2] <<
", " <<
fpREC[3] <<
".";
454 std::vector<art::Ptr<recob::Cluster>> clusters = fmc.at(
nTrks);
455 std::vector<art::Ptr<recob::Hit>>
hits = fmh.at(
nTrks);
458 std::vector<TVector3> dircos(spacepoints.size());
459 dircos[0] = TVector3(fpREC);
460 dircos.back() = TVector3(
fpRECL);
463 std::vector<TVector3> xyz(spacepoints.size());
464 size_t spStart = spcol->size();
465 for (
size_t tv = 0; tv < spacepoints.size(); ++tv) {
466 xyz[tv] = TVector3(spacepoints[tv]->XYZ());
467 spcol->push_back(*(spacepoints[tv].
get()));
469 size_t spEnd = spcol->size();
495 if (trackIter != trackIn.
end()) trackIter++;
501 evt.
put(std::move(tcol));
502 evt.
put(std::move(spcol));
503 evt.
put(std::move(tcassn));
504 evt.
put(std::move(thassn));
505 evt.
put(std::move(tspassn));
unsigned int getNDF() const
static bool sp_sort_3dz(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
const GFDetPlane & getReferencePlane() const
GFDetPlane getLastPlane() const
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
std::vector< double > fMomStart
TrackTrajectory::Flags_t Flags_t
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 >> SMatrixSym55
void Print(std::ostream &out=std::cout) const
const TMatrixT< Double_t > & getState() const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
std::vector< double > fMomErr
genf::GFAbsTrackRep * rep
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
TVector3 getNormal() const
typename data_t::const_iterator const_iterator
void push_back(Ptr< U > const &p)
A trajectory in space reconstructed from hits.
void setNumIterations(Int_t i)
Set number of iterations for Kalman Filter.
void processTrack(GFTrack *)
Performs fit on a GFTrack.
static GFFieldManager * getInstance()
std::string fGenieGenModuleLabel
genf::GFAbsTrackRep * repMC
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
const TMatrixT< Double_t > & getCov() const
const char * what() const noexcept override
standard error message handling for exceptions. use like "std::cerr << e.what();" ...
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.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
std::string fSpacePtsModuleLabel
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
TMatrixT< Double_t > * stMCT
TMatrixT< Double_t > * stREC
void init(GFAbsBField *b)
set the magntic field here. Magnetic field classes must be derived from GFAbsBField ...
EventNumber_t event() const
CLHEP::HepRandomEngine & fEngine
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
std::vector< double > fPosErr
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
TMatrixT< Double_t > * covMCT
TMatrixT< Double_t > * covREC