76 ftmatch = pset.get<
int>(
"TMatch");
81 produces<std::vector<recob::Track>>();
82 produces<std::vector<recob::SpacePoint>>();
83 produces<art::Assns<recob::Track, recob::SpacePoint>>();
84 produces<art::Assns<recob::Track, recob::Cluster>>();
85 produces<art::Assns<recob::Track, recob::Hit>>();
86 produces<art::Assns<recob::SpacePoint, recob::Hit>>();
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>>();
110 auto const& tpc = geom->
TPC(tpcid);
114 double YC = tpc.HalfHeight() * 2.;
115 double Angle = wireReadoutGeom.Plane(
geo::PlaneID{tpcid, 1}).Wire(0).ThetaZ(
false) -
118 double timetick = 0.198;
120 const double wireShift =
122 double plane_pitch = wireReadoutGeom.PlanePitch(tpcid);
123 double wire_pitch = wireReadoutGeom.FirstPlane(tpcid).WirePitch();
124 double Efield_drift = 0.5;
125 double Efield_SI = 0.7;
126 double Efield_IC = 0.9;
127 double Temperature = 90.;
129 double driftvelocity =
130 detProp.DriftVelocity(Efield_drift, Temperature);
131 double driftvelocity_SI = detProp.DriftVelocity(
132 Efield_SI, Temperature);
133 double driftvelocity_IC = detProp.DriftVelocity(
134 Efield_IC, Temperature);
135 double timepitch = driftvelocity * timetick;
136 double tSI = plane_pitch / driftvelocity_SI /
138 double tIC = plane_pitch / driftvelocity_IC /
151 for (
unsigned int i = 0; i < endpointListHandle->size(); ++i) {
158 std::vector<double> Iwirefirsts;
159 std::vector<double> Iwirelasts;
160 std::vector<double> Itimefirsts;
161 std::vector<double> Itimelasts;
162 std::vector<double> Itimefirsts_line;
163 std::vector<double> Itimelasts_line;
164 std::vector<std::vector<art::Ptr<recob::Hit>>> IclusHitlists;
165 std::vector<unsigned int> Icluster_count;
168 std::vector<double> Cwirefirsts;
169 std::vector<double> Cwirelasts;
170 std::vector<double> Ctimefirsts;
171 std::vector<double> Ctimelasts;
172 std::vector<double> Ctimefirsts_line;
173 std::vector<double> Ctimelasts_line;
174 std::vector<std::vector<art::Ptr<recob::Hit>>> CclusHitlists;
175 std::vector<unsigned int> Ccluster_count;
179 for (
unsigned int ii = 0; ii < clusterListHandle->size(); ++ii) {
185 int vtx2d_w = -99999;
186 double vtx2d_t = -99999;
187 bool found2dvtx =
false;
189 for (
unsigned int j = 0; j < endpointlist.
size(); j++) {
190 if (endpointlist[j]->View() == cl->
View()) {
191 vtx2d_w = endpointlist[j]->WireID().Wire;
192 vtx2d_t = endpointlist[j]->DriftTime();
201 double t_vtx = t + dtdw * (vtx2d_w -
w);
202 double dis =
std::abs(vtx2d_t - t_vtx);
210 std::vector<art::Ptr<recob::Hit>> hitlist = fm.at(ii);
213 TGraph* the2Dtrack =
new TGraph(hitlist.size());
215 std::vector<double> wires;
216 std::vector<double> times;
221 theHit != hitlist.end();
225 time = (*theHit)->PeakTime();
227 time -= presamplings;
233 if (wireReadoutGeom.SignalType((*theHit)->Channel()) ==
geo::kInduction)
234 wire_cm = (double)(((*theHit)->WireID().Wire + 3.95) * wire_pitch);
236 wire_cm = (double)(((*theHit)->WireID().Wire + 1.84) * wire_pitch);
241 time_cm = (double)((time - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
243 time_cm = time * driftvelocity_SI * timetick;
245 wires.push_back(wire_cm);
246 times.push_back(time_cm);
248 the2Dtrack->SetPoint(np, wire_cm, time_cm);
254 the2Dtrack->Fit(
"pol1",
"Q");
257 std::cout <<
"The 2D track fit failed" << std::endl;
261 TF1* pol1 = (TF1*)the2Dtrack->GetFunction(
"pol1");
263 pol1->GetParameters(par);
264 double intercept = par[0];
265 double slope = par[1];
267 double w0 = wires.front();
268 double w1 = wires.back();
269 double t0 = times.front();
270 double t1 = times.back();
271 double t0_line = intercept + (w0)*slope;
272 double t1_line = intercept + (w1)*slope;
275 switch (wireReadoutGeom.SignalType((*hitlist.begin())->Channel())) {
277 Iwirefirsts.push_back(w0);
278 Iwirelasts.push_back(w1);
279 Itimefirsts.push_back(t0);
280 Itimelasts.push_back(t1);
281 Itimefirsts_line.push_back(t0_line);
282 Itimelasts_line.push_back(t1_line);
283 IclusHitlists.push_back(hitlist);
284 Icluster_count.push_back(ii);
287 Cwirefirsts.push_back(w0);
288 Cwirelasts.push_back(w1);
289 Ctimefirsts.push_back(t0);
290 Ctimelasts.push_back(t1);
291 Ctimefirsts_line.push_back(t0_line);
292 Ctimelasts_line.push_back(t1_line);
293 CclusHitlists.push_back(hitlist);
294 Ccluster_count.push_back(ii);
305 for (
unsigned int collectionIter = 0; collectionIter < CclusHitlists.size();
308 double Cw0 = Cwirefirsts[collectionIter];
309 double Cw1 = Cwirelasts[collectionIter];
310 double Ct0_line = Ctimefirsts_line[collectionIter];
311 double Ct1_line = Ctimelasts_line[collectionIter];
312 std::vector<art::Ptr<recob::Hit>> hitsCtrk = CclusHitlists[collectionIter];
315 TMath::Sqrt(TMath::Power(Ct1_line - Ct0_line, 2) + TMath::Power(Cw1 - Cw0, 2));
317 for (
unsigned int inductionIter = 0; inductionIter < IclusHitlists.size();
320 double Iw0 = Iwirefirsts[inductionIter];
321 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++) {
426 auto const hitSigType = minhits[imin]->SignalType();
431 w1 = (double)((hit1WireID.
Wire + 3.95) * wire_pitch);
433 w1 = (double)((hit1WireID.
Wire + 1.84) * wire_pitch);
435 double temptime1 = minhits[imin]->PeakTime() - presamplings;
439 t1 = (double)((temptime1 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
441 t1 = temptime1 * driftvelocity_SI * timetick;
446 minVtx2D.Set(indVtx.X(), indVtx.Y());
449 maxVtx2D.Set(collVtx.X(), collVtx.Y());
452 (collLength > indLength) ? collLength / indLength : indLength / collLength;
455 double minDistance = ratio * TMath::Sqrt(TMath::Power(t1 - minVtx2D.X(), 2) +
456 TMath::Power(w1 - minVtx2D.Y(), 2));
459 double difference = 9999999.;
461 for (
unsigned int imax = 0; imax < maxhits.size();
463 if (!maxhitsMatch[imax]) {
466 auto const hit2SigType = maxhits[imax]->SignalType();
469 w2 = (double)((hit2WireID.
Wire + 3.95) * wire_pitch);
471 w2 = (double)((hit2WireID.
Wire + 1.84) * wire_pitch);
473 double temptime2 = maxhits[imax]->PeakTime() - presamplings;
477 t2 = (double)((temptime2 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
479 t2 = temptime2 * driftvelocity_SI * timetick;
482 bool wirematch = (
std::abs(w1 - w2) < wireShift * wire_pitch);
484 double maxDistance = TMath::Sqrt(TMath::Power(t2 - maxVtx2D.X(), 2) +
485 TMath::Power(w2 - maxVtx2D.Y(), 2));
486 if (wirematch && timematch &&
std::abs(maxDistance - minDistance) < difference) {
487 difference =
std::abs(maxDistance - minDistance);
492 maxhitsMatch[imaximum] =
true;
495 if (difference != 9999999.) {
501 geo::WireID hit2WireID = maxhits[imaximum]->WireID();
502 auto const hit2SigType = maxhits[imaximum]->SignalType();
504 double w1_match = 0.;
506 w1_match = (double)((hit2WireID.
Wire + 3.95) * wire_pitch);
508 w1_match = (double)((hit2WireID.
Wire + 1.84) * wire_pitch);
510 double temptime3 = maxhits[imaximum]->PeakTime() - presamplings;
515 (double)((temptime3 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
517 t1_match = temptime3 * driftvelocity_SI * timetick;
524 const TVector3 hit3d(Ct,
525 (Cw - Iw) / (2. * TMath::Sin(Angle)),
526 (Cw + Iw) / (2. * TMath::Cos(Angle)) -
527 YC / 2. * TMath::Tan(Angle));
529 Double_t hitcoord[3];
530 hitcoord[0] = hit3d.X();
531 hitcoord[1] = hit3d.Y();
532 hitcoord[2] = hit3d.Z();
538 spStart + spacepoints.size());
540 const double eps(0.1);
541 if (spacepoints.size() >= 1) {
542 TVector3 magNew(mysp.XYZ()[0], mysp.XYZ()[1], mysp.XYZ()[2]);
543 TVector3 magLast(spacepoints.back().XYZ()[0],
544 spacepoints.back().XYZ()[1],
545 spacepoints.back().XYZ()[2]);
546 if (!(magNew.Mag() >= magLast.Mag() + eps || magNew.Mag() <= magLast.Mag() - eps))
549 spacepoints.push_back(mysp);
550 spcol->push_back(mysp);
555 size_t spEnd = spcol->size();
558 if (spacepoints.size() > 0) {
561 std::vector<TVector3> xyz(spacepoints.size());
562 for (
size_t s = 0; s < spacepoints.size(); ++s) {
563 xyz[s] = TVector3(spacepoints[s].XYZ());
567 std::vector<TVector3> dircos(spacepoints.size(), DirCos);
589 for (
size_t cpt = 0; cpt < clustersPerTrack.
size(); ++cpt)
598 mf::LogVerbatim(
"Summary") << std::setfill(
'-') << std::setw(175) <<
"-" << std::setfill(
' ');
600 for (
unsigned int i = 0; i < tcol->size(); ++i)
603 evt.
put(std::move(tcol));
604 evt.
put(std::move(spcol));
605 evt.
put(std::move(tspassn));
606 evt.
put(std::move(tcassn));
607 evt.
put(std::move(thassn));
608 evt.
put(std::move(shassn));
code to link reconstructed objects back to the MC truth information
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
double fvertexclusterWindow
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
The data type to uniquely identify a Plane.
SpacePts(fhicl::ParameterSet const &pset)
TrackTrajectory::Flags_t Flags_t
constexpr auto abs(T v)
Returns the absolute value of the argument.
float StartWire() const
Returns the wire coordinate of the start of the cluster.
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 >> SMatrixSym55
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
WireID_t Wire
Index of the wire within its plane.
float StartAngle() const
Returns the starting angle of the cluster.
std::string fClusterModuleLabel
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
void produce(art::Event &evt)
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)
#define DEFINE_ART_MODULE(klass)
void push_back(Ptr< U > const &p)
Signal from induction planes.
A trajectory in space reconstructed from hits.
const TGeoVolume * ActiveVolume() const
Returns the expected drift direction based on geometry.
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
The data type to uniquely identify a TPC.
Declaration of cluster object.
Encapsulate the geometry of a wire .
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 operator()(art::Ptr< recob::Hit > const &h1, art::Ptr< recob::Hit > const &h2) const
geo::View_t View() const
Returns the view for this cluster.
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 .
Provides recob::Track data product.
constexpr double kBogusD
obviously bogus double value
std::string fEndPoint2DModuleLabel
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
float StartTick() const
Returns the tick coordinate of the start of the cluster.
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
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:
Signal from collection planes.