SingleCRTMatchingProducer_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: SingleCRTMatchingProducer
3 // Plugin Type: producer (art v2_10_03)
4 // File: SingleCRTMatchingProducer_module.cc
5 //
6 // Generated at Wed Jun 27 04:09:39 2018 by Andrew Olivier using cetskelgen
7 // from cetlib version v3_02_00.
8 ////////////////////////////////////////////////////////////////////////
9 
10 //Framework includes
23 #include "fhiclcpp/ParameterSet.h"
26 #include "canvas/Persistency/Common/FindManyP.h"
27 #include "art_root_io/TFileService.h"
29 
30 //LArSoft includes
31 
35 #include "nug4/ParticleNavigation/ParticleList.h"
36 
37 
46 
47 
48 
50 
53 
54 //Local includes
56 
57 //ROOT includes
58 #include "TH1.h"
59 #include "TH2.h"
60 #include "TCanvas.h"
61 #include "TImage.h"
62 #include "TTree.h"
63 #include "TH1D.h"
64 #include "TStyle.h"
65 #include "TString.h"
66 
67 //c++ includes
68 #include <numeric> //std::accumulate was moved from <algorithm> to <numeric> in c++14
69 #include <iostream>
70 #include <cmath>
71 using namespace std; // Namespaces established to make life easier
72 using namespace ROOT::Math;
73 namespace CRT {
74  class SingleCRTMatchingProducer;
75 }
76 
77 
79 public:
81 
82  // The compiler-generated destructor is fine for non-base
83  // classes without bare pointers or other resource use.
84 
85  // Plugins should not be copied or assigned.
88  SingleCRTMatchingProducer & operator = (SingleCRTMatchingProducer const &) = delete;
89  SingleCRTMatchingProducer & operator = (SingleCRTMatchingProducer &&) = delete;
90 
91  bool moduleMatcher(int module1, int module2);
92  void produce(art::Event & event) override;
93  std::string fTrackModuleLabel = "pandoraTrack";
94 
95  int nEvents = 0;
96  int nHaloMuons = 0;
97  private: ofstream logFile;
98  const std::string reco="reco";
99  //Parameters for reading in CRT::Triggers and associated AuxDetSimChannels.
100  art::InputTag fCRTLabel; //Label for the module that produced
108 
109 
110  double dotCos;
111  int adcX, adcY;
112  double X_CRT, Y_CRT, Z_CRT;
113  double trackX1, trackX2, trackY1, trackY2, trackZ1, trackZ2;
114  int moduleX, moduleY;
115  int stripX, stripY;
116  double deltaX, deltaY;
117  double CRTT0;
118  double flashTime;
119  double opCRTTDiff;
120  typedef struct // Structures for arrays to move hits from raw to reco to validation
121  {
122 
123  int channel;
124  int module;
126  int adc;
129  }
130  tempHits;
131 
132  typedef struct {
133  int tempId;
134  int adcX;
135  int adcY;
136  int stripX, stripY, moduleX, moduleY;
137  double hitPositionX;
138  double hitPositionY;
139  double hitPositionZ;
140  double timeAvg;
141  int trigNumberX, trigNumberY;
142  }
143  recoHits;
144 
145  typedef struct // These are displacement metrics for track and hit reco
146  {
147  int tempId;
149  int recoId;
150  int adcX1;
151  int adcY1;
152  double deltaX;
153  double deltaY;
155  double X1;
156  double Y1;
157  double Z1;
158  int trackID;
161  int moduleX1, moduleY1;
162  int stripX1, stripY1;
163  double flashTDiff;
164  double timeAvg;
165  int trigNumberX, trigNumberY;
166  }
167  tracksPair;
168 
169  struct removePairIndex // Struct to remove tracks after sorting
170  {
172  removePairIndex(const tracksPair & tracksPair0): tracksPair1(tracksPair0) {}
173 
174  bool operator()(const tracksPair & tracksPair2) {
175  return (tracksPair1.recoId == tracksPair2.recoId || tracksPair1.CRTTrackId == tracksPair2.CRTTrackId || (tracksPair1.stripX1==tracksPair2.stripX1 && tracksPair1.moduleX1==tracksPair2.moduleX1) || (tracksPair1.stripY1==tracksPair2.stripY1 && tracksPair1.moduleY1==tracksPair2.moduleY1));
176  }
177  };
178 
179  struct sortPair // Struct to sort to find best CRT track for TPC track
180  {
181  bool operator()(const tracksPair & pair1,
182  const tracksPair & pair2) {
183 
184  //return (fabs(pair1.dotProductCos)>fabs(pair2.dotProductCos));
185  return (fabs(pair1.deltaX)+fabs(pair1.deltaY)<fabs(pair2.deltaX)+fabs(pair2.deltaY));
186 
187  }
188  };
189  std::vector < recoHits > primaryHits_F;
190  std::vector < recoHits > primaryHits_B;
191 
192  std::vector < tempHits > tempHits_F;
193  std::vector < tempHits > tempHits_B;
194  std::vector < tracksPair > tracksPair_F;
195  std::vector < tracksPair > tracksPair_B;
196 };
197 
199  const & p):
200  EDProducer{p}, fCRTLabel(p.get < art::InputTag > ("CRTLabel")) {
201  consumes < std::vector < CRT::Trigger >> (fCRTLabel);
202  consumes < std::vector < art::Assns < sim::AuxDetSimChannel, CRT::Trigger >>> (fCRTLabel); // CRT art consumables
203  fMCCSwitch=(p.get<bool>("MCC"));
204  produces< std::vector<anab::T0> >();
205  produces< std::vector<anab::CosmicTag> >();
206  produces< art::Assns<recob::Track, anab::T0> >();
207  produces< art::Assns<recob::Track, anab::CosmicTag> >();
208  produces< art::Assns<anab::CosmicTag, anab::T0> >();
209  produces< art::Assns<CRT::Trigger, anab::CosmicTag> >();
210  fSCECorrection=(p.get<bool>("SCECorrection"));
211  }
212 
213 
214 
215 // v6 Geo Channel Map
216 bool CRT::SingleCRTMatchingProducer::moduleMatcher(int module1, int module2) {
217  // Function checking if two hits could reasonably be matched into a 2D hit
218  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;
219  else return 0;
220 
221 }
222 
223 
224 //Turn sim::AuxDetSimChannels into CRT::Hits.
226 {
227 
228  std::unique_ptr< std::vector<anab::T0> > T0col( new std::vector<anab::T0>);
229 
230  auto CRTTrack=std::make_unique< std::vector< anab::CosmicTag > > ();
231 
232  std::unique_ptr< art::Assns<anab::CosmicTag, anab::T0> > CRTT0assn( new art::Assns<anab::CosmicTag, anab::T0>);
233 
234  std::unique_ptr< art::Assns<recob::Track, anab::CosmicTag> > TPCCRTassn( new art::Assns<recob::Track, anab::CosmicTag>);
235  std::unique_ptr< art::Assns<recob::Track, anab::T0> > TPCT0assn( new art::Assns<recob::Track, anab::T0>);
236 
237  std::unique_ptr< art::Assns<CRT::Trigger, anab::CosmicTag>> CRTTriggerassn( new art::Assns<CRT::Trigger, anab::CosmicTag>);
238 
239  if (fMCCSwitch){
240  fModuleSwitch=1;
241  fADCThreshold=800;
244  fOpCRTTDiffCut=200;
245 
246 }
247  else {
248  fModuleSwitch=0;
249  fADCThreshold=10;
252  fOpCRTTDiffCut=1000000;
253 
254 
255 }
256 
257 /* if(!fMCCSwitch){
258  art::ValidHandle<std::vector<raw::RDTimeStamp>> timingHandle = event.getValidHandle<std::vector<raw::RDTimeStamp>>("timingrawdecoder:daq");
259  //const auto& pdspctbs = event.getValidHandle<std::vector<raw::ctb::pdspctb>>(fCTB_tag);
260 
261  //const raw::RDTimeStamp& timeStamp = timingHandle->at(0);
262  //if(timeStamp.GetFlags()!= 13) {event.put(std::move(CRTTrack)); event.put(std::move(T0col)); event.put(std::move(TPCCRTassn)); event.put(std::move(CRTT0assn)); event.put(std::move(TPCT0assn)); return;}
263  }*/
264 
265  int nHits = 0;
268 
269  primaryHits_F.clear();
270  primaryHits_B.clear();
271  tracksPair_F.clear();
272  tracksPair_B.clear();
273  tempHits_F.clear();
274  tempHits_B.clear(); // Arrays to compile hits and move them through
275  primaryHits_F.clear();
276  //allTracksPair.clear();
277  logFile.open("ProtoDUNE.log"); // Logfile I don't use right now
278 
279  //Get triggers
280  //cout << "Getting triggers" << endl;
281  const auto & triggers = event.getValidHandle < std::vector < CRT::Trigger >> (fCRTLabel);
282 
283  art::FindManyP < sim::AuxDetSimChannel > trigToSim(triggers, event, fCRTLabel);
284 
285  vector < art::Ptr < CRT::Trigger > > crtList;
286  auto crtListHandle = event.getHandle < vector < CRT::Trigger > >(fCRTLabel);
287  if (crtListHandle) {
288  art::fill_ptr_vector(crtList, crtListHandle);
289  }
290 
291  //Get a handle to the Geometry service to look up AuxDetGeos from module numbers
293  auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
294 
295  //Mapping from channel to trigger
296  std::unordered_map < size_t, double > prevTimes;
297  int hitID = 0;
298  int trigID=0;
299  for (const auto & trigger: * triggers) {
300  const auto & hits = trigger.Hits();
301  for (const auto & hit: hits) { // Collect hits on all modules
302  if (hit.ADC() > fADCThreshold) { // Keep if they are above threshold
303 
304  tempHits tHits;
305  if (!fMCCSwitch){
306  art::ValidHandle<std::vector<raw::RDTimeStamp>> timingHandle = event.getValidHandle<std::vector<raw::RDTimeStamp>>("timingrawdecoder:daq");
307 
308  tHits.module = trigger.Channel(); // Values to add to array
309  tHits.channelGeo = hit.Channel();
310  tHits.channel=hit.Channel();
311  tHits.adc = hit.ADC();
312  tHits.triggerTime=trigger.Timestamp()-timingHandle->at(0).GetTimeStamp();
313  tHits.triggerNumber=trigID;
314  }
315  else{
316  tHits.module = trigger.Channel(); // Values to add to array
317  tHits.channelGeo = hit.Channel();
318  tHits.channel=hit.Channel();
319  tHits.adc = hit.ADC();
320  tHits.triggerTime=trigger.Timestamp();
321  tHits.triggerNumber=trigID;
322  }
323  //cout<<trigger.Channel()<<','<<hit.Channel()<<','<<hit.ADC()<<endl;
324  nHits++;
325 
326  const auto & trigGeo = geom -> AuxDet(trigger.Channel()); // Get geo
327  const auto & csens = trigGeo.SensitiveVolume(hit.Channel());
328  const auto center = csens.GetCenter();
329  if (center.Z() < 100) tempHits_F.push_back(tHits); // Sort F/B from Z
330  else tempHits_B.push_back(tHits);
331  hitID++;
332  }
333  }
334  trigID++;
335  }
336 
337  //cout << "Hits compiled for event: " << nEvents << endl;
338  //cout << "Number of Hits above Threshold: " << hitID << endl;
339 
340  for (unsigned int f = 0; f < tempHits_F.size(); f++) {
341  for (unsigned int f_test = 0; f_test < tempHits_F.size(); f_test++) {
342  if (fabs(tempHits_F[f_test].triggerTime-tempHits_F[f].triggerTime)>fModuletoModuleTimingCut) continue;
343  const auto & trigGeo = geom -> AuxDet(tempHits_F[f].module);
344  const auto & trigGeo2 = geom -> AuxDet(tempHits_F[f_test].module);
345 
346  const auto & hit1Geo = trigGeo.SensitiveVolume(tempHits_F[f].channelGeo);
347  const auto hit1Center = hit1Geo.GetCenter();
348  // Create 2D hits from geo of the Y and X modules
349  const auto & hit2Geo = trigGeo2.SensitiveVolume(tempHits_F[f_test].channelGeo);
350  const auto hit2Center = hit2Geo.GetCenter();
351  bool moduleMatched;
352  moduleMatched=moduleMatcher(tempHits_F[f_test].module, tempHits_F[f].module);
353 
354  if (moduleMatched) {
355  // Get the center of the hits (CRT_Res=2.5 cm)
356  double hitX = hit1Center.X();
357  for (unsigned int a = 0; a < tempHits_F.size(); a++)
358  {
359  if(tempHits_F[a].module==tempHits_F[f].module && (tempHits_F[a].channelGeo-1)==tempHits_F[f].channelGeo) hitX=hit1Center.X()+1.25;
360  }
361  double hitYPrelim=hit2Center.Y();
362  for (unsigned int a = 0; a < tempHits_F.size(); a++)
363  {
364  if(tempHits_F[a].module==tempHits_F[f_test].module && (tempHits_F[a].channelGeo-1)==tempHits_F[f_test].channelGeo) hitYPrelim=hit2Center.Y()+1.25;
365  }
366 
367 
368 
369  double hitY=hitYPrelim;
370  double hitZ = (hit1Center.Z() + hit2Center.Z()) / 2.f;
371 
372  recoHits rHits;
373  rHits.adcX=tempHits_F[f].adc;
374  rHits.adcY=tempHits_F[f_test].adc;
375  rHits.hitPositionX = hitX;
376  rHits.hitPositionY = hitY;
377  rHits.hitPositionZ = hitZ;
378  rHits.moduleX=tempHits_F[f].module;
379  rHits.moduleY=tempHits_F[f_test].module;
380  rHits.trigNumberX=tempHits_F[f].triggerNumber;
381  rHits.trigNumberY=tempHits_F[f_test].triggerNumber;
382  rHits.stripX=tempHits_F[f].channel;
383  rHits.stripY=tempHits_F[f_test].channel;
384  rHits.timeAvg = (tempHits_F[f_test].triggerTime+tempHits_F[f].triggerTime)/2.0;
385  primaryHits_F.push_back(rHits); // Add array
386  }
387  }
388  }
389  for (unsigned int f = 0; f < tempHits_B.size(); f++) {
390  for (unsigned int f_test = 0; f_test < tempHits_B.size(); f_test++) { // Same as above but for back CRT
391  if (fabs(tempHits_B[f_test].triggerTime-tempHits_B[f].triggerTime)>fModuletoModuleTimingCut) continue;
392 
393  const auto & trigGeo = geom -> AuxDet(tempHits_B[f].module);
394  const auto & trigGeo2 = geom -> AuxDet(tempHits_B[f_test].module);
395  const auto & hit1Geo = trigGeo.SensitiveVolume(tempHits_B[f].channelGeo);
396  const auto hit1Center = hit1Geo.GetCenter();
397 
398  const auto & hit2Geo = trigGeo2.SensitiveVolume(tempHits_B[f_test].channelGeo);
399  const auto hit2Center = hit2Geo.GetCenter();
400  bool moduleMatched;
401  moduleMatched=moduleMatcher(tempHits_B[f_test].module, tempHits_B[f].module);
402 
403 
404  if (moduleMatched) {
405  double hitX = hit1Center.X();
406 
407 
408  for (unsigned int a = 0; a < tempHits_B.size(); a++)
409  {
410  if(tempHits_B[a].module==tempHits_B[f].module && (tempHits_B[a].channelGeo-1)==tempHits_B[f].channelGeo) hitX=hit1Center.X()+1.25;
411  }
412 
413  double hitYPrelim = hit2Center.Y();
414 
415  for (unsigned int a = 0; a < tempHits_B.size(); a++)
416  {
417  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;
418  }
419  double hitY=hitYPrelim;
420 
421 
422  double hitZ = (hit1Center.Z() + hit2Center.Z()) / 2.f;
423 
424  recoHits rHits;
425  rHits.adcX=tempHits_B[f].adc;
426  rHits.adcY=tempHits_B[f_test].adc;
427  rHits.hitPositionX = hitX;
428  rHits.hitPositionY = hitY;
429  rHits.hitPositionZ = hitZ;
430  rHits.moduleX=tempHits_B[f].module;
431  rHits.moduleY=tempHits_B[f_test].module;
432  rHits.stripX=tempHits_B[f].channel;
433  rHits.trigNumberX=tempHits_B[f].triggerNumber;
434  rHits.trigNumberY=tempHits_B[f_test].triggerNumber;
435  rHits.stripY=tempHits_B[f_test].channel;
436  rHits.timeAvg = (tempHits_B[f_test].triggerTime+tempHits_B[f].triggerTime)/2.0;
437  primaryHits_B.push_back(rHits);
438 
439  }
440  }
441  }
442 
443  auto const t0CandPtr = art::PtrMaker<anab::T0>(event);
444  auto const crtPtr = art::PtrMaker<anab::CosmicTag>(event);
445  // Reconstruciton information
446 
447 
448  vector < art::Ptr < recob::Track > > trackList;
449  auto trackListHandle = event.getHandle < vector < recob::Track > >(fTrackModuleLabel);
450  if (trackListHandle) {
451  art::fill_ptr_vector(trackList, trackListHandle);
452  }
453 
454  vector<art::Ptr<recob::PFParticle> > pfplist;
455  auto PFPListHandle = event.getHandle< std::vector<recob::PFParticle> >("pandora");
456  if (PFPListHandle) art::fill_ptr_vector(pfplist, PFPListHandle);
457  if (pfplist.size()<1) return;
458  art::FindManyP<anab::T0> trk_t0_assn_v(PFPListHandle, event ,"pandora");
459  art::FindManyP<recob::PFParticle> pfp_trk_assn(trackListHandle,event,"pandoraTrack");
460 
461  int nTracksReco = trackList.size();
462  art::FindManyP < recob::Hit > hitsFromTrack(trackListHandle, event, fTrackModuleLabel);
463  int tempId = 0;
464 
465  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(event);
466 
467  for (int iRecoTrack = 0; iRecoTrack < nTracksReco; ++iRecoTrack) {
468  if (primaryHits_F.size()+primaryHits_B.size()<1) break;
469  std::vector< art::Ptr<recob::Hit> > allHits = hitsFromTrack.at(iRecoTrack);
470 
471  art::Ptr<recob::Track> ptrack(trackListHandle, iRecoTrack);
472 
473  std::vector<art::Ptr<recob::PFParticle>> pfps=pfp_trk_assn.at(iRecoTrack);
474  if(!pfps.size()) continue;
475  std::vector<art::Ptr<anab::T0>> t0s=trk_t0_assn_v.at(pfps[0].key());
476  if(t0s.size()){
477  //auto t0=t0s.at(0);
478  //int t_zero=t0->Time();
479  //cout<<"Pandora T0: "<<t_zero<<endl;
480  }
481  int firstHit=0;
482  int lastHit=allHits.size()-2;
483 
484 
485 // Get track positions and find angles
486  // if (fMCCSwith==1){
487  double trackStartPositionZ_noSCE = trackList[iRecoTrack]->Vertex().Z();
488  double trackEndPositionZ_noSCE = trackList[iRecoTrack] -> End().Z();
489 
490  double trackStartPositionX_noSCE = trackList[iRecoTrack]->Vertex().X();
491  double trackStartPositionY_noSCE = trackList[iRecoTrack]->Vertex().Y();
492 
493 
494  double trackEndPositionX_noSCE = trackList[iRecoTrack] -> End().X();
495  double trackEndPositionY_noSCE = trackList[iRecoTrack] -> End().Y();
496 
497  if (trackStartPositionZ_noSCE>trackEndPositionZ_noSCE){
498  trackEndPositionZ_noSCE = trackList[iRecoTrack]->Vertex().Z();
499  trackStartPositionZ_noSCE = trackList[iRecoTrack] -> End().Z();
500  trackEndPositionX_noSCE = trackList[iRecoTrack]->Vertex().X();
501  trackEndPositionY_noSCE = trackList[iRecoTrack]->Vertex().Y();
502 
503 
504  trackStartPositionX_noSCE = trackList[iRecoTrack] -> End().X();
505  trackStartPositionY_noSCE = trackList[iRecoTrack] -> End().Y();
506  firstHit=lastHit;
507  lastHit=0;
508  }
509  if ((trackEndPositionZ_noSCE>90 && trackEndPositionZ_noSCE < 660 && trackStartPositionZ_noSCE <50 && trackStartPositionZ_noSCE<660) || (trackStartPositionZ_noSCE>90 && trackStartPositionZ_noSCE < 660 && trackEndPositionZ_noSCE <50 && trackEndPositionZ_noSCE<660)) {
510 
511  double min_delta = DBL_MAX;
512 
513  double best_dotProductCos = DBL_MAX;
514  double best_deltaXF = DBL_MAX;
515  double best_deltaYF = DBL_MAX;
516  int bestHitIndex_F=0;
517  double best_trackX1=DBL_MAX;
518  double best_trackX2=DBL_MAX;
519  for (unsigned int iHit_F = 0; iHit_F < primaryHits_F.size(); iHit_F++) {
520  double xOffset=0;
521 
522  double trackStartPositionX_notCorrected=trackStartPositionX_noSCE;
523  double trackEndPositionX_notCorrected=trackEndPositionX_noSCE;
524  if (!t0s.empty()){
525  if (event.isRealData() && fabs(t0s.at(0)->Time()-(primaryHits_F[iHit_F].timeAvg*20.f))>100000) continue;
526  if (!event.isRealData() && fabs(t0s.at(0)->Time()-primaryHits_F[iHit_F].timeAvg)>100000) continue;
527  }
528  if (t0s.empty()){
529  int RDOffset=0;
530  if (!fMCCSwitch) RDOffset=111;
531  double ticksOffset=0;
532  //cout<<(primaryHits_F[iHit_F].timeAvg+RDOffset)<<endl;
533  //cout<<detProp.GetXTicksOffset(allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat)<<endl;
534  if (!fMCCSwitch) ticksOffset = (primaryHits_F[iHit_F].timeAvg+RDOffset)/25.f+detProp.GetXTicksOffset(allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
535 
536  else if (fMCCSwitch) ticksOffset = (primaryHits_F[iHit_F].timeAvg/500.f)+detProp.GetXTicksOffset(allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
537 
538  xOffset=detProp.ConvertTicksToX(ticksOffset,allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
539  //double xOffset=.08*ticksOffset
540 
541  trackStartPositionX_noSCE=trackStartPositionX_notCorrected-xOffset;
542  trackEndPositionX_noSCE=trackEndPositionX_notCorrected-xOffset;
543  if (fabs(xOffset)>300 || ((trackStartPositionX_notCorrected<0 && trackStartPositionX_noSCE>0) || (trackEndPositionX_notCorrected<0 && trackEndPositionX_noSCE>0)) || ((trackStartPositionX_notCorrected>0 && trackStartPositionX_noSCE<0) || (trackEndPositionX_notCorrected>0 && trackEndPositionX_noSCE<0)) ) continue;
544  }
545 
546  double trackStartPositionX=trackStartPositionX_noSCE;
547  double trackStartPositionY=trackStartPositionY_noSCE;
548  double trackStartPositionZ=trackStartPositionZ_noSCE;
549 
550  double trackEndPositionX=trackEndPositionX_noSCE;
551  double trackEndPositionY=trackEndPositionY_noSCE;
552  double trackEndPositionZ=trackEndPositionZ_noSCE;
553 
554  if (fSCECorrection && SCE->EnableCalSpatialSCE()){
555 if(geom->PositionToTPCID(geo::Point_t(trackEndPositionX, trackEndPositionY, trackEndPositionZ)).deepestIndex()<13 && geom->PositionToTPCID(geo::Point_t(trackStartPositionX, trackStartPositionY, trackStartPositionZ)).deepestIndex()<13){
556  auto const & posOffsets_F = SCE->GetCalPosOffsets(geo::Point_t(trackStartPositionX, trackStartPositionY, trackStartPositionZ), geom->PositionToTPCID(geo::Point_t(trackStartPositionX, trackStartPositionY, trackStartPositionZ)).deepestIndex());
557  trackStartPositionX -= posOffsets_F.X();
558  trackStartPositionY += posOffsets_F.Y();
559  trackStartPositionZ += posOffsets_F.Z();
560  auto const & posOffsets_B = SCE->GetCalPosOffsets(geo::Point_t(trackEndPositionX, trackEndPositionY, trackEndPositionZ), geom->PositionToTPCID(geo::Point_t(trackEndPositionX, trackEndPositionY, trackEndPositionZ)).deepestIndex());
561  trackEndPositionX -= posOffsets_B.X();
562  trackEndPositionY += posOffsets_B.Y();
563  trackEndPositionZ += posOffsets_B.Z();
564  }
565  }
566 
567  double X1 = primaryHits_F[iHit_F].hitPositionX;
568 
569  double Y1 = primaryHits_F[iHit_F].hitPositionY;
570 
571  double Z1 = primaryHits_F[iHit_F].hitPositionZ;
572 
573  // Make metrics for a CRT pair to compare later
574  TVector3 trackStart(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
575  TVector3 trackEnd(trackEndPositionX, trackEndPositionY, trackEndPositionZ);
576  TVector3 v1(X1,Y1,Z1);
577  TVector3 v2(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
578 
579  TVector3 v4(trackStartPositionX,
580  trackStartPositionY,
581  trackStartPositionZ);
582  TVector3 v5(trackEndPositionX,
583  trackEndPositionY,
584  trackEndPositionZ);
585  TVector3 trackVector = (v5-v4).Unit();
586  TVector3 hitVector=(v2-v1).Unit();
587 
588 
589 
590 
591  double predictedHitPositionY1 = (v1.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.Y()-v5.Y())+v5.Y();
592 
593 
594  double predictedHitPositionX1 = (v1.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.X()-v5.X())+v5.X();
595 
596 
597  if (predictedHitPositionX1<-200 || predictedHitPositionX1>580 || predictedHitPositionY1<-50 || predictedHitPositionY1>620) continue;
598 
599  double dotProductCos=trackVector*hitVector;
600 
601  double deltaX1 = (predictedHitPositionX1-X1);
602 
603 
604  double deltaY1 = (predictedHitPositionY1-Y1);
605 
606 
607  if (min_delta > std::abs(deltaX1) + std::abs(deltaY1) ){
608 
609  min_delta=std::abs(deltaX1)+std::abs(deltaY1);
610 
611  best_dotProductCos = dotProductCos;
612  best_trackX1=trackStartPositionX;
613  best_trackX2=trackEndPositionX;
614  best_deltaXF = deltaX1;
615  best_deltaYF = deltaY1;
616  bestHitIndex_F=iHit_F;
617  }
618 
619  }
620  int f=bestHitIndex_F;
621 
622  double deltaX=best_deltaXF; double deltaY=best_deltaYF;
623  double dotProductCos=best_dotProductCos;
624 
625  double trackStartPositionX=best_trackX1;
626  double trackStartPositionY=trackStartPositionY_noSCE;
627  double trackStartPositionZ=trackStartPositionZ_noSCE;
628 
629  double trackEndPositionX=best_trackX2;
630  double trackEndPositionY=trackEndPositionY_noSCE;
631  double trackEndPositionZ=trackEndPositionZ_noSCE;
632 
633 
634  TVector3 trackStart(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
635  TVector3 trackEnd(trackEndPositionX, trackEndPositionY, trackEndPositionZ);
636 
637  tracksPair tPair;
638  tPair.tempId = tempId;
639  tPair.CRTTrackId = f;
640  tPair.recoId = iRecoTrack;
641 
642  tPair.deltaX=deltaX;
643  tPair.deltaY=deltaY;
644  tPair.dotProductCos=dotProductCos;
645 
646  tPair.moduleX1 = primaryHits_F[f].moduleX;
647  tPair.moduleY1 = primaryHits_F[f].moduleY;
648 
649  tPair.adcX1=primaryHits_F[f].adcX;
650  tPair.adcY1=primaryHits_F[f].adcY;
651 
652  tPair.stripX1 = primaryHits_F[f].stripX;
653  tPair.stripY1 = primaryHits_F[f].stripY;
654  tPair.trigNumberX = primaryHits_F[f].trigNumberX;
655  tPair.trigNumberY = primaryHits_F[f].trigNumberY;
656  tPair.X1 = primaryHits_F[f].hitPositionX;
657  tPair.Y1 = primaryHits_F[f].hitPositionY;
658  tPair.Z1 = primaryHits_F[f].hitPositionZ;
659  tPair.timeAvg=primaryHits_F[f].timeAvg;
660  tPair.trackStartPosition=trackStart;
661  tPair.trackEndPosition=trackEnd;
662  tracksPair_B.push_back(tPair);
663 
664  }
665 
666  if ( (trackStartPositionZ_noSCE<620 && trackEndPositionZ_noSCE > 660 && trackStartPositionZ_noSCE > 50 && trackEndPositionZ_noSCE > 50) || (trackStartPositionZ_noSCE>660 && trackEndPositionZ_noSCE < 620 && trackStartPositionZ_noSCE > 50 && trackEndPositionZ_noSCE > 50)) {
667  double min_delta = DBL_MAX;
668 
669  double best_dotProductCos = DBL_MAX;
670  double best_deltaXF = DBL_MAX;
671  double best_deltaYF = DBL_MAX;
672  int bestHitIndex_B=0;
673  double best_trackX1=DBL_MAX;
674  double best_trackX2=DBL_MAX;
675  for (unsigned int iHit_B = 0; iHit_B < primaryHits_B.size(); iHit_B++) {
676 double xOffset=0;
677 
678 
679  double trackStartPositionX_notCorrected=trackStartPositionX_noSCE;
680  double trackEndPositionX_notCorrected=trackEndPositionX_noSCE;
681  if (!t0s.empty()){
682  if (event.isRealData() && fabs(t0s.at(0)->Time()-(primaryHits_B[iHit_B].timeAvg*20.f))>100000) continue;
683  if (!event.isRealData() && fabs(t0s.at(0)->Time()-primaryHits_B[iHit_B].timeAvg)>100000) continue;
684  }
685  if (t0s.empty()){
686  int RDOffset=0;
687  if (!fMCCSwitch) RDOffset=111;
688  double ticksOffset=0;
689  if (!fMCCSwitch) ticksOffset = (primaryHits_B[iHit_B].timeAvg+RDOffset)/25.f+detProp.GetXTicksOffset(allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
690 
691  else if (fMCCSwitch) ticksOffset = (primaryHits_B[iHit_B].timeAvg/500.f)+detProp.GetXTicksOffset(allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
692 
693  xOffset=detProp.ConvertTicksToX(ticksOffset,allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
694  //double xOffset=.08*ticksOffset
695 
696  trackStartPositionX_noSCE=trackStartPositionX_notCorrected-xOffset;
697  trackEndPositionX_noSCE=trackEndPositionX_notCorrected-xOffset;
698 
699  if (fabs(xOffset)>300 || ((trackStartPositionX_notCorrected<0 && trackStartPositionX_noSCE>0) || (trackEndPositionX_notCorrected<0 && trackEndPositionX_noSCE>0)) || ((trackStartPositionX_notCorrected>0 && trackStartPositionX_noSCE<0) || (trackEndPositionX_notCorrected>0 && trackEndPositionX_noSCE<0)) ) continue;
700  }
701 
702  double trackStartPositionX=trackStartPositionX_noSCE;
703  double trackStartPositionY=trackStartPositionY_noSCE;
704  double trackStartPositionZ=trackStartPositionZ_noSCE;
705 
706  double trackEndPositionX=trackEndPositionX_noSCE;
707  double trackEndPositionY=trackEndPositionY_noSCE;
708  double trackEndPositionZ=trackEndPositionZ_noSCE;
709  double X1 = primaryHits_B[iHit_B].hitPositionX;
710 
711  double Y1 = primaryHits_B[iHit_B].hitPositionY;
712 
713  double Z1 = primaryHits_B[iHit_B].hitPositionZ;
714 
715 
716 
717  // Make metrics for a CRT pair to compare later
718  TVector3 trackStart(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
719  TVector3 trackEnd(trackEndPositionX, trackEndPositionY, trackEndPositionZ);
720  TVector3 v1(X1,Y1,Z1);
721  TVector3 v2(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
722 
723  TVector3 v4(trackStartPositionX,
724  trackStartPositionY,
725  trackStartPositionZ);
726  TVector3 v5(trackEndPositionX,
727  trackEndPositionY,
728  trackEndPositionZ);
729  TVector3 trackVector = (v5-v4).Unit();
730  TVector3 hitVector=(v2-v1).Unit();
731 
732 
733 
734 
735  double predictedHitPositionY1 = (v1.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.Y()-v5.Y())+v5.Y();
736 
737 
738  double predictedHitPositionX1 = (v1.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.X()-v5.X())+v5.X();
739 
740  if (abs(predictedHitPositionX1)>340 || predictedHitPositionY1<-160 || predictedHitPositionY1>560) continue;
741 
742  double dotProductCos=trackVector*hitVector;
743 
744  double deltaX1 = (predictedHitPositionX1-X1);
745 
746  double deltaY1 = (predictedHitPositionY1-Y1);
747 
748  if (min_delta > std::abs(deltaX1) + std::abs(deltaY1) ){
749 
750  min_delta=std::abs(deltaX1)+std::abs(deltaY1);
751 
752  best_dotProductCos = dotProductCos;
753  best_trackX1=trackStartPositionX;
754  best_trackX2=trackEndPositionX;
755  best_deltaXF = deltaX1;
756  best_deltaYF = deltaY1;
757  bestHitIndex_B=iHit_B;
758  }
759 
760  }
761  int f=bestHitIndex_B;
762 
763  double deltaX=best_deltaXF; double deltaY=best_deltaYF;
764  double dotProductCos=best_dotProductCos;
765 
766 
767  double trackStartPositionX=best_trackX1;
768  double trackStartPositionY=trackStartPositionY_noSCE;
769  double trackStartPositionZ=trackStartPositionZ_noSCE;
770 
771  double trackEndPositionX=best_trackX2;
772  double trackEndPositionY=trackEndPositionY_noSCE;
773  double trackEndPositionZ=trackEndPositionZ_noSCE;
774 
775  TVector3 trackStart(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
776  TVector3 trackEnd(trackEndPositionX, trackEndPositionY, trackEndPositionZ);
777 
778 
779  tracksPair tPair;
780  tPair.tempId = tempId;
781  tPair.CRTTrackId = f;
782  tPair.recoId = iRecoTrack;
783 
784  tPair.deltaX=deltaX;
785  tPair.deltaY=deltaY;
786  tPair.dotProductCos=dotProductCos;
787 
788  tPair.moduleX1 = primaryHits_B[f].moduleX;
789  tPair.moduleY1 = primaryHits_B[f].moduleY;
790 
791  tPair.adcX1=primaryHits_B[f].adcX;
792  tPair.adcY1=primaryHits_B[f].adcY;
793 
794  tPair.stripX1 = primaryHits_B[f].stripX;
795  tPair.stripY1 = primaryHits_B[f].stripY;
796  tPair.trigNumberX = primaryHits_B[f].trigNumberX;
797  tPair.trigNumberY = primaryHits_B[f].trigNumberY;
798  tPair.X1 = primaryHits_B[f].hitPositionX;
799  tPair.Y1 = primaryHits_B[f].hitPositionY;
800  tPair.Z1 = primaryHits_B[f].hitPositionZ;
801  tPair.timeAvg=primaryHits_B[f].timeAvg;
802  tPair.trackStartPosition=trackStart;
803  tPair.trackEndPosition=trackEnd;
804  tracksPair_B.push_back(tPair);
805 
806 
807  tempId++;
808  } //iRecoTrack
809  }
810 
811  //Sort pair by ascending order of absolute distance
812  sort(tracksPair_F.begin(), tracksPair_F.end(), sortPair());
813  sort(tracksPair_B.begin(), tracksPair_B.end(), sortPair());
814 
815 
816  vector < tracksPair > allUniqueTracksPair;
817  while (tracksPair_F.size()) {
818  allUniqueTracksPair.push_back(tracksPair_F.front());
819  tracksPair_F.erase(remove_if(tracksPair_F.begin(), tracksPair_F.end(), removePairIndex(tracksPair_F.front())),
820  tracksPair_F.end());
821  }
822 
823  while (tracksPair_B.size()) {
824  allUniqueTracksPair.push_back(tracksPair_B.front());
825  tracksPair_B.erase(remove_if(tracksPair_B.begin(), tracksPair_B.end(), removePairIndex(tracksPair_B.front())),
826  tracksPair_B.end());
827  }
828 
829  //cout<<"Number of reco and CRT pairs: "<<allUniqueTracksPair.size()<<endl;
830  if (allUniqueTracksPair.size() > 0) {
831  for (unsigned int u = 0; u < allUniqueTracksPair.size(); u++) {
832 
833  int CRTTrackId=allUniqueTracksPair[u].CRTTrackId;
834  int TPCTrackId=allUniqueTracksPair[u].recoId;
835  deltaX=allUniqueTracksPair[u].deltaX;
836 
837  deltaY=allUniqueTracksPair[u].deltaY;
838  opCRTTDiff=allUniqueTracksPair[u].flashTDiff;
839 
840  dotCos=fabs(allUniqueTracksPair[u].dotProductCos);
841  trackX1=allUniqueTracksPair[u].trackStartPosition.X();
842  trackY1=allUniqueTracksPair[u].trackStartPosition.Y();
843  trackZ1=allUniqueTracksPair[u].trackStartPosition.Z();
844 
845  trackX2=allUniqueTracksPair[u].trackEndPosition.X();
846  trackY2=allUniqueTracksPair[u].trackEndPosition.Y();
847  trackZ2=allUniqueTracksPair[u].trackEndPosition.Z();
848 
849  moduleX=allUniqueTracksPair[u].moduleX1;
850  moduleY=allUniqueTracksPair[u].moduleY1;
851 
852  adcX=allUniqueTracksPair[u].adcX1;
853  adcY=allUniqueTracksPair[u].adcY1;
854 
855  CRTT0=allUniqueTracksPair[u].timeAvg;
856 
857  if (!fMCCSwitch) CRTT0=(111.f+allUniqueTracksPair[u].timeAvg)*20.f;
858  // Added 111 tick CRT-CTB Timing Offset
859  stripX=allUniqueTracksPair[u].stripX1;
860  stripY=allUniqueTracksPair[u].stripY1;
861 
862  X_CRT=allUniqueTracksPair[u].X1;
863  Y_CRT=allUniqueTracksPair[u].Y1;
864  Z_CRT=allUniqueTracksPair[u].Z1;
865 
867  if ( fabs(trackX1)<400 && fabs(trackX2)<400 && fabs(deltaX)<60 && fabs(deltaY)<60 && dotCos>0.9995 ) {
868  //cout<<"Found Matched Single CRT Tag with CRT*TPC: "<<fabs(allUniqueTracksPair[u].dotProductCos)<<endl;
869  //cout<<"Displacement of match:"<<deltaX<<','<<deltaY<<endl;
870 
871  std::vector<float> hitF;
872  std::vector<float> hitB;
873  hitF.push_back(X_CRT); hitF.push_back(Y_CRT); hitF.push_back(Z_CRT);
874  hitB.push_back(trackX1); hitB.push_back(trackY1); hitB.push_back(trackZ1);
875  CRTTrack->push_back(anab::CosmicTag(hitF,hitB, fabs(allUniqueTracksPair[u].dotProductCos),anab::CosmicTagID_t::kUnknown));
876  if (Z_CRT<100) T0col->push_back(anab::T0(CRTT0, 1, 1,CRTTrackId,fabs(allUniqueTracksPair[u].dotProductCos)));
877  else T0col->push_back(anab::T0(CRTT0, 2, 1,CRTTrackId,fabs(allUniqueTracksPair[u].dotProductCos) ));
878  auto const crtTrackPtr = crtPtr(CRTTrack->size()-1);
879  auto const t0CP = t0CandPtr(CRTTrackId);
880  CRTT0assn->addSingle(crtTrackPtr,t0CP);
881 
882  util::CreateAssn(*this, event, *T0col, trackList[TPCTrackId], *TPCT0assn);
883  util::CreateAssn(*this, event, *CRTTrack, trackList[TPCTrackId], *TPCCRTassn);
884  util::CreateAssn(*this, event, *CRTTrack, crtList[allUniqueTracksPair[u].trigNumberX], *CRTTriggerassn);
885  util::CreateAssn(*this, event, *CRTTrack, crtList[allUniqueTracksPair[u].trigNumberY], *CRTTriggerassn);
886 
887  }
888  }
889  }
890  nEvents++;
891 event.put(std::move(T0col)); event.put(std::move(CRTTrack)); event.put(std::move(TPCCRTassn)); event.put(std::move(TPCT0assn)); event.put(std::move(CRTT0assn)); event.put(std::move(CRTTriggerassn));
892 
893 }
894 
895 
896 
897 
std::string string
Definition: nybbler.cc:12
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.
SingleCRTMatchingProducer(fhicl::ParameterSet const &p)
art framework interface to geometry description
bool isRealData() const
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
void produce(art::Event &event) override
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
Declaration of signal hit object.
bool operator()(const tracksPair &pair1, const tracksPair &pair2)
void End(void)
Definition: gXSecComp.cxx:210
def center(depos, point)
Definition: depos.py:117
Provides recob::Track data product.
bool moduleMatcher(int module1, int module2)
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.