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

Classes

struct  recoHits
 
struct  tempHits
 

Public Member Functions

 TwoCRTMatchingProducer (fhicl::ParameterSet const &p)
 
 TwoCRTMatchingProducer (TwoCRTMatchingProducer const &)=delete
 
 TwoCRTMatchingProducer (TwoCRTMatchingProducer &&)=delete
 
TwoCRTMatchingProduceroperator= (TwoCRTMatchingProducer const &)=delete
 
TwoCRTMatchingProduceroperator= (TwoCRTMatchingProducer &&)=delete
 
bool moduleMatcher (int module1, int module2)
 
void produce (art::Event &event) override
 
- 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)
 

Public Attributes

int nEvents = 0
 
int nHitsPerEvent =0
 
std::string fTrackModuleLabel = "pandoraTrack"
 

Private Attributes

art::InputTag fCRTLabel
 
bool fSCECorrection
 
bool fModuleSwitch
 
int fADCThreshold
 
int fModuletoModuleTimingCut
 
int fFronttoBackTimingCut
 
long long timeStamp
 
std::vector< recoHitsprimaryHits_F
 
std::vector< recoHitsprimaryHits_B
 
std::vector< tempHitstempHits_F
 
std::vector< tempHitstempHits_B
 

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 85 of file TwoCRTMatchingProducer_module.cc.

Constructor & Destructor Documentation

CRT::TwoCRTMatchingProducer::TwoCRTMatchingProducer ( fhicl::ParameterSet const &  p)
explicit

Definition at line 149 of file TwoCRTMatchingProducer_module.cc.

149  : EDProducer{p}, fCRTLabel(p.get < art::InputTag > ("CRTLabel"))
150 {
151  consumes < std::vector < CRT::Trigger >> (fCRTLabel);
152  consumes < std::vector < art::Assns < sim::AuxDetSimChannel, CRT::Trigger >>> (fCRTLabel);
153  produces< std::vector<anab::T0> >();
154 
155  produces< std::vector<anab::CosmicTag> >();
156  produces< art::Assns<recob::Track, anab::T0> >();
157  produces< art::Assns<recob::Track, anab::CosmicTag> >();
158  produces< art::Assns<CRT::Trigger, anab::CosmicTag> >();
159 
160  fSCECorrection=(p.get<bool>("SCECorrection"));
161 }
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
p
Definition: test.py:223
CRT::TwoCRTMatchingProducer::TwoCRTMatchingProducer ( TwoCRTMatchingProducer const &  )
delete
CRT::TwoCRTMatchingProducer::TwoCRTMatchingProducer ( TwoCRTMatchingProducer &&  )
delete

Member Function Documentation

bool CRT::TwoCRTMatchingProducer::moduleMatcher ( int  module1,
int  module2 
)

Definition at line 164 of file TwoCRTMatchingProducer_module.cc.

164  {
165  // Function checking if two hits could reasonably be matched into a 2D hit
166  if ((module1 == 6 && (module2 == 10 || module2 == 11)) || (module1 == 14 && (module2 == 10 || module2 == 11)) || (module1 == 19 && (module2 == 26 || module2 == 27)) || (module1 == 31 && (module2 == 26 || module2 == 27)) || (module1 == 7 && (module2 == 12 || module2 == 13)) || (module1 == 15 && (module2 == 12 || module2 == 13)) || (module1 == 18 && (module2 == 24 || module2 == 25)) || (module1 == 30 && (module2 == 24 || module2 == 25)) || (module1 == 1 && (module2 == 4 || module2 == 5)) || (module1 == 9 && (module2 == 4 || module2 == 5)) || (module1 == 16 && (module2 == 20 || module2 == 21)) || (module1 == 28 && (module2 == 20 || module2 == 21)) || (module1 == 0 && (module2 == 2 || module2 == 3)) || (module1 == 8 && (module2 == 2 || module2 == 3)) || (module1 == 17 && (module2 == 22 || module2 == 23)) || (module1 == 29 && (module2 == 22 || module2 == 23))) return 1;
167  else return 0;
168 
169 }
TwoCRTMatchingProducer& CRT::TwoCRTMatchingProducer::operator= ( TwoCRTMatchingProducer const &  )
delete
TwoCRTMatchingProducer& CRT::TwoCRTMatchingProducer::operator= ( TwoCRTMatchingProducer &&  )
delete
void CRT::TwoCRTMatchingProducer::produce ( art::Event event)
overridevirtual

Implements art::EDProducer.

Definition at line 173 of file TwoCRTMatchingProducer_module.cc.

174 {
175  bool fMCCSwitch=!event.isRealData();
176  nEvents++;
177  std::unique_ptr< std::vector<anab::T0> > T0col( new std::vector<anab::T0>);
178  auto CRTTrack=std::make_unique< std::vector< anab::CosmicTag > > ();
179 
180 
181  std::unique_ptr< art::Assns<CRT::Trigger, anab::CosmicTag>> CRTTriggerassn( new art::Assns<CRT::Trigger, anab::CosmicTag>);
182 
183  std::unique_ptr< art::Assns<recob::Track, anab::CosmicTag> > TPCCRTassn( new art::Assns<recob::Track, anab::CosmicTag>);
184  std::unique_ptr< art::Assns<recob::Track, anab::T0> > TPCT0assn( new art::Assns<recob::Track, anab::T0>);
185 
186 
187  if (fMCCSwitch){
188  fModuleSwitch=1;
189  fADCThreshold=800;
192 
193 }
194  else {
195  fModuleSwitch=0;
196  fADCThreshold=10;
199  art::ValidHandle<std::vector<raw::RDTimeStamp>> timingHandle = event.getValidHandle<std::vector<raw::RDTimeStamp>>("timingrawdecoder:daq");
200  timeStamp=timingHandle->at(0).GetTimeStamp();
201 }
202  int nHits = 0;
203 
204  //Detector properties service
205  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(event);
206 
207  primaryHits_F.clear();
208  primaryHits_B.clear();
209  tempHits_F.clear();
210  tempHits_B.clear(); // Arrays to compile hits and move them through
211 
212 
213  //Get triggers
214  //cout << "Getting triggers" << endl;
215  vector < art::Ptr < CRT::Trigger > > crtList;
216  auto crtListHandle = event.getHandle < vector < CRT::Trigger > >(fCRTLabel);
217  if (crtListHandle) {
218  art::fill_ptr_vector(crtList, crtListHandle);
219  }
220  const auto & triggers = event.getValidHandle < std::vector < CRT::Trigger >> (fCRTLabel);
221 
222  art::FindManyP < sim::AuxDetSimChannel > trigToSim(triggers, event, fCRTLabel);
223 
224  //Get a handle to the Geometry service to look up AuxDetGeos from module numbers
226 
227  auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
228 
229  //Mapping from channel to trigger
230  std::unordered_map < size_t, double > prevTimes;
231  int hitID = 0;
232  //cout << "Looking for hits in Triggers" << endl;
233  int trigID=0;
234  for (const auto & trigger: * triggers) {
235  const auto & hits = trigger.Hits();
236  for (const auto & hit: hits) { // Collect hits on all modules
237  //cout<<hits.size()<<','<<hit.ADC()<<endl;
238  if (hit.ADC() > fADCThreshold) { // Keep if they are above threshold
239 
240  tempHits tHits;
241  if (!fMCCSwitch){
242 
243  tHits.module = trigger.Channel(); // Values to add to array
244  tHits.channel=hit.Channel();
245  tHits.adc = hit.ADC();
246  tHits.triggerTime=trigger.Timestamp()-timeStamp;
247  }
248  else{
249  tHits.module = trigger.Channel(); // Values to add to array
250  tHits.channel=hit.Channel();
251  tHits.adc = hit.ADC();
252  tHits.triggerTime=trigger.Timestamp();
253  }
254  //cout<<trigger.Channel()<<','<<hit.Channel()<<','<<hit.ADC()<<endl;
255  nHits++;
256  tHits.triggerNumber=trigID;
257  const auto & trigGeo = geom -> AuxDet(trigger.Channel()); // Get geo
258  const auto & csens = trigGeo.SensitiveVolume(hit.Channel());
259  const auto center = csens.GetCenter();
260  if (center.Z() < 100) tempHits_F.push_back(tHits); // Sort F/B from Z
261  else tempHits_B.push_back(tHits);
262  hitID++;
263  }
264  }
265  trigID++;
266  }
267  nHitsPerEvent=nHits;
268  //cout << "Hits compiled for event: " << nEvents << endl;
269  //cout << "Number of Hits above Threshold: " << hitID << endl;
270 
271  for (unsigned int f = 0; f < tempHits_F.size(); f++) {
272  for (unsigned int f_test = 0; f_test < tempHits_F.size(); f_test++) {
273  if (fabs(tempHits_F[f_test].triggerTime-tempHits_F[f].triggerTime)>fModuletoModuleTimingCut) continue;
274  const auto & trigGeo = geom -> AuxDet(tempHits_F[f].module);
275  const auto & trigGeo2 = geom -> AuxDet(tempHits_F[f_test].module);
276 
277  const auto & hit1Geo = trigGeo.SensitiveVolume(tempHits_F[f].channel);
278  const auto hit1Center = hit1Geo.GetCenter();
279  // Create 2D hits from geo of the Y and X modules
280  const auto & hit2Geo = trigGeo2.SensitiveVolume(tempHits_F[f_test].channel);
281  const auto hit2Center = hit2Geo.GetCenter();
282  bool moduleMatched;
283  moduleMatched=moduleMatcher(tempHits_F[f_test].module, tempHits_F[f].module);
284  if (moduleMatched) {
285  // Get the center of the hits (CRT_Res=2.5 cm)
286  double hitX = hit1Center.X();
287  for (unsigned int a = 0; a < tempHits_F.size(); a++)
288  {
289  if(tempHits_F[a].module==tempHits_F[f].module && (tempHits_F[a].channel-1)==tempHits_F[f].channel) hitX=hit1Center.X()+1.25;
290  }
291  double hitYPrelim=hit2Center.Y();
292  for (unsigned int a = 0; a < tempHits_F.size(); a++)
293  {
294  if(tempHits_F[a].module==tempHits_F[f_test].module && (tempHits_F[a].channel-1)==tempHits_F[f_test].channel) hitYPrelim=hit2Center.Y()+1.25;
295  }
296 
297 
298 
299  double hitY=hitYPrelim;
300  double hitZ = (hit1Center.Z() + hit2Center.Z()) / 2.f;
301 
302  recoHits rHits;
303  rHits.hitPositionX = hitX;
304  rHits.hitPositionY = hitY;
305  rHits.hitPositionZ = hitZ;
306  rHits.trigNumberX=tempHits_F[f].triggerNumber;
307  rHits.trigNumberY=tempHits_F[f_test].triggerNumber;
308  rHits.timeAvg = (tempHits_F[f_test].triggerTime+tempHits_F[f].triggerTime)/2.0;
309  primaryHits_F.push_back(rHits); // Add array
310  }
311  }
312  }
313  for (unsigned int f = 0; f < tempHits_B.size(); f++) {
314  for (unsigned int f_test = 0; f_test < tempHits_B.size(); f_test++) { // Same as above but for back CRT
315  if (fabs(tempHits_B[f_test].triggerTime-tempHits_B[f].triggerTime)>fModuletoModuleTimingCut) continue;
316 
317  const auto & trigGeo = geom -> AuxDet(tempHits_B[f].module);
318  const auto & trigGeo2 = geom -> AuxDet(tempHits_B[f_test].module);
319  const auto & hit1Geo = trigGeo.SensitiveVolume(tempHits_B[f].channel);
320  const auto hit1Center = hit1Geo.GetCenter();
321 
322  const auto & hit2Geo = trigGeo2.SensitiveVolume(tempHits_B[f_test].channel);
323  const auto hit2Center = hit2Geo.GetCenter();
324  bool moduleMatched;
325  moduleMatched=moduleMatcher(tempHits_B[f_test].module, tempHits_B[f].module);
326 
327  if (moduleMatched) {
328  double hitX = hit1Center.X();
329 
330 
331  for (unsigned int a = 0; a < tempHits_B.size(); a++)
332  {
333  if(tempHits_B[a].module==tempHits_B[f].module && (tempHits_B[a].channel-1)==tempHits_B[f].channel) hitX=hit1Center.X()+1.25;
334  }
335 
336  double hitYPrelim = hit2Center.Y();
337 
338  for (unsigned int a = 0; a < tempHits_B.size(); a++)
339  {
340  if(tempHits_B[a].module==tempHits_B[f_test].module && (tempHits_B[a].channel-1)==tempHits_B[f_test].channel) hitYPrelim=hit2Center.Y()+1.25;
341  }
342  double hitY=hitYPrelim;
343 
344 
345  double hitZ = (hit1Center.Z() + hit2Center.Z()) / 2.f;
346 
347  recoHits rHits;
348  rHits.hitPositionX = hitX;
349  rHits.hitPositionY = hitY;
350  rHits.hitPositionZ = hitZ;
351  rHits.trigNumberX=tempHits_B[f].triggerNumber;
352  rHits.trigNumberY=tempHits_B[f_test].triggerNumber;
353  rHits.timeAvg = (tempHits_B[f_test].triggerTime+tempHits_B[f].triggerTime)/2.0;
354  primaryHits_B.push_back(rHits);
355 
356  }
357  }
358  }
359  vector < art::Ptr < recob::Track > > trackList;
360  auto trackListHandle = event.getHandle < vector < recob::Track > >(fTrackModuleLabel);
361  if (trackListHandle) {
362  art::fill_ptr_vector(trackList, trackListHandle);
363  }
364  else{
365  event.put(std::move(T0col)); event.put(std::move(CRTTrack)); event.put(std::move(TPCCRTassn)); event.put(std::move(TPCT0assn)); event.put(std::move(CRTTriggerassn));
366  return;
367  }
368 
369  vector<art::Ptr<recob::PFParticle> > pfplist;
370  auto PFPListHandle = event.getHandle< std::vector<recob::PFParticle> >("pandora");
371  if (PFPListHandle) art::fill_ptr_vector(pfplist, PFPListHandle);
372  if(pfplist.size()<1) return;
373  art::FindManyP<anab::T0> trk_t0_assn_v(PFPListHandle, event ,"pandora");
374  art::FindManyP<recob::PFParticle> pfp_trk_assn(trackListHandle,event,"pandoraTrack");
375  int nTracksReco = trackList.size();
376  art::FindManyP < recob::Hit > hitsFromTrack(trackListHandle, event, fTrackModuleLabel);
377 
378  for (int iRecoTrack = 0; iRecoTrack < nTracksReco; ++iRecoTrack) {
379  std::vector< art::Ptr<recob::Hit> > allHits = hitsFromTrack.at(iRecoTrack);
380 
381  art::Ptr<recob::Track> ptrack(trackListHandle, iRecoTrack);
382 
383  std::vector<art::Ptr<recob::PFParticle>> pfps=pfp_trk_assn.at(iRecoTrack);
384  if(!pfps.size()) continue;
385  std::vector<art::Ptr<anab::T0>> t0s=trk_t0_assn_v.at(pfps[0].key());
386 
387 
388  double trackStartPositionZ_noSCE = trackList[iRecoTrack]->Vertex().Z();
389  double trackEndPositionZ_noSCE = trackList[iRecoTrack] -> End().Z();
390 
391  double trackStartPositionX_notCorrected = trackList[iRecoTrack]->Vertex().X();
392  double trackStartPositionY_noSCE = trackList[iRecoTrack]->Vertex().Y();
393 
394 
395  double trackEndPositionX_notCorrected = trackList[iRecoTrack] -> End().X();
396  double trackEndPositionY_noSCE = trackList[iRecoTrack] -> End().Y();
397 
398  int firstHit=0;
399  int lastHit=allHits.size()-2;
400  if (trackStartPositionZ_noSCE>trackEndPositionZ_noSCE){
401  trackEndPositionZ_noSCE = trackList[iRecoTrack]->Vertex().Z();
402  trackStartPositionZ_noSCE = trackList[iRecoTrack] -> End().Z();
403  trackEndPositionX_notCorrected = trackList[iRecoTrack]->Vertex().X();
404  trackEndPositionY_noSCE = trackList[iRecoTrack]->Vertex().Y();
405 
406 
407  trackStartPositionX_notCorrected=trackList[iRecoTrack] -> End().X();
408  trackStartPositionY_noSCE = trackList[iRecoTrack] -> End().Y();
409  firstHit=lastHit;
410  lastHit=0;
411  }
412 
413 
414 
415 
416 
417 
418 
419  if ((trackEndPositionZ_noSCE > 660 && trackStartPositionZ_noSCE < 50) || (trackStartPositionZ_noSCE > 660 && trackEndPositionZ_noSCE < 50)) {
420 
421  double min_delta = DBL_MAX;
422  double best_XF = DBL_MAX;
423  double best_YF = DBL_MAX;
424  double best_ZF = DBL_MAX;
425  double best_XB = DBL_MAX;
426  double best_YB = DBL_MAX;
427  double best_ZB = DBL_MAX;
428  double best_dotProductCos = DBL_MAX;
429  double best_deltaXF = DBL_MAX;
430  double best_deltaYF = DBL_MAX;
431  double best_deltaXB = DBL_MAX;
432  double best_deltaYB = DBL_MAX;
433  double best_T=DBL_MAX;
434  int best_trigXF=0;
435  int best_trigYF=0;
436  int best_trigXB=0;
437  int best_trigYB=0;
438 
439  for (unsigned int f = 0; f < primaryHits_F.size(); f++) {
440  for (unsigned int b = 0; b < primaryHits_B.size(); b++) {
441 
442  double X1 = primaryHits_F[f].hitPositionX;
443  double Y1 = primaryHits_F[f].hitPositionY;
444  double Z1 = primaryHits_F[f].hitPositionZ;
445  double X2 = primaryHits_B[b].hitPositionX;
446  double Y2 = primaryHits_B[b].hitPositionY;
447  double Z2= primaryHits_B[b].hitPositionZ;
448 
449 
450  if (fabs(primaryHits_F[f].timeAvg-primaryHits_B[b].timeAvg)>fFronttoBackTimingCut) continue;
451  double t0=(primaryHits_F[f].timeAvg+primaryHits_B[b].timeAvg)/2.f;
452  int tempId = 0;
453  double xOffset=0;
454 
455 
456  double trackStartPositionX_noSCE=trackStartPositionX_notCorrected;
457  double trackEndPositionX_noSCE=trackEndPositionX_notCorrected;
458  if (t0s.empty()){
459  int RDOffset=0;
460  if (!fMCCSwitch) RDOffset=111;
461  double ticksOffset=0;
462 
463  if (!fMCCSwitch) ticksOffset = (t0+RDOffset)/25.f+detProp.GetXTicksOffset(allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
464 
465  else if (fMCCSwitch) ticksOffset = (t0/500.f)+detProp.GetXTicksOffset(allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
466 
467  xOffset=detProp.ConvertTicksToX(ticksOffset,allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
468  //double xOffset=.08*ticksOffset
469 
470  trackStartPositionX_noSCE=trackStartPositionX_notCorrected-xOffset;
471  trackEndPositionX_noSCE=trackEndPositionX_notCorrected-xOffset;
472  }
473 
474  double trackStartPositionX=trackStartPositionX_noSCE;
475  double trackStartPositionY=trackStartPositionY_noSCE;
476  double trackStartPositionZ=trackStartPositionZ_noSCE;
477 
478  double trackEndPositionX=trackEndPositionX_noSCE;
479  double trackEndPositionY=trackEndPositionY_noSCE;
480  double trackEndPositionZ=trackEndPositionZ_noSCE;
481  if (fSCECorrection && SCE->EnableCalSpatialSCE()){
482 if(geom->PositionToTPCID(geo::Point_t(trackEndPositionX, trackEndPositionY, trackEndPositionZ)).deepestIndex()<13 && geom->PositionToTPCID(geo::Point_t(trackStartPositionX, trackStartPositionY, trackStartPositionZ)).deepestIndex()<13){
483  auto const & posOffsets_F = SCE->GetCalPosOffsets(geo::Point_t(trackStartPositionX, trackStartPositionY, trackStartPositionZ), geom->PositionToTPCID(geo::Point_t(trackStartPositionX, trackStartPositionY, trackStartPositionZ)).deepestIndex());
484  trackStartPositionX -= posOffsets_F.X();
485  trackStartPositionY += posOffsets_F.Y();
486  trackStartPositionZ += posOffsets_F.Z();
487  auto const & posOffsets_B = SCE->GetCalPosOffsets(geo::Point_t(trackEndPositionX, trackEndPositionY, trackEndPositionZ), geom->PositionToTPCID(geo::Point_t(trackEndPositionX, trackEndPositionY, trackEndPositionZ)).deepestIndex());
488  trackEndPositionX -= posOffsets_B.X();
489  trackEndPositionY += posOffsets_B.Y();
490  trackEndPositionZ += posOffsets_B.Z();
491  }
492  }
493 
494 
495 
496 
497 
498  // Make metrics for a CRT pair to compare later
499  TVector3 trackStart(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
500  TVector3 trackEnd(trackEndPositionX, trackEndPositionY, trackEndPositionZ);
501  TVector3 v1(X1,Y1,Z1);
502  TVector3 v2(X2, Y2, Z2);
503 
504  TVector3 v4(trackStartPositionX,
505  trackStartPositionY,
506  trackStartPositionZ);
507  TVector3 v5(trackEndPositionX,
508  trackEndPositionY,
509  trackEndPositionZ);
510  TVector3 trackVector = (v5-v4).Unit();
511  TVector3 hitVector=(v2-v1).Unit();
512 
513 
514 
515 
516  double predictedHitPositionY1 = (v1.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.Y()-v5.Y())+v5.Y();
517  double predictedHitPositionY2 = (v2.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.Y()-v5.Y())+v5.Y();
518 
519  double predictedHitPositionX1 = (v1.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.X()-v5.X())+v5.X();
520  double predictedHitPositionX2 = (v2.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.X()-v5.X())+v5.X();
521 
522  double dotProductCos=trackVector*hitVector;
523 
524  double deltaX1 = (predictedHitPositionX1-X1);
525  double deltaX2 = (predictedHitPositionX2-X2);
526 
527  double deltaY1 = (predictedHitPositionY1-Y1);
528 
529  double deltaY2 = (predictedHitPositionY2-Y2);
530 
531 
532 
533  tempId++;
534 
535  //iRecoTrack
536  if (min_delta > std::abs(deltaX1)+std::abs(deltaX2) + std::abs(deltaY1)+std::abs(deltaY2) ){
537 
538  min_delta=std::abs(deltaX1)+std::abs(deltaX2) + std::abs(deltaY1)+std::abs(deltaY2);
539  best_XF = X1;
540  best_YF = Y1;
541  best_ZF = Z1;
542  best_XB = X2;
543  best_YB = Y2;
544  best_ZB = Z2;
545  best_dotProductCos = dotProductCos;
546  best_deltaXF = deltaX1;
547  best_deltaYF = deltaY1;
548  best_deltaXB = deltaX2;
549  best_deltaYB = deltaY2;
550  best_trigXF=primaryHits_F[f].trigNumberX;
551  best_trigYF=primaryHits_F[f].trigNumberY;
552  best_trigXB=primaryHits_B[b].trigNumberX;
553  best_trigYB=primaryHits_B[b].trigNumberY;
554  best_T = t0;
555  if (!fMCCSwitch) best_T=(111.f+best_T)*20.f;
556  // Added 111 tick CRT-CTB offset
557  }
558  }
559  }
560  if (std::abs(best_dotProductCos)>0.99 && std::abs(best_deltaXF)+std::abs(best_deltaXB)<40 && std::abs(best_deltaYF)+std::abs(best_deltaYB)<40 ) {
561  //std::cout<<"Found match with TPC*CRT "<<best_dotProductCos<<std::endl;
562 
563 //std::cout<<"Displacement in front and back "<<best_deltaXF<<","<<best_deltaYF<<","<<best_deltaXB<<","<<best_deltaYB<<std::endl;
564  std::vector<float> hitF;
565  std::vector<float> hitB;
566  hitF.push_back(best_XF); hitF.push_back(best_YF); hitF.push_back(best_ZF);
567  hitB.push_back(best_XB); hitB.push_back(best_YB); hitB.push_back(best_ZB);
568 
569  CRTTrack->push_back(anab::CosmicTag(hitF,hitB, std::abs(best_dotProductCos),anab::CosmicTagID_t::kNotTagged));
570 
571  T0col->push_back(anab::T0(best_T, 13,2,iRecoTrack,best_dotProductCos));
572  util::CreateAssn(*this, event, *CRTTrack, trackList[iRecoTrack], *TPCCRTassn);
573  util::CreateAssn(*this, event, *T0col, trackList[iRecoTrack], *TPCT0assn);
574  util::CreateAssn(*this, event, *CRTTrack, crtList[best_trigXF], *CRTTriggerassn);
575  util::CreateAssn(*this, event, *CRTTrack, crtList[best_trigYF], *CRTTriggerassn);
576 
577  util::CreateAssn(*this, event, *CRTTrack, crtList[best_trigXB], *CRTTriggerassn);
578 
579  util::CreateAssn(*this, event, *CRTTrack, crtList[best_trigYB], *CRTTriggerassn);
580 
581  }
582  }
583  }
584 
585  event.put(std::move(T0col)); event.put(std::move(CRTTrack)); event.put(std::move(TPCCRTassn)); event.put(std::move(TPCT0assn)); event.put(std::move(CRTTriggerassn));
586 }
code to link reconstructed objects back to the MC truth information
bool moduleMatcher(int module1, int module2)
geo::TPCID PositionToTPCID(geo::Point_t const &point) const
Returns the ID of the TPC at specified location.
uint8_t channel
Definition: CRTFragment.hh:201
Definition: T0.h:16
T abs(T value)
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
def key(type, name=None)
Definition: graph.py:13
const double a
def move(depos, offset)
Definition: depos.py:107
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.
Detector simulation of raw signals on wires.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:184
void End(void)
Definition: gXSecComp.cxx:210
def center(depos, point)
Definition: depos.py:117
static bool * b
Definition: config.cpp:1043
constexpr auto const & deepestIndex() const
Returns the value of the deepest ID available (TPC&#39;s).
Definition: geo_types.h:428
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
recob::tracking::Plane Plane
Definition: TrackState.h:17

Member Data Documentation

int CRT::TwoCRTMatchingProducer::fADCThreshold
private

Definition at line 110 of file TwoCRTMatchingProducer_module.cc.

art::InputTag CRT::TwoCRTMatchingProducer::fCRTLabel
private

Definition at line 107 of file TwoCRTMatchingProducer_module.cc.

int CRT::TwoCRTMatchingProducer::fFronttoBackTimingCut
private

Definition at line 112 of file TwoCRTMatchingProducer_module.cc.

bool CRT::TwoCRTMatchingProducer::fModuleSwitch
private

Definition at line 109 of file TwoCRTMatchingProducer_module.cc.

int CRT::TwoCRTMatchingProducer::fModuletoModuleTimingCut
private

Definition at line 111 of file TwoCRTMatchingProducer_module.cc.

bool CRT::TwoCRTMatchingProducer::fSCECorrection
private

Definition at line 108 of file TwoCRTMatchingProducer_module.cc.

std::string CRT::TwoCRTMatchingProducer::fTrackModuleLabel = "pandoraTrack"

Definition at line 103 of file TwoCRTMatchingProducer_module.cc.

int CRT::TwoCRTMatchingProducer::nEvents = 0

Definition at line 101 of file TwoCRTMatchingProducer_module.cc.

int CRT::TwoCRTMatchingProducer::nHitsPerEvent =0

Definition at line 102 of file TwoCRTMatchingProducer_module.cc.

std::vector< recoHits > CRT::TwoCRTMatchingProducer::primaryHits_B
private

Definition at line 139 of file TwoCRTMatchingProducer_module.cc.

std::vector< recoHits > CRT::TwoCRTMatchingProducer::primaryHits_F
private

Definition at line 138 of file TwoCRTMatchingProducer_module.cc.

std::vector< tempHits > CRT::TwoCRTMatchingProducer::tempHits_B
private

Definition at line 142 of file TwoCRTMatchingProducer_module.cc.

std::vector< tempHits > CRT::TwoCRTMatchingProducer::tempHits_F
private

Definition at line 141 of file TwoCRTMatchingProducer_module.cc.

long long CRT::TwoCRTMatchingProducer::timeStamp
private

Definition at line 113 of file TwoCRTMatchingProducer_module.cc.


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