TwoCRTMatchingProducer_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: TwoCRTMatchingProducer
3 // Plugin Type: producer (art v2_10_03)
4 // File: TwoCRTMatchingProducer_module.cc
5 // Author: Richie Diurba (rdiurba@fnal.gov)
6 // Tingjun Yang (tjyang@fnal.gov)
7 //
8 // Generated at Wed Jun 27 04:09:39 2018 by Andrew Olivier using cetskelgen
9 // from cetlib version v3_02_00.
10 ////////////////////////////////////////////////////////////////////////
11 
12 //Framework includes
18 
26 #include "fhiclcpp/ParameterSet.h"
29 #include "canvas/Persistency/Common/FindManyP.h"
30 #include "art_root_io/TFileService.h"
31 
32 //LArSoft includes
33 
37 #include "nug4/ParticleNavigation/ParticleList.h"
38 
39 
48 
49 
52 
55 
57 
58 
59 //Local includes
61 
62 
63 
64 //ROOT includes
65 #include "TH1.h"
66 #include "TH2.h"
67 #include "TCanvas.h"
68 #include "TImage.h"
69 #include "TTree.h"
70 #include "TH1D.h"
71 #include "TStyle.h"
72 #include "TString.h"
73 
74 //c++ includes
75 #include <numeric> //std::accumulate was moved from <algorithm> to <numeric> in c++14
76 #include <iostream>
77 #include <cmath>
78 using namespace std; // Namespaces established to make life easier
79 using namespace ROOT::Math;
80 
81 namespace CRT {
82  class TwoCRTMatchingProducer;
83 }
84 
86 public:
88  // The compiler-generated destructor is fine for non-base
89  // classes without bare pointers or other resource use.
90 
91  // Plugins should not be copied or assigned.
94  TwoCRTMatchingProducer & operator = (TwoCRTMatchingProducer const &) = delete;
95  TwoCRTMatchingProducer & operator = (TwoCRTMatchingProducer &&) = delete;
96  bool moduleMatcher(int module1, int module2);
97  // Required functions.
98 
99  void produce(art::Event & event) override;
100 
101  int nEvents = 0;
102  int nHitsPerEvent=0;
103  std::string fTrackModuleLabel = "pandoraTrack";
104 
105 
106 private:
107  art::InputTag fCRTLabel; //Label for the module that analyzed
113  long long timeStamp;
114 
115 
116 typedef struct // Structures for arrays to move hits from raw to reco to validation
117  {
119  int channel;
120  int module;
121  int adc;
123  }
124  tempHits;
125 
126  typedef struct {
127  int tempId;
130  double hitPositionX;
131  double hitPositionY;
132  double hitPositionZ;
133  double timeAvg;
134  }
135  recoHits;
136 
137 
138  std::vector < recoHits > primaryHits_F;
139  std::vector < recoHits > primaryHits_B;
140 
141  std::vector < tempHits > tempHits_F;
142  std::vector < tempHits > tempHits_B;
143 
144 
145 
146 };
147 
148 
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 }
162 
163 // v6 Geo Channel Map
164 bool CRT::TwoCRTMatchingProducer::moduleMatcher(int module1, int module2) {
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 }
170 
171 
172 //Turn sim::AuxDetSimChannels into CRT::Hits.
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 }
587 
588 
589 
590 
code to link reconstructed objects back to the MC truth information
std::string string
Definition: nybbler.cc:12
bool moduleMatcher(int module1, int module2)
STL namespace.
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
Particle class.
art framework interface to geometry description
T abs(T value)
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
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
TwoCRTMatchingProducer(fhicl::ParameterSet const &p)
p
Definition: test.py:223
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 produce(art::Event &event) override
Declaration of signal hit object.
void End(void)
Definition: gXSecComp.cxx:210
def center(depos, point)
Definition: depos.py:117
Provides recob::Track data product.
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
Event finding and building.