124 std::unique_ptr<std::vector<recob::Track> > tcol(
new std::vector<recob::Track>);
125 std::unique_ptr<std::vector<recob::SpacePoint> > spacepoints(
new std::vector<recob::SpacePoint>);
138 double timetick = 0.198;
139 double presamplings = 60.;
140 const double wireShift=50.;
143 double Efield_drift = 0.5;
144 double Efield_SI = 0.7;
145 double Efield_IC = 0.9;
146 double Temperature = 87.6;
148 double driftvelocity = detprop->
DriftVelocity(Efield_drift,Temperature);
150 double driftvelocity_SI = detprop->
DriftVelocity(Efield_SI,Temperature);
152 double driftvelocity_IC = detprop->
DriftVelocity(Efield_IC,Temperature);
154 double timepitch = driftvelocity*timetick;
155 double tSI = plane_pitch/driftvelocity_SI/timetick;
157 double tIC = plane_pitch/driftvelocity_IC/timetick;
168 std::vector<double> Iwirefirsts;
169 std::vector<double> Iwirelasts;
170 std::vector<double> Itimefirsts;
171 std::vector<double> Itimelasts;
172 std::vector<double> Itimefirsts_line;
173 std::vector<double> Itimelasts_line;
174 std::vector < std::vector<art::Ptr<recob::Hit> > > IclusHitlists;
175 std::vector<unsigned int> Icluster_count;
178 std::vector<double> Cwirefirsts;
179 std::vector<double> Cwirelasts;
180 std::vector<double> Ctimefirsts;
181 std::vector<double> Ctimelasts;
182 std::vector<double> Ctimefirsts_line;
183 std::vector<double> Ctimelasts_line;
184 std::vector< std::vector< art::Ptr<recob::Hit> > > CclusHitlists;
185 std::vector<unsigned int> Ccluster_count;
189 unsigned int wire = 0;
190 unsigned int plane = 0;
192 size_t startSPIndex = spacepoints->size();
193 size_t endSPIndex = spacepoints->size();
197 for(
size_t ii = 0; ii < clusterListHandle->size(); ++ii){
208 if (cl->View() ==
geo::kZ)
continue;
211 std::vector< art::Ptr<recob::Hit> > hitlist (fmh.at(ii));
213 if(hitlist.size() == 1)
continue;
217 std::sort(hitlist.begin(), hitlist.end());
219 TGraph *the2Dtrack =
new TGraph(hitlist.size());
221 std::vector<double> wires;
222 std::vector<double> times;
229 time = (*theHit)->PeakTime() ;
231 time -= presamplings;
233 plane = (*theHit)->WireID().Plane;
234 wire = (*theHit)->WireID().Wire;
237 if(plane == 1) time -= tIC;
242 wire_cm = (double)((wire+3.95) * wire_pitch);
244 wire_cm = (double)((wire+1.84) * wire_pitch);
247 if(time > tSI) time_cm = (double)( (time-tSI)*timepitch + tSI*driftvelocity_SI*timetick);
248 else time_cm = time*driftvelocity_SI*timetick;
250 wires.push_back(wire_cm);
251 times.push_back(time_cm);
253 the2Dtrack->SetPoint(np,wire_cm,time_cm);
259 the2Dtrack->Fit(
"pol1",
"Q");
266 TF1 *pol1=(TF1*) the2Dtrack->GetFunction(
"pol1");
268 pol1->GetParameters(par);
269 double intercept = par[0];
270 double slope = par[1];
272 double w0 = wires.front();
273 double w1 = wires.back();
274 double t0 = times.front();
275 double t1 = times.back();
276 double t0_line = intercept + (w0)*slope;
277 double t1_line = intercept + (w1)*slope;
282 Iwirefirsts.push_back(w0);
283 Iwirelasts.push_back(w1);
284 Itimefirsts.push_back(t0);
285 Itimelasts.push_back(t1);
286 Itimefirsts_line.push_back(t0_line);
287 Itimelasts_line.push_back(t1_line);
288 IclusHitlists.push_back(hitlist);
289 Icluster_count.push_back(ii);
292 Cwirefirsts.push_back(w0);
293 Cwirelasts.push_back(w1);
294 Ctimefirsts.push_back(t0);
295 Ctimelasts.push_back(t1);
296 Ctimefirsts_line.push_back(t0_line);
297 Ctimelasts_line.push_back(t1_line);
298 CclusHitlists.push_back(hitlist);
299 Ccluster_count.push_back(ii);
311 for(
size_t collectionIter = 0; collectionIter < CclusHitlists.size(); ++collectionIter){
313 double Cw0 = Cwirefirsts[collectionIter];
314 double Cw1 = Cwirelasts[collectionIter];
317 double Ct0_line = Ctimefirsts_line[collectionIter];
318 double Ct1_line = Ctimelasts_line[collectionIter];
319 std::vector< art::Ptr<recob::Hit> > hitsCtrk = CclusHitlists[collectionIter];
321 double collLength = TMath::Sqrt( TMath::Power(Ct1_line - Ct0_line,2) + TMath::Power(Cw1 - Cw0,2));
324 for(
size_t inductionIter = 0; inductionIter < IclusHitlists.size(); ++inductionIter){
326 double Iw0 = Iwirefirsts[inductionIter];
327 double Iw1 = Iwirelasts[inductionIter];
330 double It0_line = Itimefirsts_line[inductionIter];
331 double It1_line = Itimelasts_line[inductionIter];
332 std::vector< art::Ptr<recob::Hit> > hitsItrk = IclusHitlists[inductionIter];
334 double indLength = TMath::Sqrt( TMath::Power(It1_line - It0_line,2) + TMath::Power(Iw1 - Iw0,2));
336 bool forward_match = ((std::abs(Ct0_line-It0_line)<
ftmatch*timepitch) &&
337 (std::abs(Ct1_line-It1_line)<
ftmatch*timepitch));
338 bool backward_match = ((std::abs(Ct0_line-It1_line)<
ftmatch*timepitch) &&
339 (std::abs(Ct1_line-It0_line)<
ftmatch*timepitch));
344 if(forward_match || backward_match ){
349 XYZ0.SetXYZ(Ct0_line,(Cw0-Iw0)/(2.*TMath::Sin(Angle)),
350 (Cw0+Iw0)/(2.*TMath::Cos(Angle))-YC/2.*TMath::Tan(Angle));
351 XYZ1.SetXYZ(Ct1_line,(Cw1-Iw1)/(2.*TMath::Sin(Angle)),
352 (Cw1+Iw1)/(2.*TMath::Cos(Angle))-YC/2.*TMath::Tan(Angle));
355 XYZ0.SetXYZ(Ct0_line,(Cw0-Iw1)/(2.*TMath::Sin(Angle)),
356 (Cw0+Iw1)/(2.*TMath::Cos(Angle))-YC/2.*TMath::Tan(Angle));
357 XYZ1.SetXYZ(Ct1_line,(Cw1-Iw0)/(2.*TMath::Sin(Angle)),
358 (Cw1+Iw0)/(2.*TMath::Cos(Angle))-YC/2.*TMath::Tan(Angle));
364 TVector3 startpointVec,endpointVec;
365 TVector2 collVtx, indVtx;
366 if(XYZ0.Z() <= XYZ1.Z()){
367 startpointVec.SetXYZ(XYZ0.X(),XYZ0.Y(),XYZ0.Z());
368 endpointVec.SetXYZ(XYZ1.X(),XYZ1.Y(),XYZ1.Z());
370 collVtx.Set(Ct0_line,Cw0);
371 indVtx.Set(It0_line,Iw0);
374 collVtx.Set(Ct0_line,Cw0);
375 indVtx.Set(It1_line,Iw1);
379 startpointVec.SetXYZ(XYZ1.X(),XYZ1.Y(),XYZ1.Z());
380 endpointVec.SetXYZ(XYZ0.X(),XYZ0.Y(),XYZ0.Z());
382 collVtx.Set(Ct1_line,Cw1);
383 indVtx.Set(It1_line,Iw1);
386 collVtx.Set(Ct1_line,Cw1);
387 indVtx.Set(It0_line,Iw0);
392 TVector3 DirCos = endpointVec - startpointVec;
399 mf::LogWarning(
"Track3Dreco") <<
"The Spacepoint is infinitely small";
414 std::vector< art::Ptr<recob::Hit> > minhits = hitsCtrk.size() <= hitsItrk.size() ? hitsCtrk : hitsItrk;
415 std::vector< art::Ptr<recob::Hit> > maxhits = hitsItrk.size() <= hitsCtrk.size() ? hitsCtrk : hitsItrk;
418 std::vector<bool> maxhitsMatch(maxhits.size());
419 for(
size_t it = 0; it < maxhits.size(); ++it) maxhitsMatch[it] =
false;
421 std::vector<recob::Hit*> hits3Dmatched;
423 unsigned int imaximum = 0;
426 startSPIndex = spacepoints->size();
428 for(
size_t imin = 0; imin < minhits.size(); ++imin){
430 unsigned int wire,plane1,plane2;
431 wire = minhits[imin]->WireID().Wire;
432 plane1 = minhits[imin]->WireID().Plane;
437 w1 = (double)((wire+3.95)*wire_pitch);
439 w1 = (double)((wire+1.84) * wire_pitch);
440 double temptime1 = minhits[imin]->PeakTime()-presamplings;
441 if(plane1 == 1) temptime1 -= tIC;
443 if(temptime1>tSI) t1 = (double)( (temptime1-tSI)*timepitch + tSI*driftvelocity_SI*timetick);
444 else t1 = temptime1*driftvelocity_SI*timetick;
448 plane1==1 ? minVtx2D.Set(collVtx.X(),collVtx.Y()): minVtx2D.Set(indVtx.X(),indVtx.Y());
450 plane1==1 ? maxVtx2D.Set(indVtx.X(),indVtx.Y()): maxVtx2D.Set(collVtx.X(),collVtx.Y());
452 double ratio = (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.;
462 for(
size_t imax = 0; imax < maxhits.size(); ++imax){
463 if(!maxhitsMatch[imax]){
465 wire = maxhits[imax]->WireID().Wire;
466 plane2 = maxhits[imax]->WireID().Plane;
470 w2 = (double)((wire+3.95)*wire_pitch);
472 w2 = (double)((wire+1.84)*wire_pitch);
473 double temptime2 = maxhits[imax]->PeakTime()-presamplings;
474 if(plane2 == 1) temptime2 -= tIC;
476 if(temptime2 > tSI) t2 = (double)( (temptime2-tSI)*timepitch + tSI*driftvelocity_SI*timetick);
477 else t2 = temptime2*driftvelocity_SI*timetick;
480 bool timematch = (std::abs(t1-t2)<
ftmatch*timepitch);
481 bool wirematch = (std::abs(w1-w2)<wireShift*wire_pitch);
483 double maxDistance = TMath::Sqrt(TMath::Power(t2-maxVtx2D.X(),2)+TMath::Power(w2-maxVtx2D.Y(),2));
484 if (wirematch && timematch && std::abs(maxDistance-minDistance)<difference) {
485 difference = std::abs(maxDistance-minDistance);
490 maxhitsMatch[imaximum]=
true;
493 if(difference!= 9999999.){
499 wire = maxhits[imaximum]->WireID().Wire;
500 plane2 = maxhits[imaximum]->WireID().Plane;
504 w1_match = (double)((wire+3.95)*wire_pitch);
506 w1_match = (double)((wire+1.84)*wire_pitch);
507 double temptime3 = maxhits[imaximum]->PeakTime()-presamplings;
508 if(plane2 == 1) temptime3 -= tIC;
510 if(temptime3 > tSI) t1_match = (double)( (temptime3-tSI)*timepitch + tSI*driftvelocity_SI*timetick);
511 else t1_match = temptime3*driftvelocity_SI*timetick;
515 double Ct = plane1==1?t1:t1_match;
516 double Cw = plane1==1?w1:w1_match;
517 double Iw = plane1==1?w1_match:w1;
519 const TVector3 hit3d(Ct,(Cw-Iw)/(2.*TMath::Sin(Angle)),(Cw+Iw)/(2.*TMath::Cos(Angle))-YC/2.*TMath::Tan(Angle));
520 Double_t hitcoord[3];
521 hitcoord[0] = hit3d.X();
522 hitcoord[1] = hit3d.Y();
523 hitcoord[2] = hit3d.Z();
525 Double_t hitcoord_errs[3];
526 for (
int i=0; i<3; i++) hitcoord_errs[i]=-1.000;
531 spacepoints->push_back(mysp);
538 endSPIndex = spacepoints->size();
541 if(spacepoints->size() > startSPIndex || clustersPerTrack.
size()>0){
543 std::vector<TVector3> xyz;
544 xyz.push_back(startpointVec);
545 xyz.push_back(endpointVec);
546 std::vector<TVector3> dir_xyz;
547 dir_xyz.push_back(DirCos);
548 dir_xyz.push_back(DirCos);
554 tcol->push_back(the3DTrack);
557 util::CreateAssn(*
this, evt, *tcol, *spacepoints, *sassn, startSPIndex, endSPIndex);
565 for(
size_t p = 0; p < clustersPerTrack.
size(); ++p){
566 const std::vector< art::Ptr<recob::Hit> >&
hits = fmhc.at(p);
579 mf::LogVerbatim(
"Summary") << std::setfill(
'-') << std::setw(175) <<
"-" << std::setfill(
' ');
581 for(
unsigned int i = 0; i<tcol->size(); ++i){
585 evt.
put(std::move(tcol));
586 evt.
put(std::move(spacepoints));
587 evt.
put(std::move(cassn));
588 evt.
put(std::move(sassn));
589 evt.
put(std::move(hassn));
590 evt.
put(std::move(shassn));
code to link reconstructed objects back to the MC truth information
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
geo::Length_t PlanePitch(geo::TPCID const &tpcid, geo::PlaneID::PlaneID_t p1=0, geo::PlaneID::PlaneID_t p2=1) const
Returns the distance between two planes.
WireGeo const & Wire(unsigned int iwire) const
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 > > SMatrixSym55
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
std::string fClusterModuleLabel
label for input cluster collection
TrackTrajectory::Flags_t Flags_t
Planes which measure Z direction.
double ThetaZ() const
Returns angle of wire with respect to z axis in the Y-Z plane in radians.
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
ProductID put(std::unique_ptr< PROD > &&product)
geo::Length_t DetHalfHeight(geo::TPCID const &tpcid) const
Returns the half height of the active volume of the specified TPC.
void push_back(Ptr< U > const &p)
A trajectory in space reconstructed from hits.
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
int ftmatch
tolerance for time matching (in time samples)
data_t::const_iterator const_iterator
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
virtual double DriftVelocity(double efield=0., double temperature=0.) const =0
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
std::string GetLArTPCVolumeName(geo::TPCID const &tpcid) const
Return the name of specified LAr TPC volume.