Public Member Functions | Private Member Functions | Private Attributes | List of all members
trkf::SpacePts Class Reference
Inheritance diagram for trkf::SpacePts:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Member Functions

 SpacePts (fhicl::ParameterSet const &pset)
 
- Public Member Functions inherited from art::EDProducer
 EDProducer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDProducer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Producer
virtual ~Producer () noexcept
 
 Producer (fhicl::ParameterSet const &)
 
 Producer (Producer const &)=delete
 
 Producer (Producer &&)=delete
 
Produceroperator= (Producer const &)=delete
 
Produceroperator= (Producer &&)=delete
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Modifier
 ~Modifier () noexcept
 
 Modifier ()
 
 Modifier (Modifier const &)=delete
 
 Modifier (Modifier &&)=delete
 
Modifieroperator= (Modifier const &)=delete
 
Modifieroperator= (Modifier &&)=delete
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Private Member Functions

void produce (art::Event &evt)
 

Private Attributes

int ftmatch
 
double fPreSamplings
 
double fvertexclusterWindow
 
std::string fClusterModuleLabel
 
std::string fEndPoint2DModuleLabel
 

Additional Inherited Members

- Public Types inherited from art::EDProducer
using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
- Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 
- Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
 
- Static Public Member Functions inherited from art::EDProducer
static void commitEvent (EventPrincipal &ep, Event &e)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollectorconsumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

Definition at line 50 of file SpacePts_module.cc.

Constructor & Destructor Documentation

trkf::SpacePts::SpacePts ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 73 of file SpacePts_module.cc.

73  : EDProducer{pset}
74  {
75  fPreSamplings = pset.get<double>("TicksOffset");
76  ftmatch = pset.get<int>("TMatch");
77  fClusterModuleLabel = pset.get<std::string>("ClusterModuleLabel");
78  fEndPoint2DModuleLabel = pset.get<std::string>("EndPoint2DModuleLabel");
79  fvertexclusterWindow = pset.get<double>("vertexclusterWindow");
80 
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>>();
87  }
double fvertexclusterWindow
std::string string
Definition: nybbler.cc:12
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
std::string fClusterModuleLabel
std::string fEndPoint2DModuleLabel

Member Function Documentation

void trkf::SpacePts::produce ( art::Event evt)
privatevirtual
Todo:
really should fill the direction cosines with unique values

Implements art::EDProducer.

Definition at line 91 of file SpacePts_module.cc.

92  {
94  auto const detProp =
96 
97  //////////////////////////////////////////////////////
98  // Make a std::unique_ptr<> for the thing you want to put into the event
99  // because that handles the memory management for you
100  //////////////////////////////////////////////////////
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>>();
107 
108  // define TPC parameters
109  TString tpcName = geom->GetLArTPCVolumeName();
110 
111  //TPC dimensions
112  double YC = (geom->DetHalfHeight()) * 2.; // TPC height in cm
113  double Angle = geom->Plane(1).Wire(0).ThetaZ(false) -
114  TMath::Pi() / 2.; // wire angle with respect to the vertical direction
115  // Parameters temporary defined here, but possibly to be retrieved somewhere in the code
116  double timetick = 0.198; //time sample in us
117  double presamplings = fPreSamplings; // 60.;
118  const double wireShift =
119  50.; // half the number of wires from the Induction(Collection) plane intersecting with a wire from the Collection(Induction) plane.
120  double plane_pitch = geom->PlanePitch(0, 1); //wire plane pitch in cm
121  double wire_pitch = geom->WirePitch(); //wire pitch in cm
122  double Efield_drift = 0.5; // Electric Field in the drift region in kV/cm
123  double Efield_SI = 0.7; // Electric Field between Shield and Induction planes in kV/cm
124  double Efield_IC = 0.9; // Electric Field between Induction and Collection planes in kV/cm
125  double Temperature = 90.; // LAr Temperature in K
126 
127  double driftvelocity =
128  detProp.DriftVelocity(Efield_drift, Temperature); //drift velocity in the drift region (cm/us)
129  double driftvelocity_SI = detProp.DriftVelocity(
130  Efield_SI, Temperature); //drift velocity between shield and induction (cm/us)
131  double driftvelocity_IC = detProp.DriftVelocity(
132  Efield_IC, Temperature); //drift velocity between induction and collection (cm/us)
133  double timepitch = driftvelocity * timetick; //time sample (cm)
134  double tSI = plane_pitch / driftvelocity_SI /
135  timetick; //drift time between Shield and Collection planes (time samples)
136  double tIC = plane_pitch / driftvelocity_IC /
137  timetick; //drift time between Induction and Collection planes (time samples)
138 
139  // get input Cluster object(s).
140  art::Handle<std::vector<recob::Cluster>> clusterListHandle;
141  evt.getByLabel(fClusterModuleLabel, clusterListHandle);
142 
143  // get input EndPoint2D object(s).
144  art::Handle<std::vector<recob::EndPoint2D>> endpointListHandle;
145  evt.getByLabel(fEndPoint2DModuleLabel, endpointListHandle);
146 
148  if (evt.getByLabel(fEndPoint2DModuleLabel, endpointListHandle))
149  for (unsigned int i = 0; i < endpointListHandle->size(); ++i) {
150  art::Ptr<recob::EndPoint2D> endpointHolder(endpointListHandle, i);
151  endpointlist.push_back(endpointHolder);
152  }
153 
154  // Declare some vectors..
155  // Induction
156  std::vector<double> Iwirefirsts; // in cm
157  std::vector<double> Iwirelasts; // in cm
158  std::vector<double> Itimefirsts; // in cm
159  std::vector<double> Itimelasts; // in cm
160  std::vector<double> Itimefirsts_line; // in cm
161  std::vector<double> Itimelasts_line; // in cm
162  std::vector<std::vector<art::Ptr<recob::Hit>>> IclusHitlists;
163  std::vector<unsigned int> Icluster_count;
164 
165  // Collection
166  std::vector<double> Cwirefirsts; // in cm
167  std::vector<double> Cwirelasts; // in cm
168  std::vector<double> Ctimefirsts; // in cm
169  std::vector<double> Ctimelasts; // in cm
170  std::vector<double> Ctimefirsts_line; // in cm
171  std::vector<double> Ctimelasts_line; // in cm
172  std::vector<std::vector<art::Ptr<recob::Hit>>> CclusHitlists;
173  std::vector<unsigned int> Ccluster_count;
174 
175  art::FindManyP<recob::Hit> fm(clusterListHandle, evt, fClusterModuleLabel);
176 
177  for (unsigned int ii = 0; ii < clusterListHandle->size(); ++ii) {
178  art::Ptr<recob::Cluster> cl(clusterListHandle, ii);
179 
180  // Figure out which View the cluster belongs to
181  //only consider merged-lines that are associated with the vertex.
182  //this helps get rid of through-going muon background -spitz
183  int vtx2d_w = -99999;
184  double vtx2d_t = -99999;
185  bool found2dvtx = false;
186 
187  for (unsigned int j = 0; j < endpointlist.size(); j++) {
188  if (endpointlist[j]->View() == cl->View()) {
189  vtx2d_w = endpointlist[j]->WireID().Wire; //for update to EndPoint2D ... WK 4/22/13
190  vtx2d_t = endpointlist[j]->DriftTime();
191  found2dvtx = true;
192  break;
193  }
194  }
195  if (found2dvtx) {
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);
201  if (dis > fvertexclusterWindow) continue;
202  }
203  //else continue; //what to do if a 2D vertex is not found? perhaps vertex finder was not even run.
204 
205  // Some variables for the hit
206  float time; //hit time at maximum
207 
208  std::vector<art::Ptr<recob::Hit>> hitlist = fm.at(ii);
209  std::sort(hitlist.begin(), hitlist.end(), trkf::SortByWire());
210 
211  TGraph* the2Dtrack = new TGraph(hitlist.size());
212 
213  std::vector<double> wires;
214  std::vector<double> times;
215 
216  int np = 0;
217  //loop over cluster hits
218  for (std::vector<art::Ptr<recob::Hit>>::const_iterator theHit = hitlist.begin();
219  theHit != hitlist.end();
220  theHit++) {
221  //recover the Hit
222  // recob::Hit* theHit = (recob::Hit*)(*hitIter);
223  time = (*theHit)->PeakTime();
224 
225  time -= presamplings;
226 
227  if (geom->SignalType((*theHit)->Channel()) == geo::kCollection) time -= tIC; // Collection
228  //transform hit wire and time into cm
229  double wire_cm = 0.;
230  if (geom->SignalType((*theHit)->Channel()) == geo::kInduction)
231  wire_cm = (double)(((*theHit)->WireID().Wire + 3.95) * wire_pitch);
232  else
233  wire_cm = (double)(((*theHit)->WireID().Wire + 1.84) * wire_pitch);
234 
235  //double time_cm = (double)(time * timepitch);
236  double time_cm;
237  if (time > tSI)
238  time_cm = (double)((time - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
239  else
240  time_cm = time * driftvelocity_SI * timetick;
241 
242  wires.push_back(wire_cm);
243  times.push_back(time_cm);
244 
245  the2Dtrack->SetPoint(np, wire_cm, time_cm);
246  np++;
247  } //end of loop over cluster hits
248 
249  // fit the 2Dtrack and get some info to store
250  try {
251  the2Dtrack->Fit("pol1", "Q");
252  }
253  catch (...) {
254  std::cout << "The 2D track fit failed" << std::endl;
255  continue;
256  }
257 
258  TF1* pol1 = (TF1*)the2Dtrack->GetFunction("pol1");
259  double par[2];
260  pol1->GetParameters(par);
261  double intercept = par[0];
262  double slope = par[1];
263 
264  double w0 = wires.front(); // first hit wire (cm)
265  double w1 = wires.back(); // last hit wire (cm)
266  double t0 = times.front(); // first hit time (cm)
267  double t1 = times.back(); // last hit time (cm)
268  double t0_line = intercept + (w0)*slope; // time coordinate at wire w0 on the fit line (cm)
269  double t1_line = intercept + (w1)*slope; // time coordinate at wire w1 on the fit line (cm)
270 
271  // actually store the 2Dtrack info
272  switch (geom->SignalType((*hitlist.begin())->Channel())) {
273  case geo::kInduction:
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);
282  break;
283  case geo::kCollection:
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);
292  break;
293  case geo::kMysteryType: break;
294  }
295  delete pol1;
296  } // end of loop over all input clusters
297 
298  /////////////////////////////////////////////////////
299  /////// 2D Track Matching and 3D Track Reconstruction
300  /////////////////////////////////////////////////////
301 
302  for (unsigned int collectionIter = 0; collectionIter < CclusHitlists.size();
303  collectionIter++) { //loop over Collection view 2D tracks
304  // Recover previously stored info
305  double Cw0 = Cwirefirsts[collectionIter];
306  double Cw1 = Cwirelasts[collectionIter];
307  //double Ct0 = Ctimefirsts[collectionIter];
308  //double Ct1 = Ctimelasts[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];
312 
313  double collLength =
314  TMath::Sqrt(TMath::Power(Ct1_line - Ct0_line, 2) + TMath::Power(Cw1 - Cw0, 2));
315 
316  for (unsigned int inductionIter = 0; inductionIter < IclusHitlists.size();
317  inductionIter++) { //loop over Induction view 2D tracks
318  // Recover previously stored info
319  double Iw0 = Iwirefirsts[inductionIter];
320  double Iw1 = Iwirelasts[inductionIter];
321  //double It0 = Itimefirsts[inductionIter];
322  //double It1 = Itimelasts[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];
326 
327  double indLength =
328  TMath::Sqrt(TMath::Power(It1_line - It0_line, 2) + TMath::Power(Iw1 - Iw0, 2));
329 
330  bool forward_match = ((std::abs(Ct0_line - It0_line) < ftmatch * timepitch) &&
331  (std::abs(Ct1_line - It1_line) < ftmatch * timepitch));
332  bool backward_match = ((std::abs(Ct0_line - It1_line) < ftmatch * timepitch) &&
333  (std::abs(Ct1_line - It0_line) < ftmatch * timepitch));
334 
335  if (forward_match || backward_match) {
336 
337  // Reconstruct the 3D track
338  TVector3 XYZ0, XYZ1; // track endpoints
339  if (forward_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));
346  }
347  else {
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));
354  }
355 
356  //compute track direction in Local co-ordinate system
357  //WARNING: There is an ambiguity introduced here for the case of backwards-going tracks.
358  //If available, vertex info. could sort this out.
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());
364  if (forward_match) {
365  collVtx.Set(Ct0_line, Cw0);
366  indVtx.Set(It0_line, Iw0);
367  }
368  else {
369  collVtx.Set(Ct0_line, Cw0);
370  indVtx.Set(It1_line, Iw1);
371  }
372  }
373  else {
374  startpointVec.SetXYZ(XYZ1.X(), XYZ1.Y(), XYZ1.Z());
375  endpointVec.SetXYZ(XYZ0.X(), XYZ0.Y(), XYZ0.Z());
376  if (forward_match) {
377  collVtx.Set(Ct1_line, Cw1);
378  indVtx.Set(It1_line, Iw1);
379  }
380  else {
381  collVtx.Set(Ct1_line, Cw1);
382  indVtx.Set(It0_line, Iw0);
383  }
384  }
385 
386  //compute track (normalized) cosine directions in the TPC co-ordinate system
387  TVector3 DirCos = endpointVec - startpointVec;
388 
389  //SetMag casues a crash if the magnitude of the vector is zero
390  try {
391  DirCos.SetMag(1.0); //normalize vector
392  }
393  catch (...) {
394  std::cout << "The Spacepoint is infinitely small" << std::endl;
395  continue;
396  }
397 
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);
403 
404  /////////////////////////////
405  // Match hits
406  ////////////////////////////
407 
408  //create collection of spacepoints that will be used when creating the Track object
409  std::vector<recob::SpacePoint> spacepoints;
410 
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;
415 
416  std::vector<bool> maxhitsMatch(maxhits.size());
417  for (unsigned int it = 0; it < maxhits.size(); it++)
418  maxhitsMatch[it] = false;
419 
420  std::vector<recob::Hit*> hits3Dmatched;
421  // For the matching start from the view where the track projection presents less hits
422  unsigned int imaximum = 0;
423  size_t spStart = spcol->size();
424  for (unsigned int imin = 0; imin < minhits.size(); imin++) { //loop over hits
425  //get wire - time coordinate of the hit
426  //unsigned int channel,wire,plane1,plane2,tpc,cstat;
427  geo::WireID hit1WireID = minhits[imin]->WireID();
428  auto const hitSigType = minhits[imin]->SignalType();
429  double w1 = 0;
430 
431  //the 3.95 and 1.84 below are the ArgoNeuT TPC offsets for the induction and collection plane, respectively and are in units of wire pitch.
432  if (hitSigType == geo::kInduction)
433  w1 = (double)((hit1WireID.Wire + 3.95) * wire_pitch);
434  else
435  w1 = (double)((hit1WireID.Wire + 1.84) * wire_pitch);
436 
437  double temptime1 = minhits[imin]->PeakTime() - presamplings;
438  if (hitSigType == geo::kCollection) temptime1 -= tIC;
439  double
440  t1; // = plane1==1?(double)((minhits[imin]->PeakTime()-presamplings-tIC)*timepitch):(double)((minhits[imin]->PeakTime()-presamplings)*timepitch); //in cm
441  if (temptime1 > tSI)
442  t1 = (double)((temptime1 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
443  else
444  t1 = temptime1 * driftvelocity_SI * timetick;
445 
446  //get the track origin co-ordinates in the two views
447  TVector2 minVtx2D;
448  (hitSigType == geo::kCollection) ? minVtx2D.Set(collVtx.X(), collVtx.Y()) :
449  minVtx2D.Set(indVtx.X(), indVtx.Y());
450  TVector2 maxVtx2D;
451  (hitSigType == geo::kCollection) ? maxVtx2D.Set(indVtx.X(), indVtx.Y()) :
452  maxVtx2D.Set(collVtx.X(), collVtx.Y());
453 
454  double ratio =
455  (collLength > indLength) ? collLength / indLength : indLength / collLength;
456 
457  //compute the distance of the hit (imin) from the relative track origin
458  double minDistance = ratio * TMath::Sqrt(TMath::Power(t1 - minVtx2D.X(), 2) +
459  TMath::Power(w1 - minVtx2D.Y(), 2));
460 
461  //core matching algorithm
462  double difference = 9999999.;
463 
464  for (unsigned int imax = 0; imax < maxhits.size();
465  imax++) { //loop over hits of the other view
466  if (!maxhitsMatch[imax]) {
467  //get wire - time coordinate of the hit
468  geo::WireID hit2WireID = maxhits[imax]->WireID();
469  auto const hit2SigType = maxhits[imax]->SignalType();
470  double w2 = 0.;
471  if (hit2SigType == geo::kInduction)
472  w2 = (double)((hit2WireID.Wire + 3.95) * wire_pitch);
473  else
474  w2 = (double)((hit2WireID.Wire + 1.84) * wire_pitch);
475 
476  double temptime2 = maxhits[imax]->PeakTime() - presamplings;
477  if (hit2SigType == geo::kCollection) temptime2 -= tIC;
478  double t2;
479  if (temptime2 > tSI)
480  t2 = (double)((temptime2 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
481  else
482  t2 = temptime2 * driftvelocity_SI * timetick;
483 
484  bool timematch = (std::abs(t1 - t2) < ftmatch * timepitch);
485  bool wirematch = (std::abs(w1 - w2) < wireShift * wire_pitch);
486 
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);
491  imaximum = imax;
492  }
493  }
494  }
495  maxhitsMatch[imaximum] = true;
496 
498  if (difference != 9999999.) {
499  sp_hits.push_back(minhits[imin]);
500  sp_hits.push_back(maxhits[imaximum]);
501  }
502 
503  // Get the time-wire co-ordinates of the matched hit
504  geo::WireID hit2WireID = maxhits[imaximum]->WireID();
505  auto const hit2SigType = maxhits[imaximum]->SignalType();
506 
507  //double w1_match = (double)((wire+1)*wire_pitch);
508  double w1_match = 0.;
509  if (hit2SigType == geo::kInduction)
510  w1_match = (double)((hit2WireID.Wire + 3.95) * wire_pitch);
511  else
512  w1_match = (double)((hit2WireID.Wire + 1.84) * wire_pitch);
513 
514  double temptime3 = maxhits[imaximum]->PeakTime() - presamplings;
515  if (hit2SigType == geo::kCollection) temptime3 -= tIC;
516  double t1_match;
517  if (temptime3 > tSI)
518  t1_match =
519  (double)((temptime3 - tSI) * timepitch + tSI * driftvelocity_SI * timetick);
520  else
521  t1_match = temptime3 * driftvelocity_SI * timetick;
522 
523  // create the 3D hit, compute its co-ordinates and add it to the 3D hits list
524  double Ct = hitSigType == geo::kCollection ? t1 : t1_match;
525  double Cw = hit2SigType == geo::kCollection ? w1 : w1_match;
526  double Iw = hit2SigType == geo::kCollection ? w1_match : w1;
527 
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));
532 
533  Double_t hitcoord[3];
534  hitcoord[0] = hit3d.X();
535  hitcoord[1] = hit3d.Y();
536  hitcoord[2] = hit3d.Z();
537 
538  /*
539  double yy,zz;
540  if(geom->ChannelsIntersect(geom->PlaneWireToChannel(0,(int)((Iw/wire_pitch)-3.95)), geom->PlaneWireToChannel(1,(int)((Cw/wire_pitch)-1.84)),yy,zz))
541  {
542  //channelsintersect provides a slightly more accurate set of y and z coordinates. use channelsintersect in case the wires in question do cross.
543  hitcoord[1] = yy;
544  hitcoord[2] = zz;
545  mf::LogInfo("SpacePts: ") << "SpacePoint adding xyz ..." << hitcoord[0] <<","<< hitcoord[1] <<","<< hitcoord[2];
546  // std::cout<<"wire 1: "<<(Iw/wire_pitch)-3.95<<" "<<(Cw/wire_pitch)-1.84<<std::endl;
547  // std::cout<<"Intersect: "<<yy<<" "<<zz<<std::endl;
548  }
549  else
550  continue;
551  */
552 
553  double err[6] = {util::kBogusD};
554  recob::SpacePoint mysp(hitcoord,
555  err,
557  spStart + spacepoints.size()); //3d point at end of track
558  // Don't add a spacepoint right on top of the last one.
559  const double eps(0.1); // 1mm
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))
566  continue;
567  }
568  spacepoints.push_back(mysp);
569  spcol->push_back(mysp);
570  util::CreateAssn(*this, evt, *spcol, sp_hits, *shassn);
571 
572  } //loop over min-hits
573 
574  size_t spEnd = spcol->size();
575 
576  // Add the 3D track to the vector of the reconstructed tracks
577  if (spacepoints.size() > 0) {
578 
579  // make a vector of the trajectory points along the track
580  std::vector<TVector3> xyz(spacepoints.size());
581  for (size_t s = 0; s < spacepoints.size(); ++s) {
582  xyz[s] = TVector3(spacepoints[s].XYZ());
583  }
584 
585  ///\todo really should fill the direction cosines with unique values
586  std::vector<TVector3> dircos(spacepoints.size(), DirCos);
587 
588  tcol->push_back(
591  recob::Track::Flags_t(xyz.size()),
592  false),
593  0,
594  -1.,
595  0,
598  tcol->size()));
599 
600  // make associations between the track and space points
601  util::CreateAssn(*this, evt, *tcol, *spcol, *tspassn, spStart, spEnd);
602 
603  // now the track and clusters
604  util::CreateAssn(*this, evt, *tcol, clustersPerTrack, *tcassn);
605 
606  // and the hits and track
607  art::FindManyP<recob::Hit> fmh(clustersPerTrack, evt, fClusterModuleLabel);
608  for (size_t cpt = 0; cpt < clustersPerTrack.size(); ++cpt)
609  util::CreateAssn(*this, evt, *tcol, fmh.at(cpt), *thassn);
610  }
611  } //close match 2D tracks
612 
613  } //close loop over Induction view 2D tracks
614 
615  } //close loop over Collection xxview 2D tracks
616 
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);
621 
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));
628 
629  } // end SpacePts::produce()
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.
Who knows?
Definition: geo_types.h:147
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.
AdcChannelData::View View
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 > > SMatrixSym55
Definition: TrackingTypes.h:85
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Definition: TrackingTypes.h:58
struct vector vector
TrackTrajectory::Flags_t Flags_t
Definition: Track.h:68
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
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.
Definition: WireGeo.h:250
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
QAsciiDict< Entry > cl
T abs(T value)
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:633
bool SortByWire(art::Ptr< recob::Hit > const &h1, art::Ptr< recob::Hit > const &h2)
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)
Definition: PtrVector.h:435
Signal from induction planes.
Definition: geo_types.h:145
A trajectory in space reconstructed from hits.
def move(depos, offset)
Definition: depos.py:107
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
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)
Definition: TrackingTypes.h:55
void err(const char *fmt,...)
Definition: message.cpp:226
size_type size() const
Definition: PtrVector.h:302
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
int imax
Definition: tracks.py:195
static constexpr double fm
Definition: Units.h:75
constexpr double kBogusD
obviously bogus double value
std::string fEndPoint2DModuleLabel
Q_EXPORT QTSManip setfill(int f)
Definition: qtextstream.h:337
Definition: fwd.h:31
static QCString * s
Definition: config.cpp:1042
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
QTextStream & endl(QTextStream &s)
std::string GetLArTPCVolumeName(geo::TPCID const &tpcid) const
Return the name of specified LAr TPC volume.
Signal from collection planes.
Definition: geo_types.h:146

Member Data Documentation

std::string trkf::SpacePts::fClusterModuleLabel
private

Definition at line 60 of file SpacePts_module.cc.

std::string trkf::SpacePts::fEndPoint2DModuleLabel
private

Definition at line 61 of file SpacePts_module.cc.

double trkf::SpacePts::fPreSamplings
private

Definition at line 58 of file SpacePts_module.cc.

int trkf::SpacePts::ftmatch
private

Definition at line 57 of file SpacePts_module.cc.

double trkf::SpacePts::fvertexclusterWindow
private

Definition at line 59 of file SpacePts_module.cc.


The documentation for this class was generated from the following file: