LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
trkf::SpacePts Class Reference
Inheritance diagram for trkf::SpacePts:
art::EDProducer art::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = ProducerBase::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 SpacePts (fhicl::ParameterSet const &pset)
 
 ~SpacePts ()
 
void reconfigure (fhicl::ParameterSet const &p)
 
void produce (art::Event &evt)
 
void beginJob ()
 
void endJob ()
 
template<typename PROD , BranchType B = InEvent>
ProductID getProductID (std::string const &instanceName={}) const
 
template<typename PROD , BranchType B>
ProductID getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
 
bool modifiesEvent () const
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< T > consumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< T > consumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< T > mayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< T > mayConsumeView (InputTag const &it)
 
base_engine_tcreateEngine (seed_t seed)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumernon_module_context ()
 

Protected Member Functions

CurrentProcessingContext const * currentContext () const
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

Private Attributes

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

Detailed Description

Definition at line 53 of file SpacePts_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 43 of file EDProducer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>
inherited

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

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

Definition at line 90 of file SpacePts_module.cc.

References reconfigure().

91 {
92  this->reconfigure(pset);
93 
94  produces< std::vector<recob::Track> >();
95  produces< std::vector<recob::SpacePoint> >();
96  produces< art::Assns<recob::Track, recob::SpacePoint> >();
97  produces< art::Assns<recob::Track, recob::Cluster> >();
98  produces< art::Assns<recob::Track, recob::Hit> >();
99  produces< art::Assns<recob::SpacePoint, recob::Hit> >();
100 }
void reconfigure(fhicl::ParameterSet const &p)
trkf::SpacePts::~SpacePts ( )

Definition at line 103 of file SpacePts_module.cc.

104 {
105 }

Member Function Documentation

void trkf::SpacePts::beginJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 117 of file SpacePts_module.cc.

118 {
119 }
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::consumes ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::consumes ( InputTag const &  it)
inherited

Definition at line 147 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

148 {
149  if (!moduleContext_)
150  return ProductToken<T>::invalid();
151 
152  consumables_[BT].emplace_back(ConsumableType::Product,
153  TypeID{typeid(T)},
154  it.label(),
155  it.instance(),
156  it.process());
157  return ProductToken<T>{it};
158 }
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 162 of file Consumer.h.

163 {
164  if (!moduleContext_)
165  return;
166 
167  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
168 }
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 172 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

173 {
174  if (!moduleContext_)
175  return ViewToken<T>::invalid();
176 
177  consumables_[BT].emplace_back(ConsumableType::ViewElement,
178  TypeID{typeid(T)},
179  it.label(),
180  it.instance(),
181  it.process());
182  return ViewToken<T>{it};
183 }
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited

Definition at line 32 of file EngineCreator.cc.

References art::EngineCreator::rng().

34 {
35  return rng()->createEngine(
36  placeholder_schedule_id(), seed, kind_of_engine_to_make);
37 }
long seed
Definition: chem4.cc:68
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited

Definition at line 40 of file EngineCreator.cc.

References art::EngineCreator::rng().

43 {
44  return rng()->createEngine(
45  placeholder_schedule_id(), seed, kind_of_engine_to_make, engine_label);
46 }
long seed
Definition: chem4.cc:68
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
CurrentProcessingContext const * art::EDProducer::currentContext ( ) const
protectedinherited

Definition at line 120 of file EDProducer.cc.

References art::EDProducer::current_context_.

121  {
122  return current_context_.get();
123  }
CPC_exempt_ptr current_context_
Definition: EDProducer.h:116
void trkf::SpacePts::endJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 121 of file SpacePts_module.cc.

122 {
123 }
EngineCreator::seed_t EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited

Definition at line 49 of file EngineCreator.cc.

References fhicl::ParameterSet::get().

Referenced by art::MixFilter< T >::initEngine_().

52 {
53  auto const& explicit_seeds = pset.get<std::vector<int>>(key, {});
54  return explicit_seeds.empty() ? implicit_seed : explicit_seeds.front();
55 }
template<typename PROD , BranchType B>
ProductID art::EDProducer::getProductID ( std::string const &  instanceName = {}) const
inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

124  {
125  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
126  instanceName);
127  }
ModuleDescription moduleDescription_
Definition: EDProducer.h:115
template<typename PROD , BranchType B>
ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription,
std::string const &  instanceName 
) const
inherited

Definition at line 56 of file ProducerBase.h.

References B, and art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
Int_t B
Definition: plot.C:25
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::mayConsume ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::mayConsume ( InputTag const &  it)
inherited

Definition at line 190 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

191 {
192  if (!moduleContext_)
193  return ProductToken<T>::invalid();
194 
195  consumables_[BT].emplace_back(ConsumableType::Product,
196  TypeID{typeid(T)},
197  it.label(),
198  it.instance(),
199  it.process());
200  return ProductToken<T>{it};
201 }
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 205 of file Consumer.h.

206 {
207  if (!moduleContext_)
208  return;
209 
210  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
211 }
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it)
inherited

Definition at line 215 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

216 {
217  if (!moduleContext_)
218  return ViewToken<T>::invalid();
219 
220  consumables_[BT].emplace_back(ConsumableType::ViewElement,
221  TypeID{typeid(T)},
222  it.label(),
223  it.instance(),
224  it.process());
225  return ViewToken<T>{it};
226 }
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
bool art::ProducerBase::modifiesEvent ( ) const
inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID().

41  {
42  return true;
43  }
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited

Definition at line 89 of file Consumer.cc.

References fhicl::ParameterSet::get_if_present().

Referenced by art::EDProducer::doBeginJob(), art::EDFilter::doBeginJob(), and art::EDAnalyzer::doBeginJob().

90 {
91  if (!moduleContext_)
92  return;
93 
94  pset.get_if_present("errorOnMissingConsumes", requireConsumes_);
95  for (auto& consumablesPerBranch : consumables_) {
96  cet::sort_all(consumablesPerBranch);
97  }
98 }
bool requireConsumes_
Definition: Consumer.h:137
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
void trkf::SpacePts::produce ( art::Event evt)
virtual
Todo:
really should fill the direction cosines with unique values

Implements art::EDProducer.

Definition at line 126 of file SpacePts_module.cc.

References recob::tracking::convertCollToPoint(), recob::tracking::convertCollToVector(), util::CreateAssn(), DEFINE_ART_MODULE, geo::GeometryCore::DetHalfHeight(), detinfo::DetectorProperties::DriftVelocity(), fClusterModuleLabel, fEndPoint2DModuleLabel, fPreSamplings, ftmatch, fvertexclusterWindow, art::DataViewImpl::getByLabel(), geo::GeometryCore::GetLArTPCVolumeName(), util::kBogusD, geo::kCollection, geo::kInduction, geo::kMysteryType, geo::GeometryCore::Plane(), geo::GeometryCore::PlanePitch(), art::PtrVector< T >::push_back(), art::Event::put(), s, geo::GeometryCore::SignalType(), art::PtrVector< T >::size(), SortByWire(), recob::Cluster::StartAngle(), recob::Cluster::StartTick(), recob::Cluster::StartWire(), t1, t2, geo::WireGeo::ThetaZ(), lar::dump::vector(), recob::Cluster::View(), w, geo::PlaneGeo::Wire(), geo::WireID::Wire, and geo::GeometryCore::WirePitch().

127 {
128 
129 
130  // get services
132  const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
133 
135  // Make a std::unique_ptr<> for the thing you want to put into the event
136  // because that handles the memory management for you
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>);
140  std::unique_ptr<art::Assns<recob::Track, recob::SpacePoint> > tspassn(new art::Assns<recob::Track, recob::SpacePoint>);
141  std::unique_ptr<art::Assns<recob::Track, recob::Cluster> > tcassn(new art::Assns<recob::Track, recob::Cluster>);
142  std::unique_ptr<art::Assns<recob::Track, recob::Hit> > thassn(new art::Assns<recob::Track, recob::Hit>);
143  std::unique_ptr<art::Assns<recob::SpacePoint, recob::Hit> > shassn(new art::Assns<recob::SpacePoint, recob::Hit>);
144  // define TPC parameters
145  TString tpcName = geom->GetLArTPCVolumeName();
146 
147  //TPC dimensions
148  double YC = (geom->DetHalfHeight())*2.; // TPC height in cm
149  double Angle = geom->Plane(1).Wire(0).ThetaZ(false)-TMath::Pi()/2.; // wire angle with respect to the vertical direction
150  // Parameters temporary defined here, but possibly to be retrieved somewhere in the code
151  double timetick = 0.198; //time sample in us
152  double presamplings = fPreSamplings; // 60.;
153  const double wireShift=50.; // half the number of wires from the Induction(Collection) plane intersecting with a wire from the Collection(Induction) plane.
154  double plane_pitch = geom->PlanePitch(0,1); //wire plane pitch in cm
155  double wire_pitch = geom->WirePitch(); //wire pitch in cm
156  double Efield_drift = 0.5; // Electric Field in the drift region in kV/cm
157  double Efield_SI = 0.7; // Electric Field between Shield and Induction planes in kV/cm
158  double Efield_IC = 0.9; // Electric Field between Induction and Collection planes in kV/cm
159  double Temperature = 90.; // LAr Temperature in K
160 
161  double driftvelocity = detprop->DriftVelocity(Efield_drift,Temperature); //drift velocity in the drift region (cm/us)
162  double driftvelocity_SI = detprop->DriftVelocity(Efield_SI,Temperature); //drift velocity between shield and induction (cm/us)
163  double driftvelocity_IC = detprop->DriftVelocity(Efield_IC,Temperature); //drift velocity between induction and collection (cm/us)
164  double timepitch = driftvelocity*timetick; //time sample (cm)
165  double tSI = plane_pitch/driftvelocity_SI/timetick; //drift time between Shield and Collection planes (time samples)
166  double tIC = plane_pitch/driftvelocity_IC/timetick; //drift time between Induction and Collection planes (time samples)
167 
168 
169  // get input Cluster object(s).
170  art::Handle< std::vector<recob::Cluster> > clusterListHandle;
171  evt.getByLabel(fClusterModuleLabel,clusterListHandle);
172 
173  // get input EndPoint2D object(s).
174  art::Handle< std::vector<recob::EndPoint2D> > endpointListHandle;
175  evt.getByLabel(fEndPoint2DModuleLabel,endpointListHandle);
176 
178  if(evt.getByLabel(fEndPoint2DModuleLabel,endpointListHandle))
179  for (unsigned int i = 0; i < endpointListHandle->size(); ++i){
180  art::Ptr<recob::EndPoint2D> endpointHolder(endpointListHandle,i);
181  endpointlist.push_back(endpointHolder);
182  }
183 
184  // Declare some vectors..
185  // Induction
186  std::vector<double> Iwirefirsts; // in cm
187  std::vector<double> Iwirelasts; // in cm
188  std::vector<double> Itimefirsts; // in cm
189  std::vector<double> Itimelasts; // in cm
190  std::vector<double> Itimefirsts_line; // in cm
191  std::vector<double> Itimelasts_line; // in cm
192  std::vector < std::vector< art::Ptr<recob::Hit> > > IclusHitlists;
193  std::vector<unsigned int> Icluster_count;
194 
195  // Collection
196  std::vector<double> Cwirefirsts; // in cm
197  std::vector<double> Cwirelasts; // in cm
198  std::vector<double> Ctimefirsts; // in cm
199  std::vector<double> Ctimelasts; // in cm
200  std::vector<double> Ctimefirsts_line; // in cm
201  std::vector<double> Ctimelasts_line; // in cm
202  std::vector< std::vector< art::Ptr<recob::Hit> > > CclusHitlists;
203  std::vector<unsigned int> Ccluster_count;
204 
205  art::FindManyP<recob::Hit> fm(clusterListHandle, evt, fClusterModuleLabel);
206 
207  for(unsigned int ii = 0; ii < clusterListHandle->size(); ++ii){
208 
209  art::Ptr<recob::Cluster> cl(clusterListHandle, ii);
210 
211  // Figure out which View the cluster belongs to
212  //only consider merged-lines that are associated with the vertex.
213  //this helps get rid of through-going muon background -spitz
214  int vtx2d_w = -99999;
215  double vtx2d_t = -99999;
216  bool found2dvtx = false;
217 
218  for (unsigned int j = 0; j<endpointlist.size();j++){
219  if (endpointlist[j]->View() == cl->View()){
220  vtx2d_w = endpointlist[j]->WireID().Wire; //for update to EndPoint2D ... WK 4/22/13
221  vtx2d_t = endpointlist[j]->DriftTime();
222  found2dvtx = true;
223  break;
224  }
225  }
226  if (found2dvtx){
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);
232  if (dis>fvertexclusterWindow) continue;
233  }
234  //else continue; //what to do if a 2D vertex is not found? perhaps vertex finder was not even run.
235 
236  // Some variables for the hit
237  float time; //hit time at maximum
238 
239  std::vector< art::Ptr<recob::Hit> > hitlist = fm.at(ii);
240  std::sort(hitlist.begin(), hitlist.end(), trkf::SortByWire());
241 
242  TGraph *the2Dtrack = new TGraph(hitlist.size());
243 
244  std::vector<double> wires;
245  std::vector<double> times;
246 
247 
248  int np=0;
249  //loop over cluster hits
250  for(std::vector< art::Ptr<recob::Hit> >::const_iterator theHit = hitlist.begin(); theHit != hitlist.end(); theHit++){
251  //recover the Hit
252  // recob::Hit* theHit = (recob::Hit*)(*hitIter);
253  time = (*theHit)->PeakTime() ;
254 
255  time -= presamplings;
256 
257 
258  if(geom->SignalType((*theHit)->Channel()) == geo::kCollection)
259  time -= tIC; // Collection
260  //transform hit wire and time into cm
261  double wire_cm = 0.;
262  if(geom->SignalType((*theHit)->Channel()) == geo::kInduction)
263  wire_cm = (double)(((*theHit)->WireID().Wire+3.95) * wire_pitch);
264  else
265  wire_cm = (double)(((*theHit)->WireID().Wire+1.84) * wire_pitch);
266 
267  //double time_cm = (double)(time * timepitch);
268  double time_cm;
269  if(time>tSI) time_cm = (double)( (time-tSI)*timepitch + tSI*driftvelocity_SI*timetick);
270  else time_cm = time*driftvelocity_SI*timetick;
271 
272  wires.push_back(wire_cm);
273  times.push_back(time_cm);
274 
275  the2Dtrack->SetPoint(np,wire_cm,time_cm);
276  np++;
277  }//end of loop over cluster hits
278 
279  // fit the 2Dtrack and get some info to store
280  try{
281  the2Dtrack->Fit("pol1","Q");
282  }
283  catch(...){
284  std::cout<<"The 2D track fit failed"<<std::endl;
285  continue;
286  }
287 
288  TF1 *pol1=(TF1*) the2Dtrack->GetFunction("pol1");
289  double par[2];
290  pol1->GetParameters(par);
291  double intercept = par[0];
292  double slope = par[1];
293 
294 
295  double w0 = wires.front(); // first hit wire (cm)
296  double w1 = wires.back(); // last hit wire (cm)
297  double t0 = times.front(); // first hit time (cm)
298  double t1 = times.back(); // last hit time (cm)
299  double t0_line = intercept + (w0)*slope;// time coordinate at wire w0 on the fit line (cm)
300  double t1_line = intercept + (w1)*slope;// time coordinate at wire w1 on the fit line (cm)
301 
302 
303 
304  // actually store the 2Dtrack info
305  switch(geom->SignalType((*hitlist.begin())->Channel())){
306  case geo::kInduction:
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);
315  break;
316  case geo::kCollection:
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);
325  break;
326  case geo::kMysteryType:
327  break;
328  }
329  delete pol1;
330  }// end of loop over all input clusters
331 
335 
336  for(unsigned int collectionIter=0; collectionIter < CclusHitlists.size();collectionIter++){ //loop over Collection view 2D tracks
337  // Recover previously stored info
338  double Cw0 = Cwirefirsts[collectionIter];
339  double Cw1 = Cwirelasts[collectionIter];
340  //double Ct0 = Ctimefirsts[collectionIter];
341  //double Ct1 = Ctimelasts[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];
345 
346  double collLength = TMath::Sqrt( TMath::Power(Ct1_line - Ct0_line,2) + TMath::Power(Cw1 - Cw0,2));
347 
348  for(unsigned int inductionIter=0;inductionIter<IclusHitlists.size();inductionIter++){ //loop over Induction view 2D tracks
349  // Recover previously stored info
350  double Iw0 = Iwirefirsts[inductionIter];
351  double Iw1 = Iwirelasts[inductionIter];
352  //double It0 = Itimefirsts[inductionIter];
353  //double It1 = Itimelasts[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];
357 
358  double indLength = TMath::Sqrt( TMath::Power(It1_line - It0_line,2) + TMath::Power(Iw1 - Iw0,2));
359 
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));
362 
363 
364  if(forward_match || backward_match ){
365 
366  // Reconstruct the 3D track
367  TVector3 XYZ0, XYZ1; // track endpoints
368  if(forward_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));
371  }
372  else{
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));
375  }
376 
377  //compute track direction in Local co-ordinate system
378  //WARNING: There is an ambiguity introduced here for the case of backwards-going tracks.
379  //If available, vertex info. could sort this out.
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());
385  if(forward_match){
386  collVtx.Set(Ct0_line,Cw0);
387  indVtx.Set(It0_line,Iw0);
388  }else{
389  collVtx.Set(Ct0_line,Cw0);
390  indVtx.Set(It1_line,Iw1);
391  }
392  }
393  else{
394  startpointVec.SetXYZ(XYZ1.X(),XYZ1.Y(),XYZ1.Z());
395  endpointVec.SetXYZ(XYZ0.X(),XYZ0.Y(),XYZ0.Z());
396  if(forward_match){
397  collVtx.Set(Ct1_line,Cw1);
398  indVtx.Set(It1_line,Iw1);
399  }else{
400  collVtx.Set(Ct1_line,Cw1);
401  indVtx.Set(It0_line,Iw0);
402  }
403  }
404 
405  //compute track (normalized) cosine directions in the TPC co-ordinate system
406  TVector3 DirCos = endpointVec - startpointVec;
407 
408  //SetMag casues a crash if the magnitude of the vector is zero
409  try
410  {
411  DirCos.SetMag(1.0);//normalize vector
412  }
413  catch(...){std::cout<<"The Spacepoint is infinitely small"<<std::endl;
414  continue;
415  }
416 
417  art::Ptr <recob::Cluster> cl1(clusterListHandle,Icluster_count[inductionIter]);
418  art::Ptr <recob::Cluster> cl2(clusterListHandle,Ccluster_count[collectionIter]);
419  art::PtrVector<recob::Cluster> clustersPerTrack;
420  clustersPerTrack.push_back(cl1);
421  clustersPerTrack.push_back(cl2);
422 
423 
425  // Match hits
427 
428  //create collection of spacepoints that will be used when creating the Track object
429  std::vector<recob::SpacePoint> spacepoints;
430 
431 
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;
434 
435 
436  std::vector<bool> maxhitsMatch(maxhits.size());
437  for(unsigned int it=0;it<maxhits.size();it++) maxhitsMatch[it] = false;
438 
439  std::vector<recob::Hit*> hits3Dmatched;
440  // For the matching start from the view where the track projection presents less hits
441  unsigned int imaximum = 0;
442  size_t spStart = spcol->size();
443  for(unsigned int imin=0;imin<minhits.size();imin++){ //loop over hits
444  //get wire - time coordinate of the hit
445  //unsigned int channel,wire,plane1,plane2,tpc,cstat;
446  geo::WireID hit1WireID = minhits[imin]->WireID();
447  auto const hitSigType = minhits[imin]->SignalType();
448  double w1=0;
449 
450  //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.
451  if(hitSigType == geo::kInduction)
452  w1 = (double)((hit1WireID.Wire+3.95) * wire_pitch);
453  else
454  w1 = (double)((hit1WireID.Wire+1.84) * wire_pitch);
455 
456  double temptime1 = minhits[imin]->PeakTime()-presamplings;
457  if(hitSigType == geo::kCollection) temptime1 -= tIC;
458  double t1;// = plane1==1?(double)((minhits[imin]->PeakTime()-presamplings-tIC)*timepitch):(double)((minhits[imin]->PeakTime()-presamplings)*timepitch); //in cm
459  if(temptime1>tSI) t1 = (double)( (temptime1-tSI)*timepitch + tSI*driftvelocity_SI*timetick);
460  else t1 = temptime1*driftvelocity_SI*timetick;
461 
462  //get the track origin co-ordinates in the two views
463  TVector2 minVtx2D;
464  (hitSigType == geo::kCollection) ? minVtx2D.Set(collVtx.X(),collVtx.Y()): minVtx2D.Set(indVtx.X(),indVtx.Y());
465  TVector2 maxVtx2D;
466  (hitSigType == geo::kCollection) ? maxVtx2D.Set(indVtx.X(),indVtx.Y()): maxVtx2D.Set(collVtx.X(),collVtx.Y());
467 
468  double ratio = (collLength>indLength) ? collLength/indLength : indLength/collLength;
469 
470  //compute the distance of the hit (imin) from the relative track origin
471  double minDistance = ratio*TMath::Sqrt(TMath::Power(t1-minVtx2D.X(),2) + TMath::Power(w1-minVtx2D.Y(),2));
472 
473 
474  //core matching algorithm
475  double difference = 9999999.;
476 
477  for(unsigned int imax = 0; imax < maxhits.size(); imax++){ //loop over hits of the other view
478  if(!maxhitsMatch[imax]){
479  //get wire - time coordinate of the hit
480  geo::WireID hit2WireID = maxhits[imax]->WireID();
481  auto const hit2SigType = maxhits[imax]->SignalType();
482  double w2=0.;
483  if(hit2SigType == geo::kInduction)
484  w2 = (double)((hit2WireID.Wire+3.95) * wire_pitch);
485  else
486  w2 = (double)((hit2WireID.Wire+1.84) * wire_pitch);
487 
488  double temptime2 = maxhits[imax]->PeakTime()-presamplings;
489  if(hit2SigType == geo::kCollection) temptime2 -= tIC;
490  double t2;
491  if(temptime2>tSI) t2 = (double)( (temptime2-tSI)*timepitch + tSI*driftvelocity_SI*timetick);
492  else t2 = temptime2*driftvelocity_SI*timetick;
493 
494  bool timematch = (std::abs(t1-t2)<ftmatch*timepitch);
495  bool wirematch = (std::abs(w1-w2)<wireShift*wire_pitch);
496 
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);
500  imaximum = imax;
501  }
502  }
503  }
504  maxhitsMatch[imaximum]=true;
505 
507  if(difference!= 9999999.){
508  sp_hits.push_back(minhits[imin]);
509  sp_hits.push_back(maxhits[imaximum]);
510  }
511 
512  // Get the time-wire co-ordinates of the matched hit
513  geo::WireID hit2WireID = maxhits[imaximum]->WireID();
514  auto const hit2SigType = maxhits[imaximum]->SignalType();
515 
516  //double w1_match = (double)((wire+1)*wire_pitch);
517  double w1_match=0.;
518  if(hit2SigType == geo::kInduction)
519  w1_match = (double)((hit2WireID.Wire+3.95) * wire_pitch);
520  else
521  w1_match = (double)((hit2WireID.Wire+1.84) * wire_pitch);
522 
523  double temptime3 = maxhits[imaximum]->PeakTime()-presamplings;
524  if(hit2SigType == geo::kCollection) temptime3 -= tIC;
525  double t1_match;
526  if(temptime3>tSI) t1_match = (double)( (temptime3-tSI)*timepitch + tSI*driftvelocity_SI*timetick);
527  else t1_match = temptime3*driftvelocity_SI*timetick;
528 
529  // create the 3D hit, compute its co-ordinates and add it to the 3D hits list
530  double Ct = hitSigType==geo::kCollection?t1:t1_match;
531  double Cw = hit2SigType==geo::kCollection?w1:w1_match;
532  double Iw = hit2SigType==geo::kCollection?w1_match:w1;
533 
534  const TVector3 hit3d(Ct,(Cw-Iw)/(2.*TMath::Sin(Angle)),(Cw+Iw)/(2.*TMath::Cos(Angle))-YC/2.*TMath::Tan(Angle));
535 
536 
537  Double_t hitcoord[3];
538  hitcoord[0] = hit3d.X();
539  hitcoord[1] = hit3d.Y();
540  hitcoord[2] = hit3d.Z();
541 
542  /*
543  double yy,zz;
544  if(geom->ChannelsIntersect(geom->PlaneWireToChannel(0,(int)((Iw/wire_pitch)-3.95)), geom->PlaneWireToChannel(1,(int)((Cw/wire_pitch)-1.84)),yy,zz))
545  {
546  //channelsintersect provides a slightly more accurate set of y and z coordinates. use channelsintersect in case the wires in question do cross.
547  hitcoord[1] = yy;
548  hitcoord[2] = zz;
549  mf::LogInfo("SpacePts: ") << "SpacePoint adding xyz ..." << hitcoord[0] <<","<< hitcoord[1] <<","<< hitcoord[2];
550  // std::cout<<"wire 1: "<<(Iw/wire_pitch)-3.95<<" "<<(Cw/wire_pitch)-1.84<<std::endl;
551  // std::cout<<"Intersect: "<<yy<<" "<<zz<<std::endl;
552  }
553  else
554  continue;
555  */
556 
557  double err[6] = {util::kBogusD};
558  recob::SpacePoint mysp(hitcoord, err, util::kBogusD, spStart + spacepoints.size());//3d point at end of track
559  // Don't add a spacepoint right on top of the last one.
560  const double eps(0.1); // 1mm
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) )
568  continue;
569  }
570  spacepoints.push_back(mysp);
571  spcol->push_back(mysp);
572  util::CreateAssn(*this, evt, *spcol, sp_hits, *shassn);
573 
574  }//loop over min-hits
575 
576  size_t spEnd = spcol->size();
577 
578  // Add the 3D track to the vector of the reconstructed tracks
579  if(spacepoints.size()>0){
580 
581  // make a vector of the trajectory points along the track
582  std::vector<TVector3> xyz(spacepoints.size());
583  for(size_t s = 0; s < spacepoints.size(); ++s){
584  xyz[s] = TVector3(spacepoints[s].XYZ());
585  }
586 
588  std::vector<TVector3> dircos(spacepoints.size(), DirCos);
589 
592  recob::Track::Flags_t(xyz.size()), false),
593  0, -1., 0, recob::tracking::SMatrixSym55(), recob::tracking::SMatrixSym55(), tcol->size()));
594 
595  // make associations between the track and space points
596  util::CreateAssn(*this, evt, *tcol, *spcol, *tspassn, spStart, spEnd);
597 
598  // now the track and clusters
599  util::CreateAssn(*this, evt, *tcol, clustersPerTrack, *tcassn);
600 
601  // and the hits and track
602  art::FindManyP<recob::Hit> fmh(clustersPerTrack, evt, fClusterModuleLabel);
603  for(size_t cpt = 0; cpt < clustersPerTrack.size(); ++cpt)
604  util::CreateAssn(*this, evt, *tcol, fmh.at(cpt), *thassn);
605 
606  }
607  } //close match 2D tracks
608 
609  }//close loop over Induction view 2D tracks
610 
611  }//close loop over Collection xxview 2D tracks
612 
613  mf::LogVerbatim("Summary") << std::setfill('-') << std::setw(175) << "-" << std::setfill(' ');
614  mf::LogVerbatim("Summary") << "SpacePts Summary:";
615  for(unsigned int i = 0; i<tcol->size(); ++i) mf::LogVerbatim("Summary") << tcol->at(i) ;
616 
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));
623 
624 } // end SpacePts::produce()
code to link reconstructed objects back to the MC truth information
Float_t s
Definition: plot.C:23
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.
TTree * t1
Definition: plottest35.C:26
Who knows?
Definition: geo_types.h:94
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
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
TrackTrajectory::Flags_t Flags_t
Definition: Track.h:71
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:313
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:192
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
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:441
Signal from induction planes.
Definition: geo_types.h:92
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)
Definition: TrackingTypes.h:55
bool SortByWire(art::Ptr< recob::Hit > const &h1, art::Ptr< recob::Hit > const &h2)
size_type size() const
Definition: PtrVector.h:308
TTree * t2
Definition: plottest35.C:36
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
virtual double DriftVelocity(double efield=0., double temperature=0.) const =0
constexpr double kBogusD
obviously bogus double value
std::string fEndPoint2DModuleLabel
Float_t w
Definition: plot.C:23
Definition: fwd.h:25
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:52
std::string GetLArTPCVolumeName(geo::TPCID const &tpcid) const
Return the name of specified LAr TPC volume.
Signal from collection planes.
Definition: geo_types.h:93
void trkf::SpacePts::reconfigure ( fhicl::ParameterSet const &  p)

Definition at line 107 of file SpacePts_module.cc.

References fClusterModuleLabel, fEndPoint2DModuleLabel, fPreSamplings, ftmatch, fvertexclusterWindow, and fhicl::ParameterSet::get().

Referenced by SpacePts().

108 {
109  fPreSamplings = pset.get< double >("TicksOffset");
110  ftmatch = pset.get< int >("TMatch");
111  fClusterModuleLabel = pset.get< std::string >("ClusterModuleLabel");
112  fEndPoint2DModuleLabel = pset.get< std::string >("EndPoint2DModuleLabel");
113  fvertexclusterWindow = pset.get< double >("vertexclusterWindow");
114 }
double fvertexclusterWindow
std::string fClusterModuleLabel
std::string fEndPoint2DModuleLabel
void art::Consumer::showMissingConsumes ( ) const
protectedinherited

Definition at line 125 of file Consumer.cc.

Referenced by art::EDProducer::doEndJob(), art::EDFilter::doEndJob(), art::EDAnalyzer::doEndJob(), and art::RootOutput::endJob().

126 {
127  if (!moduleContext_)
128  return;
129 
130  // If none of the branches have missing consumes statements, exit early.
131  if (std::all_of(cbegin(missingConsumes_),
132  cend(missingConsumes_),
133  [](auto const& perBranch) { return perBranch.empty(); }))
134  return;
135 
136  constexpr cet::HorizontalRule rule{60};
137  mf::LogPrint log{"MTdiagnostics"};
138  log << '\n'
139  << rule('=') << '\n'
140  << "The following consumes (or mayConsume) statements are missing from\n"
141  << module_context(moduleDescription_) << '\n'
142  << rule('-') << '\n';
143 
144  cet::for_all_with_index(
145  missingConsumes_, [&log](std::size_t const i, auto const& perBranch) {
146  for (auto const& pi : perBranch) {
147  log << " "
148  << assemble_consumes_statement(static_cast<BranchType>(i), pi)
149  << '\n';
150  }
151  });
152  log << rule('=');
153 }
cet::exempt_ptr< ModuleDescription const > moduleDescription_
Definition: Consumer.h:140
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
bool moduleContext_
Definition: Consumer.h:136
ConsumableProductSets missingConsumes_
Definition: Consumer.h:139
void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited

Definition at line 101 of file Consumer.cc.

References art::errors::ProductRegistrationFailure.

103 {
104  // Early exits if consumes tracking has been disabled or if the
105  // consumed product is an allowed consumable.
106  if (!moduleContext_)
107  return;
108 
109  if (cet::binary_search_all(consumables_[bt], pi))
110  return;
111 
112  if (requireConsumes_) {
114  "Consumer: an error occurred during validation of a "
115  "retrieved product\n\n")
116  << "The following consumes (or mayConsume) statement is missing from\n"
117  << module_context(moduleDescription_) << ":\n\n"
118  << " " << assemble_consumes_statement(bt, pi) << "\n\n";
119  }
120 
121  missingConsumes_[bt].insert(pi);
122 }
cet::exempt_ptr< ModuleDescription const > moduleDescription_
Definition: Consumer.h:140
bool requireConsumes_
Definition: Consumer.h:137
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
ConsumableProductSets missingConsumes_
Definition: Consumer.h:139

Member Data Documentation

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

Definition at line 71 of file SpacePts_module.cc.

Referenced by produce(), and reconfigure().

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

Definition at line 72 of file SpacePts_module.cc.

Referenced by produce(), and reconfigure().

double trkf::SpacePts::fPreSamplings
private

Definition at line 69 of file SpacePts_module.cc.

Referenced by produce(), and reconfigure().

int trkf::SpacePts::ftmatch
private

Definition at line 68 of file SpacePts_module.cc.

Referenced by produce(), and reconfigure().

double trkf::SpacePts::fvertexclusterWindow
private

Definition at line 70 of file SpacePts_module.cc.

Referenced by produce(), and reconfigure().


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