138 std::unique_ptr<std::vector<recob::Track> > tcol (
new std::vector<recob::Track>);
139 std::unique_ptr<std::vector<recob::SpacePoint> > spcol(
new std::vector<recob::SpacePoint>);
151 double timetick = 0.198;
153 const double wireShift=50.;
156 double Efield_drift = 0.5;
157 double Efield_SI = 0.7;
158 double Efield_IC = 0.9;
159 double Temperature = 90.;
161 double driftvelocity = detprop->
DriftVelocity(Efield_drift,Temperature);
162 double driftvelocity_SI = detprop->
DriftVelocity(Efield_SI,Temperature);
163 double driftvelocity_IC = detprop->
DriftVelocity(Efield_IC,Temperature);
164 double timepitch = driftvelocity*timetick;
165 double tSI = plane_pitch/driftvelocity_SI/timetick;
166 double tIC = plane_pitch/driftvelocity_IC/timetick;
179 for (
unsigned int i = 0; i < endpointListHandle->size(); ++i){
186 std::vector<double> Iwirefirsts;
187 std::vector<double> Iwirelasts;
188 std::vector<double> Itimefirsts;
189 std::vector<double> Itimelasts;
190 std::vector<double> Itimefirsts_line;
191 std::vector<double> Itimelasts_line;
192 std::vector < std::vector< art::Ptr<recob::Hit> > > IclusHitlists;
193 std::vector<unsigned int> Icluster_count;
196 std::vector<double> Cwirefirsts;
197 std::vector<double> Cwirelasts;
198 std::vector<double> Ctimefirsts;
199 std::vector<double> Ctimelasts;
200 std::vector<double> Ctimefirsts_line;
201 std::vector<double> Ctimelasts_line;
202 std::vector< std::vector< art::Ptr<recob::Hit> > > CclusHitlists;
203 std::vector<unsigned int> Ccluster_count;
207 for(
unsigned int ii = 0; ii < clusterListHandle->size(); ++ii){
214 int vtx2d_w = -99999;
215 double vtx2d_t = -99999;
216 bool found2dvtx =
false;
218 for (
unsigned int j = 0; j<endpointlist.
size();j++){
219 if (endpointlist[j]->View() == cl->View()){
220 vtx2d_w = endpointlist[j]->WireID().Wire;
221 vtx2d_t = endpointlist[j]->DriftTime();
227 double w = cl->StartWire();
228 double t = cl->StartTick();
229 double dtdw = std::tan(cl->StartAngle());
230 double t_vtx = t+dtdw*(vtx2d_w-
w);
231 double dis = std::abs(vtx2d_t-t_vtx);
239 std::vector< art::Ptr<recob::Hit> > hitlist = fm.at(ii);
242 TGraph *the2Dtrack =
new TGraph(hitlist.size());
244 std::vector<double> wires;
245 std::vector<double> times;
253 time = (*theHit)->PeakTime() ;
255 time -= presamplings;
263 wire_cm = (double)(((*theHit)->WireID().Wire+3.95) * wire_pitch);
265 wire_cm = (double)(((*theHit)->WireID().Wire+1.84) * wire_pitch);
269 if(time>tSI) time_cm = (double)( (time-tSI)*timepitch + tSI*driftvelocity_SI*timetick);
270 else time_cm = time*driftvelocity_SI*timetick;
272 wires.push_back(wire_cm);
273 times.push_back(time_cm);
275 the2Dtrack->SetPoint(np,wire_cm,time_cm);
281 the2Dtrack->Fit(
"pol1",
"Q");
284 std::cout<<
"The 2D track fit failed"<<std::endl;
288 TF1 *pol1=(TF1*) the2Dtrack->GetFunction(
"pol1");
290 pol1->GetParameters(par);
291 double intercept = par[0];
292 double slope = par[1];
295 double w0 = wires.front();
296 double w1 = wires.back();
297 double t0 = times.front();
298 double t1 = times.back();
299 double t0_line = intercept + (w0)*slope;
300 double t1_line = intercept + (w1)*slope;
305 switch(geom->
SignalType((*hitlist.begin())->Channel())){
307 Iwirefirsts.push_back(w0);
308 Iwirelasts.push_back(w1);
309 Itimefirsts.push_back(t0);
310 Itimelasts.push_back(t1);
311 Itimefirsts_line.push_back(t0_line);
312 Itimelasts_line.push_back(t1_line);
313 IclusHitlists.push_back(hitlist);
314 Icluster_count.push_back(ii);
317 Cwirefirsts.push_back(w0);
318 Cwirelasts.push_back(w1);
319 Ctimefirsts.push_back(t0);
320 Ctimelasts.push_back(t1);
321 Ctimefirsts_line.push_back(t0_line);
322 Ctimelasts_line.push_back(t1_line);
323 CclusHitlists.push_back(hitlist);
324 Ccluster_count.push_back(ii);
336 for(
unsigned int collectionIter=0; collectionIter < CclusHitlists.size();collectionIter++){
338 double Cw0 = Cwirefirsts[collectionIter];
339 double Cw1 = Cwirelasts[collectionIter];
342 double Ct0_line = Ctimefirsts_line[collectionIter];
343 double Ct1_line = Ctimelasts_line[collectionIter];
344 std::vector< art::Ptr<recob::Hit> > hitsCtrk = CclusHitlists[collectionIter];
346 double collLength = TMath::Sqrt( TMath::Power(Ct1_line - Ct0_line,2) + TMath::Power(Cw1 - Cw0,2));
348 for(
unsigned int inductionIter=0;inductionIter<IclusHitlists.size();inductionIter++){
350 double Iw0 = Iwirefirsts[inductionIter];
351 double Iw1 = Iwirelasts[inductionIter];
354 double It0_line = Itimefirsts_line[inductionIter];
355 double It1_line = Itimelasts_line[inductionIter];
356 std::vector< art::Ptr<recob::Hit> > hitsItrk = IclusHitlists[inductionIter];
358 double indLength = TMath::Sqrt( TMath::Power(It1_line - It0_line,2) + TMath::Power(Iw1 - Iw0,2));
360 bool forward_match = ((std::abs(Ct0_line-It0_line)<
ftmatch*timepitch) && (std::abs(Ct1_line-It1_line)<
ftmatch*timepitch));
361 bool backward_match = ((std::abs(Ct0_line-It1_line)<
ftmatch*timepitch) && (std::abs(Ct1_line-It0_line)<
ftmatch*timepitch));
364 if(forward_match || backward_match ){
369 XYZ0.SetXYZ(Ct0_line,(Cw0-Iw0)/(2.*TMath::Sin(Angle)),(Cw0+Iw0)/(2.*TMath::Cos(Angle))-YC/2.*TMath::Tan(Angle));
370 XYZ1.SetXYZ(Ct1_line,(Cw1-Iw1)/(2.*TMath::Sin(Angle)),(Cw1+Iw1)/(2.*TMath::Cos(Angle))-YC/2.*TMath::Tan(Angle));
373 XYZ0.SetXYZ(Ct0_line,(Cw0-Iw1)/(2.*TMath::Sin(Angle)),(Cw0+Iw1)/(2.*TMath::Cos(Angle))-YC/2.*TMath::Tan(Angle));
374 XYZ1.SetXYZ(Ct1_line,(Cw1-Iw0)/(2.*TMath::Sin(Angle)),(Cw1+Iw0)/(2.*TMath::Cos(Angle))-YC/2.*TMath::Tan(Angle));
380 TVector3 startpointVec,endpointVec;
381 TVector2 collVtx, indVtx;
382 if(XYZ0.Z() <= XYZ1.Z()){
383 startpointVec.SetXYZ(XYZ0.X(),XYZ0.Y(),XYZ0.Z());
384 endpointVec.SetXYZ(XYZ1.X(),XYZ1.Y(),XYZ1.Z());
386 collVtx.Set(Ct0_line,Cw0);
387 indVtx.Set(It0_line,Iw0);
389 collVtx.Set(Ct0_line,Cw0);
390 indVtx.Set(It1_line,Iw1);
394 startpointVec.SetXYZ(XYZ1.X(),XYZ1.Y(),XYZ1.Z());
395 endpointVec.SetXYZ(XYZ0.X(),XYZ0.Y(),XYZ0.Z());
397 collVtx.Set(Ct1_line,Cw1);
398 indVtx.Set(It1_line,Iw1);
400 collVtx.Set(Ct1_line,Cw1);
401 indVtx.Set(It0_line,Iw0);
406 TVector3 DirCos = endpointVec - startpointVec;
413 catch(...){std::cout<<
"The Spacepoint is infinitely small"<<std::endl;
429 std::vector<recob::SpacePoint> spacepoints;
432 std::vector< art::Ptr<recob::Hit> > minhits = hitsCtrk.size() <= hitsItrk.size() ? hitsCtrk : hitsItrk;
433 std::vector< art::Ptr<recob::Hit> > maxhits = hitsItrk.size() < hitsCtrk.size() ? hitsCtrk : hitsItrk;
436 std::vector<bool> maxhitsMatch(maxhits.size());
437 for(
unsigned int it=0;it<maxhits.size();it++) maxhitsMatch[it] =
false;
439 std::vector<recob::Hit*> hits3Dmatched;
441 unsigned int imaximum = 0;
442 size_t spStart = spcol->size();
443 for(
unsigned int imin=0;imin<minhits.size();imin++){
447 auto const hitSigType = minhits[imin]->SignalType();
452 w1 = (double)((hit1WireID.
Wire+3.95) * wire_pitch);
454 w1 = (double)((hit1WireID.
Wire+1.84) * wire_pitch);
456 double temptime1 = minhits[imin]->PeakTime()-presamplings;
459 if(temptime1>tSI) t1 = (double)( (temptime1-tSI)*timepitch + tSI*driftvelocity_SI*timetick);
460 else t1 = temptime1*driftvelocity_SI*timetick;
464 (hitSigType ==
geo::kCollection) ? minVtx2D.Set(collVtx.X(),collVtx.Y()): minVtx2D.Set(indVtx.X(),indVtx.Y());
466 (hitSigType ==
geo::kCollection) ? maxVtx2D.Set(indVtx.X(),indVtx.Y()): maxVtx2D.Set(collVtx.X(),collVtx.Y());
468 double ratio = (collLength>indLength) ? collLength/indLength : indLength/collLength;
471 double minDistance = ratio*TMath::Sqrt(TMath::Power(t1-minVtx2D.X(),2) + TMath::Power(w1-minVtx2D.Y(),2));
475 double difference = 9999999.;
477 for(
unsigned int imax = 0; imax < maxhits.size(); imax++){
478 if(!maxhitsMatch[imax]){
481 auto const hit2SigType = maxhits[imax]->SignalType();
484 w2 = (double)((hit2WireID.
Wire+3.95) * wire_pitch);
486 w2 = (double)((hit2WireID.
Wire+1.84) * wire_pitch);
488 double temptime2 = maxhits[imax]->PeakTime()-presamplings;
491 if(temptime2>tSI) t2 = (double)( (temptime2-tSI)*timepitch + tSI*driftvelocity_SI*timetick);
492 else t2 = temptime2*driftvelocity_SI*timetick;
494 bool timematch = (std::abs(t1-t2)<
ftmatch*timepitch);
495 bool wirematch = (std::abs(w1-w2)<wireShift*wire_pitch);
497 double maxDistance = TMath::Sqrt(TMath::Power(t2-maxVtx2D.X(),2)+TMath::Power(w2-maxVtx2D.Y(),2));
498 if (wirematch && timematch && std::abs(maxDistance-minDistance)<difference) {
499 difference = std::abs(maxDistance-minDistance);
504 maxhitsMatch[imaximum]=
true;
507 if(difference!= 9999999.){
513 geo::WireID hit2WireID = maxhits[imaximum]->WireID();
514 auto const hit2SigType = maxhits[imaximum]->SignalType();
519 w1_match = (double)((hit2WireID.
Wire+3.95) * wire_pitch);
521 w1_match = (double)((hit2WireID.
Wire+1.84) * wire_pitch);
523 double temptime3 = maxhits[imaximum]->PeakTime()-presamplings;
526 if(temptime3>tSI) t1_match = (double)( (temptime3-tSI)*timepitch + tSI*driftvelocity_SI*timetick);
527 else t1_match = temptime3*driftvelocity_SI*timetick;
534 const TVector3 hit3d(Ct,(Cw-Iw)/(2.*TMath::Sin(Angle)),(Cw+Iw)/(2.*TMath::Cos(Angle))-YC/2.*TMath::Tan(Angle));
537 Double_t hitcoord[3];
538 hitcoord[0] = hit3d.X();
539 hitcoord[1] = hit3d.Y();
540 hitcoord[2] = hit3d.Z();
560 const double eps(0.1);
561 if (spacepoints.size()>=1){
562 TVector3 magNew(mysp.XYZ()[0],mysp.XYZ()[1],mysp.XYZ()[2]);
563 TVector3 magLast(spacepoints.back().XYZ()[0],
564 spacepoints.back().XYZ()[1],
565 spacepoints.back().XYZ()[2]);
566 if (!(magNew.Mag()>=magLast.Mag()+eps ||
567 magNew.Mag()<=magLast.Mag()-eps) )
570 spacepoints.push_back(mysp);
571 spcol->push_back(mysp);
576 size_t spEnd = spcol->size();
579 if(spacepoints.size()>0){
582 std::vector<TVector3> xyz(spacepoints.size());
583 for(
size_t s = 0;
s < spacepoints.size(); ++
s){
584 xyz[
s] = TVector3(spacepoints[
s].XYZ());
588 std::vector<TVector3> dircos(spacepoints.size(), DirCos);
603 for(
size_t cpt = 0; cpt < clustersPerTrack.
size(); ++cpt)
613 mf::LogVerbatim(
"Summary") << std::setfill(
'-') << std::setw(175) <<
"-" << std::setfill(
' ');
615 for(
unsigned int i = 0; i<tcol->size(); ++i)
mf::LogVerbatim(
"Summary") << tcol->at(i) ;
617 evt.
put(std::move(tcol));
618 evt.
put(std::move(spcol));
619 evt.
put(std::move(tspassn));
620 evt.
put(std::move(tcassn));
621 evt.
put(std::move(thassn));
622 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.
double fvertexclusterWindow
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)
TrackTrajectory::Flags_t Flags_t
WireID_t Wire
Index of the wire within its plane.
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
std::string fClusterModuleLabel
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)
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
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)
Signal from induction planes.
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)
bool SortByWire(art::Ptr< recob::Hit > const &h1, art::Ptr< recob::Hit > const &h2)
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
constexpr double kBogusD
obviously bogus double value
std::string fEndPoint2DModuleLabel
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.
Signal from collection planes.