99 auto tcol = std::make_unique<std::vector<recob::Track>>();
100 auto spcol = std::make_unique<std::vector<recob::SpacePoint>>();
101 auto tspassn = std::make_unique<art::Assns<recob::Track, recob::SpacePoint>>();
102 auto tcassn = std::make_unique<art::Assns<recob::Track, recob::Cluster>>();
103 auto thassn = std::make_unique<art::Assns<recob::Track, recob::Hit>>();
104 auto shassn = std::make_unique<art::Assns<recob::SpacePoint, recob::Hit>>();
115 double timetick = 0.198;
117 const double wireShift =
119 double plane_pitch = geom->
PlanePitch(tpcid, 0, 1);
121 double Efield_drift = 0.5;
122 double Efield_SI = 0.7;
123 double Efield_IC = 0.9;
124 double Temperature = 90.;
126 double driftvelocity =
127 detProp.DriftVelocity(Efield_drift, Temperature);
128 double driftvelocity_SI = detProp.DriftVelocity(
129 Efield_SI, Temperature);
130 double driftvelocity_IC = detProp.DriftVelocity(
131 Efield_IC, Temperature);
132 double timepitch = driftvelocity * timetick;
133 double tSI = plane_pitch / driftvelocity_SI /
135 double tIC = plane_pitch / driftvelocity_IC /
148 for (
unsigned int i = 0; i < endpointListHandle->size(); ++i) {
155 std::vector<double> Iwirefirsts;
156 std::vector<double> Iwirelasts;
157 std::vector<double> Itimefirsts;
158 std::vector<double> Itimelasts;
159 std::vector<double> Itimefirsts_line;
160 std::vector<double> Itimelasts_line;
161 std::vector<std::vector<art::Ptr<recob::Hit>>> IclusHitlists;
162 std::vector<unsigned int> Icluster_count;
165 std::vector<double> Cwirefirsts;
166 std::vector<double> Cwirelasts;
167 std::vector<double> Ctimefirsts;
168 std::vector<double> Ctimelasts;
169 std::vector<double> Ctimefirsts_line;
170 std::vector<double> Ctimelasts_line;
171 std::vector<std::vector<art::Ptr<recob::Hit>>> CclusHitlists;
172 std::vector<unsigned int> Ccluster_count;
176 for (
unsigned int ii = 0; ii < clusterListHandle->size(); ++ii) {
182 int vtx2d_w = -99999;
183 double vtx2d_t = -99999;
184 bool found2dvtx =
false;
186 for (
unsigned int j = 0; j < endpointlist.
size(); j++) {
187 if (endpointlist[j]->View() == cl->View()) {
188 vtx2d_w = endpointlist[j]->WireID().Wire;
189 vtx2d_t = endpointlist[j]->DriftTime();
195 double w = cl->StartWire();
196 double t = cl->StartTick();
197 double dtdw = std::tan(cl->StartAngle());
198 double t_vtx = t + dtdw * (vtx2d_w -
w);
199 double dis =
std::abs(vtx2d_t - t_vtx);
207 std::vector<art::Ptr<recob::Hit>> hitlist = fm.at(ii);
210 TGraph* the2Dtrack =
new TGraph(hitlist.size());
212 std::vector<double> wires;
213 std::vector<double> times;
218 theHit != hitlist.end();
222 time = (*theHit)->PeakTime();
224 time -= presamplings;
230 wire_cm = (double)(((*theHit)->WireID().Wire + 3.95) * wire_pitch);
232 wire_cm = (double)(((*theHit)->WireID().Wire + 1.84) * wire_pitch);
237 time_cm = (double)((time - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
239 time_cm = time * driftvelocity_SI * timetick;
241 wires.push_back(wire_cm);
242 times.push_back(time_cm);
244 the2Dtrack->SetPoint(np, wire_cm, time_cm);
250 the2Dtrack->Fit(
"pol1",
"Q");
253 std::cout <<
"The 2D track fit failed" << std::endl;
257 TF1* pol1 = (TF1*)the2Dtrack->GetFunction(
"pol1");
259 pol1->GetParameters(par);
260 double intercept = par[0];
261 double slope = par[1];
263 double w0 = wires.front();
264 double w1 = wires.back();
265 double t0 = times.front();
266 double t1 = times.back();
267 double t0_line = intercept + (w0)*slope;
268 double t1_line = intercept + (w1)*slope;
271 switch (geom->
SignalType((*hitlist.begin())->Channel())) {
273 Iwirefirsts.push_back(w0);
274 Iwirelasts.push_back(w1);
275 Itimefirsts.push_back(t0);
276 Itimelasts.push_back(t1);
277 Itimefirsts_line.push_back(t0_line);
278 Itimelasts_line.push_back(t1_line);
279 IclusHitlists.push_back(hitlist);
280 Icluster_count.push_back(ii);
283 Cwirefirsts.push_back(w0);
284 Cwirelasts.push_back(w1);
285 Ctimefirsts.push_back(t0);
286 Ctimelasts.push_back(t1);
287 Ctimefirsts_line.push_back(t0_line);
288 Ctimelasts_line.push_back(t1_line);
289 CclusHitlists.push_back(hitlist);
290 Ccluster_count.push_back(ii);
301 for (
unsigned int collectionIter = 0; collectionIter < CclusHitlists.size();
304 double Cw0 = Cwirefirsts[collectionIter];
305 double Cw1 = Cwirelasts[collectionIter];
308 double Ct0_line = Ctimefirsts_line[collectionIter];
309 double Ct1_line = Ctimelasts_line[collectionIter];
310 std::vector<art::Ptr<recob::Hit>> hitsCtrk = CclusHitlists[collectionIter];
313 TMath::Sqrt(TMath::Power(Ct1_line - Ct0_line, 2) + TMath::Power(Cw1 - Cw0, 2));
315 for (
unsigned int inductionIter = 0; inductionIter < IclusHitlists.size();
318 double Iw0 = Iwirefirsts[inductionIter];
319 double Iw1 = Iwirelasts[inductionIter];
322 double It0_line = Itimefirsts_line[inductionIter];
323 double It1_line = Itimelasts_line[inductionIter];
324 std::vector<art::Ptr<recob::Hit>> hitsItrk = IclusHitlists[inductionIter];
327 TMath::Sqrt(TMath::Power(It1_line - It0_line, 2) + TMath::Power(Iw1 - Iw0, 2));
329 bool forward_match = ((
std::abs(Ct0_line - It0_line) <
ftmatch * timepitch) &&
331 bool backward_match = ((
std::abs(Ct0_line - It1_line) <
ftmatch * timepitch) &&
334 if (forward_match || backward_match) {
339 XYZ0.SetXYZ(Ct0_line,
340 (Cw0 - Iw0) / (2. * TMath::Sin(Angle)),
341 (Cw0 + Iw0) / (2. * TMath::Cos(Angle)) - YC / 2. * TMath::Tan(Angle));
342 XYZ1.SetXYZ(Ct1_line,
343 (Cw1 - Iw1) / (2. * TMath::Sin(Angle)),
344 (Cw1 + Iw1) / (2. * TMath::Cos(Angle)) - YC / 2. * TMath::Tan(Angle));
347 XYZ0.SetXYZ(Ct0_line,
348 (Cw0 - Iw1) / (2. * TMath::Sin(Angle)),
349 (Cw0 + Iw1) / (2. * TMath::Cos(Angle)) - YC / 2. * TMath::Tan(Angle));
350 XYZ1.SetXYZ(Ct1_line,
351 (Cw1 - Iw0) / (2. * TMath::Sin(Angle)),
352 (Cw1 + Iw0) / (2. * TMath::Cos(Angle)) - YC / 2. * TMath::Tan(Angle));
358 TVector3 startpointVec, endpointVec;
359 TVector2 collVtx, indVtx;
360 if (XYZ0.Z() <= XYZ1.Z()) {
361 startpointVec.SetXYZ(XYZ0.X(), XYZ0.Y(), XYZ0.Z());
362 endpointVec.SetXYZ(XYZ1.X(), XYZ1.Y(), XYZ1.Z());
364 collVtx.Set(Ct0_line, Cw0);
365 indVtx.Set(It0_line, Iw0);
368 collVtx.Set(Ct0_line, Cw0);
369 indVtx.Set(It1_line, Iw1);
373 startpointVec.SetXYZ(XYZ1.X(), XYZ1.Y(), XYZ1.Z());
374 endpointVec.SetXYZ(XYZ0.X(), XYZ0.Y(), XYZ0.Z());
376 collVtx.Set(Ct1_line, Cw1);
377 indVtx.Set(It1_line, Iw1);
380 collVtx.Set(Ct1_line, Cw1);
381 indVtx.Set(It0_line, Iw0);
386 TVector3 DirCos = endpointVec - startpointVec;
393 std::cout <<
"The Spacepoint is infinitely small" << std::endl;
408 std::vector<recob::SpacePoint> spacepoints;
410 std::vector<art::Ptr<recob::Hit>> minhits =
411 hitsCtrk.size() <= hitsItrk.size() ? hitsCtrk : hitsItrk;
412 std::vector<art::Ptr<recob::Hit>> maxhits =
413 hitsItrk.size() < hitsCtrk.size() ? hitsCtrk : hitsItrk;
415 std::vector<bool> maxhitsMatch(maxhits.size());
416 for (
unsigned int it = 0; it < maxhits.size(); it++)
417 maxhitsMatch[it] =
false;
419 std::vector<recob::Hit*> hits3Dmatched;
421 unsigned int imaximum = 0;
422 size_t spStart = spcol->size();
423 for (
unsigned int imin = 0; imin < minhits.size(); imin++) {
427 auto const hitSigType = minhits[imin]->SignalType();
432 w1 = (double)((hit1WireID.
Wire + 3.95) * wire_pitch);
434 w1 = (double)((hit1WireID.
Wire + 1.84) * wire_pitch);
436 double temptime1 = minhits[imin]->PeakTime() - presamplings;
441 t1 = (double)((temptime1 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
443 t1 = temptime1 * driftvelocity_SI * timetick;
448 minVtx2D.Set(indVtx.X(), indVtx.Y());
451 maxVtx2D.Set(collVtx.X(), collVtx.Y());
454 (collLength > indLength) ? collLength / indLength : indLength / collLength;
457 double minDistance = ratio * TMath::Sqrt(TMath::Power(t1 - minVtx2D.X(), 2) +
458 TMath::Power(w1 - minVtx2D.Y(), 2));
461 double difference = 9999999.;
463 for (
unsigned int imax = 0; imax < maxhits.size();
465 if (!maxhitsMatch[imax]) {
468 auto const hit2SigType = maxhits[imax]->SignalType();
471 w2 = (double)((hit2WireID.
Wire + 3.95) * wire_pitch);
473 w2 = (double)((hit2WireID.
Wire + 1.84) * wire_pitch);
475 double temptime2 = maxhits[imax]->PeakTime() - presamplings;
479 t2 = (double)((temptime2 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
481 t2 = temptime2 * driftvelocity_SI * timetick;
484 bool wirematch = (
std::abs(w1 - w2) < wireShift * wire_pitch);
486 double maxDistance = TMath::Sqrt(TMath::Power(t2 - maxVtx2D.X(), 2) +
487 TMath::Power(w2 - maxVtx2D.Y(), 2));
488 if (wirematch && timematch &&
std::abs(maxDistance - minDistance) < difference) {
489 difference =
std::abs(maxDistance - minDistance);
494 maxhitsMatch[imaximum] =
true;
497 if (difference != 9999999.) {
503 geo::WireID hit2WireID = maxhits[imaximum]->WireID();
504 auto const hit2SigType = maxhits[imaximum]->SignalType();
507 double w1_match = 0.;
509 w1_match = (double)((hit2WireID.
Wire + 3.95) * wire_pitch);
511 w1_match = (double)((hit2WireID.
Wire + 1.84) * wire_pitch);
513 double temptime3 = maxhits[imaximum]->PeakTime() - presamplings;
518 (double)((temptime3 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
520 t1_match = temptime3 * driftvelocity_SI * timetick;
527 const TVector3 hit3d(Ct,
528 (Cw - Iw) / (2. * TMath::Sin(Angle)),
529 (Cw + Iw) / (2. * TMath::Cos(Angle)) -
530 YC / 2. * TMath::Tan(Angle));
532 Double_t hitcoord[3];
533 hitcoord[0] = hit3d.X();
534 hitcoord[1] = hit3d.Y();
535 hitcoord[2] = hit3d.Z();
556 spStart + spacepoints.size());
558 const double eps(0.1);
559 if (spacepoints.size() >= 1) {
560 TVector3 magNew(mysp.XYZ()[0], mysp.XYZ()[1], mysp.XYZ()[2]);
561 TVector3 magLast(spacepoints.back().XYZ()[0],
562 spacepoints.back().XYZ()[1],
563 spacepoints.back().XYZ()[2]);
564 if (!(magNew.Mag() >= magLast.Mag() + eps || magNew.Mag() <= magLast.Mag() - eps))
567 spacepoints.push_back(mysp);
568 spcol->push_back(mysp);
573 size_t spEnd = spcol->size();
576 if (spacepoints.size() > 0) {
579 std::vector<TVector3> xyz(spacepoints.size());
580 for (
size_t s = 0; s < spacepoints.size(); ++s) {
581 xyz[s] = TVector3(spacepoints[s].XYZ());
585 std::vector<TVector3> dircos(spacepoints.size(), DirCos);
607 for (
size_t cpt = 0; cpt < clustersPerTrack.
size(); ++cpt)
616 mf::LogVerbatim(
"Summary") << std::setfill(
'-') << std::setw(175) <<
"-" << std::setfill(
' ');
618 for (
unsigned int i = 0; i < tcol->size(); ++i)
621 evt.
put(std::move(tcol));
622 evt.
put(std::move(spcol));
623 evt.
put(std::move(tspassn));
624 evt.
put(std::move(tcassn));
625 evt.
put(std::move(thassn));
626 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.
double fvertexclusterWindow
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
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.
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 >> SMatrixSym55
WireID_t Wire
Index of the wire within its plane.
std::string fClusterModuleLabel
double ThetaZ() const
Angle of the wires from positive z axis; .
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.
bool SortByWire(art::Ptr< recob::Hit > const &h1, art::Ptr< recob::Hit > const &h2)
void push_back(Ptr< U > const &p)
Signal from induction planes.
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.
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.
SigType_t SignalType(PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
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.
constexpr double kBogusD
obviously bogus double value
Length_t DetHalfHeight(TPCID const &tpcid=tpc_zero) const
Returns the half height of the active volume of the specified TPC.
std::string fEndPoint2DModuleLabel
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:
Signal from collection planes.
1.8.11 
