89 auto tcol = std::make_unique<std::vector<recob::Track>>();
90 auto spacepoints = std::make_unique<std::vector<recob::SpacePoint>>();
91 auto cassn = std::make_unique<art::Assns<recob::Track, recob::Cluster>>();
92 auto sassn = std::make_unique<art::Assns<recob::Track, recob::SpacePoint>>();
93 auto shassn = std::make_unique<art::Assns<recob::SpacePoint, recob::Hit>>();
94 auto hassn = std::make_unique<art::Assns<recob::Track, recob::Hit>>();
104 double timetick = 0.198;
105 double presamplings = 60.;
106 const double wireShift =
108 double plane_pitch = geom->
PlanePitch(tpcid, 0, 1);
110 double Efield_drift = 0.5;
111 double Efield_SI = 0.7;
112 double Efield_IC = 0.9;
113 double Temperature = 87.6;
115 double driftvelocity =
116 detProp.DriftVelocity(Efield_drift, Temperature);
118 double driftvelocity_SI =
119 detProp.DriftVelocity(Efield_SI, Temperature);
121 double driftvelocity_IC =
122 detProp.DriftVelocity(Efield_IC, Temperature);
124 double timepitch = driftvelocity * timetick;
125 double tSI = plane_pitch / driftvelocity_SI / timetick;
127 double tIC = plane_pitch / driftvelocity_IC / timetick;
136 std::vector<double> Iwirefirsts;
137 std::vector<double> Iwirelasts;
138 std::vector<double> Itimefirsts;
139 std::vector<double> Itimelasts;
140 std::vector<double> Itimefirsts_line;
141 std::vector<double> Itimelasts_line;
142 std::vector<std::vector<art::Ptr<recob::Hit>>> IclusHitlists;
143 std::vector<unsigned int> Icluster_count;
146 std::vector<double> Cwirefirsts;
147 std::vector<double> Cwirelasts;
148 std::vector<double> Ctimefirsts;
149 std::vector<double> Ctimelasts;
150 std::vector<double> Ctimefirsts_line;
151 std::vector<double> Ctimelasts_line;
152 std::vector<std::vector<art::Ptr<recob::Hit>>> CclusHitlists;
153 std::vector<unsigned int> Ccluster_count;
157 unsigned int wire = 0;
158 unsigned int plane = 0;
160 size_t startSPIndex =
167 for (
size_t ii = 0; ii < clusterListHandle->size(); ++ii) {
177 if (cl->View() ==
geo::kZ)
continue;
180 std::vector<art::Ptr<recob::Hit>> hitlist(fmh.at(ii));
182 if (hitlist.size() == 1)
187 std::sort(hitlist.begin(), hitlist.end());
189 TGraph* the2Dtrack =
new TGraph(hitlist.size());
191 std::vector<double> wires;
192 std::vector<double> times;
197 theHit != hitlist.end();
200 time = (*theHit)->PeakTime();
202 time -= presamplings;
204 plane = (*theHit)->WireID().Plane;
205 wire = (*theHit)->WireID().Wire;
208 if (plane == 1) time -= tIC;
213 wire_cm = (double)((wire + 3.95) * wire_pitch);
215 wire_cm = (double)((wire + 1.84) * wire_pitch);
219 time_cm = (double)((time - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
221 time_cm = time * driftvelocity_SI * timetick;
223 wires.push_back(wire_cm);
224 times.push_back(time_cm);
226 the2Dtrack->SetPoint(np, wire_cm, time_cm);
232 the2Dtrack->Fit(
"pol1",
"Q");
239 TF1* pol1 = (TF1*)the2Dtrack->GetFunction(
"pol1");
241 pol1->GetParameters(par);
242 double intercept = par[0];
243 double slope = par[1];
245 double w0 = wires.front();
246 double w1 = wires.back();
247 double t0 = times.front();
248 double t1 = times.back();
249 double t0_line = intercept + (w0)*slope;
250 double t1_line = intercept + (w1)*slope;
255 Iwirefirsts.push_back(w0);
256 Iwirelasts.push_back(w1);
257 Itimefirsts.push_back(t0);
258 Itimelasts.push_back(t1);
259 Itimefirsts_line.push_back(t0_line);
260 Itimelasts_line.push_back(t1_line);
261 IclusHitlists.push_back(hitlist);
262 Icluster_count.push_back(ii);
265 Cwirefirsts.push_back(w0);
266 Cwirelasts.push_back(w1);
267 Ctimefirsts.push_back(t0);
268 Ctimelasts.push_back(t1);
269 Ctimefirsts_line.push_back(t0_line);
270 Ctimelasts_line.push_back(t1_line);
271 CclusHitlists.push_back(hitlist);
272 Ccluster_count.push_back(ii);
283 for (
size_t collectionIter = 0; collectionIter < CclusHitlists.size(); ++collectionIter) {
285 double Cw0 = Cwirefirsts[collectionIter];
286 double Cw1 = Cwirelasts[collectionIter];
289 double Ct0_line = Ctimefirsts_line[collectionIter];
290 double Ct1_line = Ctimelasts_line[collectionIter];
291 std::vector<art::Ptr<recob::Hit>> hitsCtrk = CclusHitlists[collectionIter];
293 double collLength = std::hypot(Ct1_line - Ct0_line, Cw1 - Cw0);
296 for (
size_t inductionIter = 0; inductionIter < IclusHitlists.size(); ++inductionIter) {
298 double Iw0 = Iwirefirsts[inductionIter];
299 double Iw1 = Iwirelasts[inductionIter];
300 double It0_line = Itimefirsts_line[inductionIter];
301 double It1_line = Itimelasts_line[inductionIter];
302 std::vector<art::Ptr<recob::Hit>> hitsItrk = IclusHitlists[inductionIter];
304 double indLength = std::hypot(It1_line - It0_line, Iw1 - Iw0);
306 bool forward_match = ((
std::abs(Ct0_line - It0_line) <
ftmatch * timepitch) &&
308 bool backward_match = ((
std::abs(Ct0_line - It1_line) <
ftmatch * timepitch) &&
312 if (forward_match || backward_match) {
317 XYZ0.SetXYZ(Ct0_line,
318 (Cw0 - Iw0) / (2. * std::sin(Angle)),
319 (Cw0 + Iw0) / (2. * std::cos(Angle)) - YC / 2. * std::tan(Angle));
320 XYZ1.SetXYZ(Ct1_line,
321 (Cw1 - Iw1) / (2. * std::sin(Angle)),
322 (Cw1 + Iw1) / (2. * std::cos(Angle)) - YC / 2. * std::tan(Angle));
325 XYZ0.SetXYZ(Ct0_line,
326 (Cw0 - Iw1) / (2. * std::sin(Angle)),
327 (Cw0 + Iw1) / (2. * std::cos(Angle)) - YC / 2. * std::tan(Angle));
328 XYZ1.SetXYZ(Ct1_line,
329 (Cw1 - Iw0) / (2. * std::sin(Angle)),
330 (Cw1 + Iw0) / (2. * std::cos(Angle)) - YC / 2. * std::tan(Angle));
336 TVector3 startpointVec, endpointVec;
337 TVector2 collVtx, indVtx;
338 if (XYZ0.Z() <= XYZ1.Z()) {
339 startpointVec.SetXYZ(XYZ0.X(), XYZ0.Y(), XYZ0.Z());
340 endpointVec.SetXYZ(XYZ1.X(), XYZ1.Y(), XYZ1.Z());
342 collVtx.Set(Ct0_line, Cw0);
343 indVtx.Set(It0_line, Iw0);
346 collVtx.Set(Ct0_line, Cw0);
347 indVtx.Set(It1_line, Iw1);
351 startpointVec.SetXYZ(XYZ1.X(), XYZ1.Y(), XYZ1.Z());
352 endpointVec.SetXYZ(XYZ0.X(), XYZ0.Y(), XYZ0.Z());
354 collVtx.Set(Ct1_line, Cw1);
355 indVtx.Set(It1_line, Iw1);
358 collVtx.Set(Ct1_line, Cw1);
359 indVtx.Set(It0_line, Iw0);
364 TVector3 DirCos = endpointVec - startpointVec;
371 mf::LogWarning(
"Track3Dreco") <<
"The Spacepoint is infinitely small";
385 std::vector<art::Ptr<recob::Hit>> minhits =
386 hitsCtrk.size() <= hitsItrk.size() ? hitsCtrk : hitsItrk;
387 std::vector<art::Ptr<recob::Hit>> maxhits =
388 hitsItrk.size() <= hitsCtrk.size() ? hitsCtrk : hitsItrk;
390 std::vector<bool> maxhitsMatch(maxhits.size());
391 for (
size_t it = 0; it < maxhits.size(); ++it)
392 maxhitsMatch[it] =
false;
394 std::vector<recob::Hit*> hits3Dmatched;
396 unsigned int imaximum = 0;
399 startSPIndex = spacepoints->size();
401 for (
size_t imin = 0; imin < minhits.size(); ++imin) {
403 unsigned int wire, plane1, plane2;
404 wire = minhits[imin]->WireID().Wire;
405 plane1 = minhits[imin]->WireID().Plane;
410 w1 = (double)((wire + 3.95) * wire_pitch);
412 w1 = (double)((wire + 1.84) * wire_pitch);
413 double temptime1 = minhits[imin]->PeakTime() - presamplings;
414 if (plane1 == 1) temptime1 -= tIC;
417 t1 = (double)((temptime1 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
419 t1 = temptime1 * driftvelocity_SI * timetick;
423 plane1 == 1 ? minVtx2D.Set(collVtx.X(), collVtx.Y()) :
424 minVtx2D.Set(indVtx.X(), indVtx.Y());
426 plane1 == 1 ? maxVtx2D.Set(indVtx.X(), indVtx.Y()) :
427 maxVtx2D.Set(collVtx.X(), collVtx.Y());
430 (collLength > indLength) ? collLength / indLength : indLength / collLength;
433 double minDistance = ratio * TMath::Sqrt(TMath::Power(t1 - minVtx2D.X(), 2) +
434 TMath::Power(w1 - minVtx2D.Y(), 2));
437 double difference = 9999999.;
440 for (
size_t imax = 0; imax < maxhits.size(); ++imax) {
441 if (!maxhitsMatch[imax]) {
443 wire = maxhits[imax]->WireID().Wire;
444 plane2 = maxhits[imax]->WireID().Plane;
448 w2 = (double)((wire + 3.95) * wire_pitch);
450 w2 = (double)((wire + 1.84) * wire_pitch);
451 double temptime2 = maxhits[imax]->PeakTime() - presamplings;
452 if (plane2 == 1) temptime2 -= tIC;
455 t2 = (double)((temptime2 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
457 t2 = temptime2 * driftvelocity_SI * timetick;
460 bool wirematch = (
std::abs(w1 - w2) < wireShift * wire_pitch);
462 double maxDistance = TMath::Sqrt(TMath::Power(t2 - maxVtx2D.X(), 2) +
463 TMath::Power(w2 - maxVtx2D.Y(), 2));
464 if (wirematch && timematch &&
std::abs(maxDistance - minDistance) < difference) {
465 difference =
std::abs(maxDistance - minDistance);
470 maxhitsMatch[imaximum] =
true;
473 if (difference != 9999999.) {
479 wire = maxhits[imaximum]->WireID().Wire;
480 plane2 = maxhits[imaximum]->WireID().Plane;
484 w1_match = (double)((wire + 3.95) * wire_pitch);
486 w1_match = (double)((wire + 1.84) * wire_pitch);
487 double temptime3 = maxhits[imaximum]->PeakTime() - presamplings;
488 if (plane2 == 1) temptime3 -= tIC;
492 (double)((temptime3 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
494 t1_match = temptime3 * driftvelocity_SI * timetick;
497 double Ct = plane1 == 1 ? t1 : t1_match;
498 double Cw = plane1 == 1 ? w1 : w1_match;
499 double Iw = plane1 == 1 ? w1_match : w1;
501 const TVector3 hit3d(Ct,
502 (Cw - Iw) / (2. * std::sin(Angle)),
503 (Cw + Iw) / (2. * std::cos(Angle)) - YC / 2. * std::tan(Angle));
504 Double_t hitcoord[3];
505 hitcoord[0] = hit3d.X();
506 hitcoord[1] = hit3d.Y();
507 hitcoord[2] = hit3d.Z();
509 Double_t hitcoord_errs[3];
510 for (
int i = 0; i < 3; i++)
511 hitcoord_errs[i] = -1.000;
516 spacepoints->push_back(mysp);
523 endSPIndex = spacepoints->size();
526 if (spacepoints->size() > startSPIndex || clustersPerTrack.
size() > 0) {
528 std::vector<TVector3> xyz;
529 xyz.push_back(startpointVec);
530 xyz.push_back(endpointVec);
531 std::vector<TVector3> dir_xyz;
532 dir_xyz.push_back(DirCos);
533 dir_xyz.push_back(DirCos);
546 tcol->push_back(the3DTrack);
549 util::CreateAssn(evt, *tcol, *spacepoints, *sassn, startSPIndex, endSPIndex);
557 for (
size_t p = 0; p < clustersPerTrack.
size(); ++p) {
558 const std::vector<art::Ptr<recob::Hit>>&
hits = fmhc.at(p);
569 mf::LogVerbatim(
"Summary") << std::setfill(
'-') << std::setw(175) <<
"-" << std::setfill(
' ');
571 for (
unsigned int i = 0; i < tcol->size(); ++i) {
575 evt.
put(std::move(tcol));
576 evt.
put(std::move(spacepoints));
577 evt.
put(std::move(cassn));
578 evt.
put(std::move(sassn));
579 evt.
put(std::move(hassn));
580 evt.
put(std::move(shassn));
code to link reconstructed objects back to the MC truth information
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::string GetLArTPCVolumeName(TPCID const &tpcid=tpc_zero) const
Return the name of specified LAr TPC volume.
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
std::string fClusterModuleLabel
label for input cluster collection
The data type to uniquely identify a Plane.
TrackTrajectory::Flags_t Flags_t
constexpr auto abs(T v)
Returns the absolute value of the argument.
Planes which measure Z direction.
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 >> SMatrixSym55
double ThetaZ() const
Angle of the wires from positive z axis; .
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
typename data_t::const_iterator const_iterator
void push_back(Ptr< U > const &p)
A trajectory in space reconstructed from hits.
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
The data type to uniquely identify a TPC.
int ftmatch
tolerance for time matching (in time samples)
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
Length_t PlanePitch(TPCID const &tpcid, PlaneID::PlaneID_t p1=0, PlaneID::PlaneID_t p2=1) const
Returns the distance between two planes.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Length_t DetHalfHeight(TPCID const &tpcid=tpc_zero) const
Returns the half height of the active volume of the specified TPC.
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
1.8.11 
