93 art::ServiceHandle<geo::Geometry const> geom;
95 art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(
evt);
101 auto tcol = std::make_unique<std::vector<recob::Track>>();
102 auto spcol = std::make_unique<std::vector<recob::SpacePoint>>();
103 auto tspassn = std::make_unique<art::Assns<recob::Track, recob::SpacePoint>>();
104 auto tcassn = std::make_unique<art::Assns<recob::Track, recob::Cluster>>();
105 auto thassn = std::make_unique<art::Assns<recob::Track, recob::Hit>>();
106 auto shassn = std::make_unique<art::Assns<recob::SpacePoint, recob::Hit>>();
109 TString tpcName = geom->GetLArTPCVolumeName();
112 double YC = (geom->DetHalfHeight()) * 2.;
113 double Angle = geom->Plane(1).Wire(0).ThetaZ(
false) -
116 double timetick = 0.198;
118 const double wireShift =
120 double plane_pitch = geom->PlanePitch(0, 1);
121 double wire_pitch = geom->WirePitch();
122 double Efield_drift = 0.5;
123 double Efield_SI = 0.7;
124 double Efield_IC = 0.9;
125 double Temperature = 90.;
127 double driftvelocity =
128 detProp.DriftVelocity(Efield_drift, Temperature);
129 double driftvelocity_SI =
detProp.DriftVelocity(
130 Efield_SI, Temperature);
131 double driftvelocity_IC =
detProp.DriftVelocity(
132 Efield_IC, Temperature);
133 double timepitch = driftvelocity * timetick;
134 double tSI = plane_pitch / driftvelocity_SI /
136 double tIC = plane_pitch / driftvelocity_IC /
140 art::Handle<std::vector<recob::Cluster>> clusterListHandle;
144 art::Handle<std::vector<recob::EndPoint2D>> endpointListHandle;
147 art::PtrVector<recob::EndPoint2D> endpointlist;
149 for (
unsigned int i = 0; i < endpointListHandle->size(); ++i) {
150 art::Ptr<recob::EndPoint2D> endpointHolder(endpointListHandle, i);
151 endpointlist.push_back(endpointHolder);
156 std::vector<double> Iwirefirsts;
157 std::vector<double> Iwirelasts;
158 std::vector<double> Itimefirsts;
159 std::vector<double> Itimelasts;
160 std::vector<double> Itimefirsts_line;
161 std::vector<double> Itimelasts_line;
162 std::vector<std::vector<art::Ptr<recob::Hit>>> IclusHitlists;
163 std::vector<unsigned int> Icluster_count;
166 std::vector<double> Cwirefirsts;
167 std::vector<double> Cwirelasts;
168 std::vector<double> Ctimefirsts;
169 std::vector<double> Ctimelasts;
170 std::vector<double> Ctimefirsts_line;
171 std::vector<double> Ctimelasts_line;
172 std::vector<std::vector<art::Ptr<recob::Hit>>> CclusHitlists;
173 std::vector<unsigned int> Ccluster_count;
177 for (
unsigned int ii = 0; ii < clusterListHandle->size(); ++ii) {
178 art::Ptr<recob::Cluster> cl(clusterListHandle, ii);
183 int vtx2d_w = -99999;
184 double vtx2d_t = -99999;
185 bool found2dvtx =
false;
187 for (
unsigned int j = 0; j < endpointlist.size(); j++) {
188 if (endpointlist[j]->View() == cl->View()) {
189 vtx2d_w = endpointlist[j]->WireID().Wire;
190 vtx2d_t = endpointlist[j]->DriftTime();
196 double w = cl->StartWire();
197 double t = cl->StartTick();
198 double dtdw = std::tan(cl->StartAngle());
199 double t_vtx = t + dtdw * (vtx2d_w -
w);
200 double dis =
std::abs(vtx2d_t - t_vtx);
208 std::vector<art::Ptr<recob::Hit>> hitlist = fm.at(ii);
211 TGraph* the2Dtrack =
new TGraph(hitlist.size());
213 std::vector<double> wires;
214 std::vector<double> times;
218 for (
std::vector<art::Ptr<recob::Hit>>::const_iterator theHit = hitlist.begin();
219 theHit != hitlist.end();
223 time = (*theHit)->PeakTime();
225 time -= presamplings;
227 if (geom->SignalType((*theHit)->Channel()) ==
geo::kCollection) time -= tIC;
231 wire_cm = (double)(((*theHit)->WireID().Wire + 3.95) * wire_pitch);
233 wire_cm = (double)(((*theHit)->WireID().Wire + 1.84) * wire_pitch);
238 time_cm = (double)((time - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
240 time_cm = time * driftvelocity_SI * timetick;
242 wires.push_back(wire_cm);
243 times.push_back(time_cm);
245 the2Dtrack->SetPoint(np, wire_cm, time_cm);
251 the2Dtrack->Fit(
"pol1",
"Q");
254 std::cout <<
"The 2D track fit failed" << std::endl;
258 TF1* pol1 = (TF1*)the2Dtrack->GetFunction(
"pol1");
260 pol1->GetParameters(par);
261 double intercept = par[0];
262 double slope = par[1];
264 double w0 = wires.front();
265 double w1 = wires.back();
266 double t0 = times.front();
267 double t1 = times.back();
268 double t0_line = intercept + (w0)*slope;
269 double t1_line = intercept + (w1)*slope;
272 switch (geom->SignalType((*hitlist.begin())->Channel())) {
274 Iwirefirsts.push_back(w0);
275 Iwirelasts.push_back(w1);
276 Itimefirsts.push_back(t0);
277 Itimelasts.push_back(t1);
278 Itimefirsts_line.push_back(t0_line);
279 Itimelasts_line.push_back(t1_line);
280 IclusHitlists.push_back(hitlist);
281 Icluster_count.push_back(ii);
284 Cwirefirsts.push_back(w0);
285 Cwirelasts.push_back(w1);
286 Ctimefirsts.push_back(t0);
287 Ctimelasts.push_back(t1);
288 Ctimefirsts_line.push_back(t0_line);
289 Ctimelasts_line.push_back(t1_line);
290 CclusHitlists.push_back(hitlist);
291 Ccluster_count.push_back(ii);
302 for (
unsigned int collectionIter = 0; collectionIter < CclusHitlists.size();
305 double Cw0 = Cwirefirsts[collectionIter];
306 double Cw1 = Cwirelasts[collectionIter];
309 double Ct0_line = Ctimefirsts_line[collectionIter];
310 double Ct1_line = Ctimelasts_line[collectionIter];
311 std::vector<art::Ptr<recob::Hit>> hitsCtrk = CclusHitlists[collectionIter];
314 TMath::Sqrt(TMath::Power(Ct1_line - Ct0_line, 2) + TMath::Power(Cw1 - Cw0, 2));
316 for (
unsigned int inductionIter = 0; inductionIter < IclusHitlists.size();
319 double Iw0 = Iwirefirsts[inductionIter];
320 double Iw1 = Iwirelasts[inductionIter];
323 double It0_line = Itimefirsts_line[inductionIter];
324 double It1_line = Itimelasts_line[inductionIter];
325 std::vector<art::Ptr<recob::Hit>> hitsItrk = IclusHitlists[inductionIter];
328 TMath::Sqrt(TMath::Power(It1_line - It0_line, 2) + TMath::Power(Iw1 - Iw0, 2));
330 bool forward_match = ((
std::abs(Ct0_line - It0_line) <
ftmatch * timepitch) &&
332 bool backward_match = ((
std::abs(Ct0_line - It1_line) <
ftmatch * timepitch) &&
335 if (forward_match || backward_match) {
340 XYZ0.SetXYZ(Ct0_line,
341 (Cw0 - Iw0) / (2. * TMath::Sin(Angle)),
342 (Cw0 + Iw0) / (2. * TMath::Cos(Angle)) - YC / 2. * TMath::Tan(Angle));
343 XYZ1.SetXYZ(Ct1_line,
344 (Cw1 - Iw1) / (2. * TMath::Sin(Angle)),
345 (Cw1 + Iw1) / (2. * TMath::Cos(Angle)) - YC / 2. * TMath::Tan(Angle));
348 XYZ0.SetXYZ(Ct0_line,
349 (Cw0 - Iw1) / (2. * TMath::Sin(Angle)),
350 (Cw0 + Iw1) / (2. * TMath::Cos(Angle)) - YC / 2. * TMath::Tan(Angle));
351 XYZ1.SetXYZ(Ct1_line,
352 (Cw1 - Iw0) / (2. * TMath::Sin(Angle)),
353 (Cw1 + Iw0) / (2. * TMath::Cos(Angle)) - YC / 2. * TMath::Tan(Angle));
359 TVector3 startpointVec, endpointVec;
360 TVector2 collVtx, indVtx;
361 if (XYZ0.Z() <= XYZ1.Z()) {
362 startpointVec.SetXYZ(XYZ0.X(), XYZ0.Y(), XYZ0.Z());
363 endpointVec.SetXYZ(XYZ1.X(), XYZ1.Y(), XYZ1.Z());
365 collVtx.Set(Ct0_line, Cw0);
366 indVtx.Set(It0_line, Iw0);
369 collVtx.Set(Ct0_line, Cw0);
370 indVtx.Set(It1_line, Iw1);
374 startpointVec.SetXYZ(XYZ1.X(), XYZ1.Y(), XYZ1.Z());
375 endpointVec.SetXYZ(XYZ0.X(), XYZ0.Y(), XYZ0.Z());
377 collVtx.Set(Ct1_line, Cw1);
378 indVtx.Set(It1_line, Iw1);
381 collVtx.Set(Ct1_line, Cw1);
382 indVtx.Set(It0_line, Iw0);
387 TVector3 DirCos = endpointVec - startpointVec;
394 std::cout <<
"The Spacepoint is infinitely small" << std::endl;
398 art::Ptr<recob::Cluster> cl1(clusterListHandle, Icluster_count[inductionIter]);
399 art::Ptr<recob::Cluster> cl2(clusterListHandle, Ccluster_count[collectionIter]);
400 art::PtrVector<recob::Cluster> clustersPerTrack;
401 clustersPerTrack.push_back(cl1);
402 clustersPerTrack.push_back(cl2);
409 std::vector<recob::SpacePoint> spacepoints;
411 std::vector<art::Ptr<recob::Hit>> minhits =
412 hitsCtrk.size() <= hitsItrk.size() ? hitsCtrk : hitsItrk;
413 std::vector<art::Ptr<recob::Hit>> maxhits =
414 hitsItrk.size() < hitsCtrk.size() ? hitsCtrk : hitsItrk;
416 std::vector<bool> maxhitsMatch(maxhits.size());
417 for (
unsigned int it = 0; it < maxhits.size(); it++)
418 maxhitsMatch[it] =
false;
420 std::vector<recob::Hit*> hits3Dmatched;
422 unsigned int imaximum = 0;
423 size_t spStart = spcol->size();
424 for (
unsigned int imin = 0; imin < minhits.size(); imin++) {
428 auto const hitSigType = minhits[imin]->SignalType();
433 w1 = (double)((hit1WireID.
Wire + 3.95) * wire_pitch);
435 w1 = (double)((hit1WireID.
Wire + 1.84) * wire_pitch);
437 double temptime1 = minhits[imin]->PeakTime() - presamplings;
442 t1 = (double)((temptime1 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
444 t1 = temptime1 * driftvelocity_SI * timetick;
449 minVtx2D.Set(indVtx.X(), indVtx.Y());
452 maxVtx2D.Set(collVtx.X(), collVtx.Y());
455 (collLength > indLength) ? collLength / indLength : indLength / collLength;
458 double minDistance = ratio *
TMath::Sqrt(TMath::Power(t1 - minVtx2D.X(), 2) +
459 TMath::Power(w1 - minVtx2D.Y(), 2));
462 double difference = 9999999.;
464 for (
unsigned int imax = 0; imax < maxhits.size();
466 if (!maxhitsMatch[imax]) {
469 auto const hit2SigType = maxhits[imax]->SignalType();
472 w2 = (double)((hit2WireID.
Wire + 3.95) * wire_pitch);
474 w2 = (double)((hit2WireID.
Wire + 1.84) * wire_pitch);
476 double temptime2 = maxhits[imax]->PeakTime() - presamplings;
480 t2 = (double)((temptime2 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
482 t2 = temptime2 * driftvelocity_SI * timetick;
485 bool wirematch = (
std::abs(w1 - w2) < wireShift * wire_pitch);
487 double maxDistance =
TMath::Sqrt(TMath::Power(t2 - maxVtx2D.X(), 2) +
488 TMath::Power(w2 - maxVtx2D.Y(), 2));
489 if (wirematch && timematch &&
std::abs(maxDistance - minDistance) < difference) {
490 difference =
std::abs(maxDistance - minDistance);
495 maxhitsMatch[imaximum] =
true;
497 art::PtrVector<recob::Hit> sp_hits;
498 if (difference != 9999999.) {
499 sp_hits.push_back(minhits[imin]);
500 sp_hits.push_back(maxhits[imaximum]);
504 geo::WireID hit2WireID = maxhits[imaximum]->WireID();
505 auto const hit2SigType = maxhits[imaximum]->SignalType();
508 double w1_match = 0.;
510 w1_match = (double)((hit2WireID.
Wire + 3.95) * wire_pitch);
512 w1_match = (double)((hit2WireID.
Wire + 1.84) * wire_pitch);
514 double temptime3 = maxhits[imaximum]->PeakTime() - presamplings;
519 (double)((temptime3 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
521 t1_match = temptime3 * driftvelocity_SI * timetick;
528 const TVector3 hit3d(Ct,
529 (Cw - Iw) / (2. * TMath::Sin(Angle)),
530 (Cw + Iw) / (2. * TMath::Cos(Angle)) -
531 YC / 2. * TMath::Tan(Angle));
533 Double_t hitcoord[3];
534 hitcoord[0] = hit3d.X();
535 hitcoord[1] = hit3d.Y();
536 hitcoord[2] = hit3d.Z();
557 spStart + spacepoints.size());
559 const double eps(0.1);
560 if (spacepoints.size() >= 1) {
561 TVector3 magNew(mysp.XYZ()[0], mysp.XYZ()[1], mysp.XYZ()[2]);
562 TVector3 magLast(spacepoints.back().XYZ()[0],
563 spacepoints.back().XYZ()[1],
564 spacepoints.back().XYZ()[2]);
565 if (!(magNew.Mag() >= magLast.Mag() + eps || magNew.Mag() <= magLast.Mag() - eps))
568 spacepoints.push_back(mysp);
569 spcol->push_back(mysp);
574 size_t spEnd = spcol->size();
577 if (spacepoints.size() > 0) {
580 std::vector<TVector3> xyz(spacepoints.size());
581 for (
size_t s = 0;
s < spacepoints.size(); ++
s) {
582 xyz[
s] = TVector3(spacepoints[
s].XYZ());
586 std::vector<TVector3> dircos(spacepoints.size(), DirCos);
608 for (
size_t cpt = 0; cpt < clustersPerTrack.size(); ++cpt)
617 mf::LogVerbatim(
"Summary") << std::setfill(
'-') << std::setw(175) <<
"-" << std::setfill(
' ');
618 mf::LogVerbatim(
"Summary") <<
"SpacePts Summary:";
619 for (
unsigned int i = 0; i < tcol->size(); ++i)
620 mf::LogVerbatim(
"Summary") << tcol->at(i);
622 evt.put(std::move(tcol));
623 evt.put(std::move(spcol));
624 evt.put(std::move(tspassn));
625 evt.put(std::move(tcassn));
626 evt.put(std::move(thassn));
627 evt.put(std::move(shassn));
double fvertexclusterWindow
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 > > SMatrixSym55
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
EResult err(const char *call)
TrackTrajectory::Flags_t Flags_t
WireID_t Wire
Index of the wire within its plane.
std::string fClusterModuleLabel
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
bool SortByWire(art::Ptr< recob::Hit > const &h1, art::Ptr< recob::Hit > const &h2)
Signal from induction planes.
A trajectory in space reconstructed from hits.
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
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.
Var Sqrt(const Var &v)
Use to take sqrt of a var.
then echo File list $list not found else cat $list while read file do echo $file sed s
constexpr double kBogusD
obviously bogus double value
std::string fEndPoint2DModuleLabel
BEGIN_PROLOG could also be cout
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.