HoughTrack_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: HoughTrack
3 // Module Type: analyzer
4 // File: HoughTrack_module.cc
5 //
6 // Generated at Thu May 5 08:55:45 2016 by Stefano Tognini using artmod
7 // from cetpkgsupport v1_08_07.
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include "art/Framework/Core/EDAnalyzer.h"
11 #include "art/Framework/Core/ModuleMacros.h"
12 #include <canvas/Persistency/Common/FindManyP.h>
13 
14 #include "art/Framework/Principal/Event.h"
15 #include "art/Framework/Principal/Handle.h"
16 #include "art/Framework/Principal/Run.h"
17 #include "art/Framework/Principal/SubRun.h"
18 
19 #include "art/Framework/Services/Optional/TFileService.h"
20 
21 #include "canvas/Utilities/InputTag.h"
22 #include "fhiclcpp/ParameterSet.h"
23 #include "messagefacility/MessageLogger/MessageLogger.h"
24 
25 #include "Geometry/Geometry.h"
26 #include "GeometryObjects/Geo.h"
27 #include "GeometryObjects/CellGeo.h"
28 
29 #include "RawData/RawTrigger.h"
30 #include "RawData/RawDigit.h"
31 #include "RawData/FlatDAQData.h"
32 
33 #include "DAQDataFormats/RawEvent.h"
34 #include "DAQDataFormats/RawTrigger.h"
35 #include "DAQDataFormats/RawTriggerMask.h"
36 #include "DAQDataFormats/RawTriggerHeader.h"
37 #include "DAQDataFormats/RawDataBlock.h"
38 #include "DAQDataFormats/RawMicroBlock.h"
39 #include "DAQDataFormats/RawMicroSlice.h"
40 #include "DAQDataFormats/RawNanoSlice.h"
41 
42 #include "NovaDAQUtilities/TimeUtils.h"
43 #include <NovaTimingUtilities/TimingUtilities.h>
44 
45 
46 #include "RecoBase/CellHit.h"
47 #include "RecoBase/Cluster.h"
48 #include "RecoBase/RecoHit.h"
49 #include "RecoBase/HoughResult.h"
50 #include "RecoBase/Track.h"
51 
52 #include "Utilities/func/MathUtil.h"
53 
54 #include "MCCheater/BackTracker.h"
55 
56 #include "Simulation/ParticleNavigator.h"
57 #include "Simulation/Particle.h"
58 #include "Simulation/FLSHit.h"
59 #include "Simulation/FLSHitList.h"
60 #include "nusimdata/SimulationBase/MCTruth.h"
61 
62 #include "Track2D.h"
63 #include "Track3D.h"
64 
65 #include <TH1.h>
66 #include <TTimeStamp.h>
67 #include <TTree.h>
68 #include <TMath.h>
69 #include "TTimeStamp.h"
70 
71 #include <iostream>
72 #include <vector>
73 
74 #ifdef __MAKECINT__
75 #pragma link C++ class vector<double>+;
76 #endif
77 
78 
79 
80 //........................................[ NAMESPACE ]........................................
81 
82 namespace htk {
83  struct Track2D;
84  struct Track_Match;
85  class HoughTrack;
86  bool sortByNCell(Track2D &track_i, Track2D &track_j);
87  bool sortByFLSHits(sim::Particle particle_i, sim::Particle particle_j);
88 
89 }
90 
91 
92 
93 
94 //........................................[ STRUCT ]........................................
95 
96 struct htk::Track_Match {
97 
98  Track_Match(std::vector<Track2D>::const_iterator x, std::vector<Track2D>::const_iterator y):
99  x_track(x), y_track(y) {}
100 
101  std::vector<Track2D>::const_iterator x_track, y_track;
102 
103 };
104 
105 
106 
107 //........................................[ CLASS ]........................................
108 
109 class htk::HoughTrack : public art::EDAnalyzer {
110 public:
111  explicit HoughTrack(fhicl::ParameterSet const & p);
112  // The destructor generated by the compiler is fine for classes
113  // without bare pointers or other resource use.
114 
115  // Plugins should not be copied or assigned.
116  HoughTrack(HoughTrack const &) = delete;
117  HoughTrack(HoughTrack &&) = delete;
118  HoughTrack & operator = (HoughTrack const &) = delete;
119  HoughTrack & operator = (HoughTrack &&) = delete;
120 
121  // Required functions.
122  void analyze(art::Event const & e) override;
123 
124 
125  // Selected optional functions.
126  void beginJob() override;
127  void endJob() override;
128 
129 private:
130 
131  std::string hough_label_;
132  std::string slice_label_;
133  unsigned int hough_peak_cut_;
134  unsigned track2D_min_hits_;
135  double max_distance_between_hits_;
136  double fake_tracks_threshold_;
137  unsigned track3D_min_hits_;
138  int hit_distance_cut_;
139  double vertex_max_distance_;
140 
141  unsigned event_counter_;
142 
143  bool hit_is_on_track
144  (art::Ptr<rb::CellHit> const& hit, geo::View_t view, double const& slope, double const& intercept) const;
145 
146  art::ServiceHandle<cheat::BackTracker> bt;
147  bool isMC;
148 
149  art::ServiceHandle<geo::Geometry> geometry_;
150 
151  unsigned int totalntracksX = 0;
152  unsigned int totalntracksY = 0;
153  unsigned int total3DTracks = 0;
154  unsigned int failed3Devents = 0;
155 
156 
157  TTree *tree_Events;
158  unsigned long run;
159  unsigned long subrun;
160  unsigned long event;
161  unsigned long long timestamp;
162  unsigned long long triggerStartUnixTime;
163  int triggerLength;
164  int year;
165  int month;
166  int day;
167  int doy;
168  int hour;
169  int minute;
170  int second;
171  long nanoSecond;
172  int UTCyear;
173  int UTCmonth;
174  int UTCday;
175  int UTCdoy;
176  int UTChour;
177  int UTCminute;
178  int UTCsecond;
179  long UTCnanoSecond;
180  int ntracksXZ;
181  int ntracksYZ;
182  int ntracks3D;
183 
184 
185  TTree *tree_Tracks;
186  std::vector<double> nhits;
187  std::vector<double> nplanesXZ;
188  std::vector<double> nplanesYZ;
189  std::vector<double> nplanes;
190  std::vector<double> startPlaneXZ;
191  std::vector<double> endPlaneXZ;
192  std::vector<double> startPlaneYZ;
193  std::vector<double> endPlaneYZ;
194  std::vector<double> len;
195  std::vector<double> startTime;
196  std::vector<double> meanTime;
197  std::vector<double> calE;
198  std::vector<double> startX;
199  std::vector<double> startY;
200  std::vector<double> startZ;
201  std::vector<double> endX;
202  std::vector<double> endY;
203  std::vector<double> endZ;
204  std::vector<double> dirX;
205  std::vector<double> dirY;
206  std::vector<double> dirZ;
207  std::vector<double> thetaXZ;
208  std::vector<double> thetaYZ;
209  std::vector<double> theta;
210  std::vector<double> phi;
211  std::vector<double> cosTheta;
212  std::vector<double> cosPhi;
213  std::vector<double> zenith;
214  std::vector<double> azimuth;
215  std::vector<double> cosZenith;
216  std::vector<double> cosAzimuth;
217 
218 
219 
220 
221  TTree *tree_Truth;
222  std::vector<double> particlePDG;
223  std::vector<double> particleMass;
224  std::vector<double> particleFLSHits;
225  std::vector<double> particleNplanesXZ;
226  std::vector<double> particleNplanesYZ;
227  std::vector<double> particlePx;
228  std::vector<double> particlePy;
229  std::vector<double> particlePz;
230  std::vector<double> particleE;
231  std::vector<double> particleStartX;
232  std::vector<double> particleStartY;
233  std::vector<double> particleStartZ;
234  std::vector<double> particleEndX;
235  std::vector<double> particleEndY;
236  std::vector<double> particleEndZ;
237  std::vector<double> particleDirX;
238  std::vector<double> particleDirY;
239  std::vector<double> particleDirZ;
240  std::vector<double> particleLen;
241  std::vector<double> particleT;
242  std::vector<double> particleTheta;
243  std::vector<double> particlePhi;
244  std::vector<double> particleCosTheta;
245  std::vector<double> particleCosPhi;
246 
247  TTree *tree_Rates;
248  int startYear;
249  int startMonth;
250  int startDay;
251  int startDoy;
252  int startHour;
253  int startMinute;
254  int startSecond;
255  int endYear;
256  int endMonth;
257  int endDay;
258  int endDoy;
259  int endHour;
260  int endMinute;
261  int endSecond;
262  int UTCstartYear;
263  int UTCstartMonth;
264  int UTCstartDay;
265  int UTCstartDoy;
266  int UTCstartHour;
267  int UTCstartMinute;
268  int UTCstartSecond;
269  int UTCendYear;
270  int UTCendMonth;
271  int UTCendDay;
272  int UTCendDoy;
273  int UTCendHour;
274  int UTCendMinute;
275  int UTCendSecond;
276  double exposureTime;
277  int nSinglemuEvents = 0;
278  int nMultimuEvents = 0;
279  double singlemuRate;
280  double multimuRate;
281 
282 
283  TTree *tree_Hough;
284  std::vector<double> houghView;
285  std::vector<double> houghRho;
286  std::vector<double> houghTheta;
287  std::vector<double> houghHeight;
288  std::vector<double> houghSlope;
289  std::vector<double> houghIntercept;
290 
291 
292  double totalTracks = 0;
293  double startTimeInSeconds;
294  double endTimeInSeconds;
295 
296 
297 };
298 
299 
300 
301 //........................................[ PARAMETER SET ]........................................
302 
303 htk::HoughTrack::HoughTrack(fhicl::ParameterSet const & p):
304 EDAnalyzer(p),
305 hough_label_ (p.get<std::string>("hough_label")),
306 slice_label_ (p.get<std::string>("slice_label")),
307 hough_peak_cut_ (p.get<unsigned int>("hough_peak_cut")),
308 track2D_min_hits_ (p.get<unsigned int>("track2D_min_hits")),
309 max_distance_between_hits_ (p.get<double>("max_distance_between_hits")),
310 fake_tracks_threshold_ (p.get<double>("fake_tracks_threshold")),
311 track3D_min_hits_ (p.get<unsigned>("track3D_min_hits")),
312 hit_distance_cut_ (p.get<int>("hit_distance_cut")),
313 vertex_max_distance_ (p.get<double>("vertex_max_distance")),
314 
315 event_counter_(0)
316 
317 {}
318 
319 
320 
321 //........................................[ BEGIN JOB ]........................................
322 
323 void htk::HoughTrack::beginJob() {
324 
325  art::ServiceHandle<art::TFileService> tfs;
326 
327  isMC = bt->HaveTruthInfo();
328 
329 
330  // Tree Events
331  tree_Events = tfs->make<TTree>("Events", "Events");
332 
333  tree_Events -> Branch("run", &run, "run/l");
334  tree_Events -> Branch("subrun", &subrun, "subrun/l");
335  tree_Events -> Branch("event", &event, "event/l");
336  tree_Events -> Branch("triggerStartUnixTime", &triggerStartUnixTime, "triggerStartUnixTime/l");
337  //tree_Events -> Branch("triggerLength", &triggerLength, "triggerLength/I");
338  tree_Events -> Branch("year", &year, "year/I");
339  tree_Events -> Branch("month", &month, "month/I");
340  tree_Events -> Branch("day", &day, "day/I");
341  tree_Events -> Branch("doy", &doy, "doy/I");
342  tree_Events -> Branch("hour", &hour, "hour/I");
343  tree_Events -> Branch("minute", &minute, "minute/I");
344  tree_Events -> Branch("second", &second, "second/I");
345  tree_Events -> Branch("nanoSecond", &nanoSecond, "nanoSecond/l");
346  tree_Events -> Branch("UTCyear", &UTCyear, "UTCyear/I");
347  tree_Events -> Branch("UTCmonth", &UTCmonth, "UTCmonth/I");
348  tree_Events -> Branch("UTCday", &UTCday, "UTCday/I");
349  tree_Events -> Branch("UTCdoy", &UTCdoy, "UTCdoy/I");
350  tree_Events -> Branch("UTChour", &UTChour, "UTChour/I");
351  tree_Events -> Branch("UTCminute", &UTCminute, "UTCminute/I");
352  tree_Events -> Branch("UTCsecond", &UTCsecond, "UTCsecond/I");
353  tree_Events -> Branch("UTCnanoSecond", &UTCnanoSecond, "UTCnanoSecond/l");
354  tree_Events -> Branch("ntracksXZ", &ntracksXZ, "ntracksXZ/I");
355  tree_Events -> Branch("ntracksYZ", &ntracksYZ, "ntracksYZ/I");
356  tree_Events -> Branch("ntracks3D", &ntracks3D, "ntracks3D/I");
357 
358  // Tree Tracks
359  tree_Tracks = tfs->make<TTree>("Tracks", "Tracks");
360  tree_Tracks -> Branch("nhits", &nhits);
361  tree_Tracks -> Branch("nplanesXZ", &nplanesXZ);
362  tree_Tracks -> Branch("nplanesYZ", &nplanesYZ);
363  tree_Tracks -> Branch("nplanes", &nplanes);
364  tree_Tracks -> Branch("startPlaneXZ", &startPlaneXZ);
365  tree_Tracks -> Branch("endPlaneXZ", &endPlaneXZ);
366  tree_Tracks -> Branch("startPlaneYZ", &startPlaneYZ);
367  tree_Tracks -> Branch("endPlaneYZ", &endPlaneYZ);
368  tree_Tracks -> Branch("len", &len);
369  tree_Tracks -> Branch("startTime", &startTime);
370  tree_Tracks -> Branch("meanTime", &meanTime);
371  tree_Tracks -> Branch("calE", &calE);
372  tree_Tracks -> Branch("startX", &startX);
373  tree_Tracks -> Branch("startY", &startY);
374  tree_Tracks -> Branch("startZ", &startZ);
375  tree_Tracks -> Branch("endX", &endX);
376  tree_Tracks -> Branch("endY", &endY);
377  tree_Tracks -> Branch("endZ", &endZ);
378  tree_Tracks -> Branch("dirX", &dirX);
379  tree_Tracks -> Branch("dirY", &dirY);
380  tree_Tracks -> Branch("dirZ", &dirZ);
381  tree_Tracks -> Branch("thetaXZ", &thetaXZ);
382  tree_Tracks -> Branch("thetaYZ", &thetaYZ);
383  tree_Tracks -> Branch("theta", &theta);
384  tree_Tracks -> Branch("phi", &phi);
385  tree_Tracks -> Branch("cosTheta", &cosTheta);
386  tree_Tracks -> Branch("cosPhi", &cosPhi);
387  tree_Tracks -> Branch("zenith", &zenith);
388  tree_Tracks -> Branch("azimuth", &azimuth);
389  tree_Tracks -> Branch("cosZenith", &cosZenith);
390  tree_Tracks -> Branch("cosAzimuth", &cosAzimuth);
391 
392  // Tree Truth (MC only)
393  tree_Truth = tfs->make<TTree>("Truth", "Truth");
394  tree_Truth -> Branch("particlePDG", &particlePDG);
395  tree_Truth -> Branch("particleMass", &particleMass);
396  tree_Truth -> Branch("particleFLSHits", &particleFLSHits);
397  tree_Truth -> Branch("particleNplanesXZ", &particleNplanesXZ);
398  tree_Truth -> Branch("particleNplanesYZ", &particleNplanesYZ);
399  tree_Truth -> Branch("particlePx", &particlePx);
400  tree_Truth -> Branch("particlePy", &particlePy);
401  tree_Truth -> Branch("particlePz", &particlePz);
402  tree_Truth -> Branch("particleE", &particleE);
403  tree_Truth -> Branch("particleStartX", &particleStartX);
404  tree_Truth -> Branch("particleStartY", &particleStartY);
405  tree_Truth -> Branch("particleStartZ", &particleStartZ);
406  tree_Truth -> Branch("particleEndX", &particleEndX);
407  tree_Truth -> Branch("particleEndY", &particleEndY);
408  tree_Truth -> Branch("particleEndZ", &particleEndZ);
409  tree_Truth -> Branch("particleDirX", &particleDirX);
410  tree_Truth -> Branch("particleDirY", &particleDirY);
411  tree_Truth -> Branch("particleDirZ", &particleDirZ);
412  tree_Truth -> Branch("particleLen", &particleLen);
413  tree_Truth -> Branch("particleT", &particleT);
414  tree_Truth -> Branch("particleTheta", &particleTheta);
415  tree_Truth -> Branch("particlePhi", &particlePhi);
416  tree_Truth -> Branch("particleCosTheta", &particleCosTheta);
417  tree_Truth -> Branch("particleCosPhi", &particleCosPhi);
418 
419  // Tree Rates
420  tree_Rates = tfs->make<TTree>("Rates", "Rates");
421  tree_Rates -> Branch("startYear", &startYear, "startYear/I");
422  tree_Rates -> Branch("startMonth", &startMonth, "startMonth/I");
423  tree_Rates -> Branch("startDay", &startDay, "startDay/I");
424  tree_Rates -> Branch("startDoy", &startDoy, "startDoy/I");
425  tree_Rates -> Branch("startHour", &startHour, "startHour/I");
426  tree_Rates -> Branch("startMinute", &startMinute, "startMinute/I");
427  tree_Rates -> Branch("startSecond", &startSecond, "startSecond/I");
428  tree_Rates -> Branch("endYear", &endYear, "endYear/I");
429  tree_Rates -> Branch("endMonth", &endMonth, "endMonth/I");
430  tree_Rates -> Branch("endDay", &endDay, "endDay/I");
431  tree_Rates -> Branch("endDoy", &endDoy, "endDoy/I");
432  tree_Rates -> Branch("endHour", &endHour, "endHour/I");
433  tree_Rates -> Branch("endMinute", &endMinute, "endMinute/I");
434  tree_Rates -> Branch("endSecond", &endSecond, "endSecond/I");
435  tree_Rates -> Branch("UTCstartYear", &UTCstartYear, "UTCstartYear/I");
436  tree_Rates -> Branch("UTCstartMonth", &UTCstartMonth, "UTCstartMonth/I");
437  tree_Rates -> Branch("UTCstartDay", &UTCstartDay, "UTCstartDay/I");
438  tree_Rates -> Branch("UTCstartDoy", &UTCstartDoy, "UTCstartDoy/I");
439  tree_Rates -> Branch("UTCstartHour", &UTCstartHour, "UTCstartHour/I");
440  tree_Rates -> Branch("UTCstartMinute", &UTCstartMinute, "UTCstartMinute/I");
441  tree_Rates -> Branch("UTCstartSecond", &UTCstartSecond, "UTCstartSecond/I");
442  tree_Rates -> Branch("UTCendYear", &UTCendYear, "UTCendYear/I");
443  tree_Rates -> Branch("UTCendMonth", &UTCendMonth, "UTCendMonth/I");
444  tree_Rates -> Branch("UTCendDay", &UTCendDay, "UTCendDay/I");
445  tree_Rates -> Branch("UTCendDoy", &UTCendDoy, "UTCendDoy/I");
446  tree_Rates -> Branch("UTCendHour", &UTCendHour, "UTCendHour/I");
447  tree_Rates -> Branch("UTCendMinute", &UTCendMinute, "UTCendMinute/I");
448  tree_Rates -> Branch("UTCendSecond", &UTCendSecond, "UTCendSecond/I");
449  tree_Rates -> Branch("exposureTime", &exposureTime, "exposureTime/D");
450  tree_Rates -> Branch("nSinglemuEvents", &nSinglemuEvents, "nSinglemuEvents/I");
451  tree_Rates -> Branch("nMultimuEvents", &nMultimuEvents, "nMultimuEvents/I");
452  tree_Rates -> Branch("singlemuRate", &singlemuRate, "singlemuRate/D");
453  tree_Rates -> Branch("multimuRate", &multimuRate, "multimuRate/D");
454 
455  // Tree HoughResults
456  tree_Hough = tfs->make<TTree>("HoughResults", "HoughResults");
457  tree_Hough -> Branch("view", &houghView);
458  tree_Hough -> Branch("rho", &houghRho);
459  tree_Hough -> Branch("theta", &houghTheta);
460  tree_Hough -> Branch("slope", &houghSlope);
461  tree_Hough -> Branch("intercept", &houghIntercept);
462  tree_Hough -> Branch("height", &houghHeight);
463 
464 
465 }
466 
467 
468 
469 //........................................[ ANALYZER ]........................................
470 
471 void htk::HoughTrack::analyze(art::Event const &e) {
472 
473 
474 
475 
476  //##########################################################
477  //
478  // TRUE INFORMATION FROM MC FILES
479  //
480  //##########################################################
481 
482 
483  isMC = bt->HaveTruthInfo();
484 
485  if (!isMC) {
486 
487  if (event_counter_ == 0) {
488 
489  particlePDG.push_back(99);
490  particleMass.push_back(99);
491  particleFLSHits.push_back(99);
492  particleNplanesXZ.push_back(99);
493  particleNplanesYZ.push_back(99);
494  particlePx.push_back(99);
495  particlePy.push_back(99);
496  particlePz.push_back(99);
497  particleE.push_back(99);
498  particleStartX.push_back(99);
499  particleStartY.push_back(99);
500  particleStartZ.push_back(99);
501  particleEndX.push_back(99);
502  particleEndY.push_back(99);
503  particleEndZ.push_back(99);
504  particleDirX.push_back(99);
505  particleDirY.push_back(99);
506  particleDirZ.push_back(99);
507  particleLen.push_back(99);
508  particleT.push_back(99);
509  particleTheta.push_back(99);
510  particlePhi.push_back(99);
511  particleCosTheta.push_back(99);
512  particleCosPhi.push_back(99);
513 
514  }
515 
516  }
517 
518  if (isMC) {
519 
520 
521  // Starting a fresh vector
522  particlePDG.clear();
523  particleMass.clear();
524  particleFLSHits.clear();
525  particleNplanesXZ.clear();
526  particleNplanesYZ.clear();
527  particlePx.clear();
528  particlePy.clear();
529  particlePz.clear();
530  particleE.clear();
531  particleStartX.clear();
532  particleStartY.clear();
533  particleStartZ.clear();
534  particleEndX.clear();
535  particleEndY.clear();
536  particleEndZ.clear();
537  particleDirX.clear();
538  particleDirY.clear();
539  particleDirZ.clear();
540  particleLen.clear();
541  particleT.clear();
542  particleTheta.clear();
543  particlePhi.clear();
544  particleCosTheta.clear();
545  particleCosPhi.clear();
546 
547 
548  const sim::ParticleNavigator &pnav = bt->ParticleNavigator();
549 
550 
551  // Looping over particles in the event
552 
553  for(sim::ParticleNavigator::const_iterator it = pnav.begin(); it != pnav.end(); ++it) {
554 
555  const sim::Particle* particle = (*it).second;
556 
557  const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(particle->TrackId());
558  const int nflshits = flshits.size();
559 
560  // Skipping particles with no FLSHits
561  if (nflshits == 0) { continue; }
562 
563  // Skipping particles that are not muons
564  if (TMath::Abs(particle->PdgCode()) != 13) { continue; }
565 
566  // Getting particles' true information
567  const TLorentzVector& fourPosition = particle->Position();
568  const TLorentzVector& fourMomentum = particle->Momentum();
569 
570 
571  // Create matrix (vector of vectors) of number of hits in each cell
572  std::vector<std::vector<double> > hitsByPlaneCell(geometry_->GetPlanesByView().size());
573  unsigned int numCells = geometry_->Plane(0)->Ncells();
574  if(geometry_->Plane(1)->Ncells() > numCells) { numCells = geometry_->Plane(1)->Ncells(); }
575 
576  // Resize each vector in the larger plane vector to have the number of cells, initialized to 0
577  for(unsigned int i_plane = 0, n_plane = hitsByPlaneCell.size(); i_plane < n_plane; ++i_plane) {
578  hitsByPlaneCell[i_plane].resize(numCells, 0.);
579  }
580 
581 
582  // Vector matching a plane to a view
583  std::vector<geo::View_t> planeView(geometry_->GetPlanesByView().size());
584 
585  double startPlane;
586  double startCell;
587  double endPlane;
588  double endCell;
589 
590  double entryPos[3] = {-9999, -9999, -9999};
591  double exitPos[3] = {9999, 9999, 9999};
592 
593  double localStartPos[3] = {-9999, -9999, -9999};
594  double localEndPos[3] = {9999, 9999, 9999};
595 
596  double entryPosTemp[3];
597  double exitPosTemp[3];
598 
599  // double pathLen = 0;
600 
601  double planesX = 0.;
602  //double cellsX = 0.;
603  double planesY = 0.;
604  //double cellsY = 0.;
605  //double hitCells = 0.;
606 
607  int FLShitsCounter = 0;
608 
609  // Loop over FLSHits for this particle
610  for (const auto& hit : flshits) {
611 
612  startPlane = hit.GetPlaneID();
613  startCell = hit.GetCellID();
614  localStartPos[0] = hit.GetEntryX();
615  localStartPos[1] = hit.GetEntryY();
616  localStartPos[2] = hit.GetEntryZ();
617 
618  endPlane = hit.GetPlaneID();
619  endCell = hit.GetCellID();
620  localEndPos[0] = hit.GetExitX();
621  localEndPos[1] = hit.GetExitY();
622  localEndPos[2] = hit.GetExitZ();
623 
624  //pathLen = pathLen + hit.GetTotalPathLength();
625 
626  hitsByPlaneCell[hit.GetPlaneID()][hit.GetCellID()] += 1.; // Add a hit to the proper cell
627  planeView[hit.GetPlaneID()] = geometry_->Plane(hit.GetPlaneID())->View();
628 
629  // Convert entry/exit locations to world coordinates
630  geometry_->Plane(startPlane)->Cell(startCell)->LocalToWorld(localStartPos, entryPosTemp);
631  geometry_->Plane(endPlane) ->Cell(endCell) ->LocalToWorld(localEndPos, exitPosTemp);
632 
633 
634  // For particleStart
635  if (entryPos[1] < entryPosTemp[1] || entryPos[1] < exitPosTemp[1]) {
636 
637  entryPos[0] = entryPosTemp[0];
638  entryPos[1] = entryPosTemp[1];
639  entryPos[2] = entryPosTemp[2];
640 
641  }
642 
643  // For particleEnd
644  if (exitPos[1] > exitPosTemp[1] || exitPos[1] > entryPosTemp[1] ) {
645 
646  exitPos[0] = entryPosTemp[0];
647  exitPos[1] = entryPosTemp[1];
648  exitPos[2] = entryPosTemp[2];
649 
650  }
651 
652  if (hit.GetPlaneID() % 2) { planesX++; }
653  else { planesY++; }
654 
655  FLShitsCounter++;
656 
657  } // end of loop over FLSHits
658 
659 
660  double pRx = exitPos[0] - entryPos[0];
661  double pRy = exitPos[1] - entryPos[1];
662  double pRz = exitPos[2] - entryPos[2];
663  double abspR = TMath::Sqrt( pRx*pRx + pRy*pRy + pRz*pRz );
664 
665  // Calculating the length
666  double deltaXsquare = (exitPos[0] - entryPos[0])*(exitPos[0] - entryPos[0]);
667  double deltaYsquare = (exitPos[1] - entryPos[1])*(exitPos[1] - entryPos[1]);
668  double deltaZsquare = (exitPos[2] - entryPos[2])*(exitPos[2] - entryPos[2]);
669  double particleLength = TMath::Sqrt( deltaXsquare + deltaYsquare + deltaZsquare );
670 
671  // Calculating zenith and azimuth
672  double pTheta = TMath::ACos( pRy / abspR );
673  double pZenith = TMath::Pi() - pTheta;
674  double pPhi = TMath::ATan(pRz / pRx); // Starts at the X+ axis
675  double pAzimuth = 0;
676 
677  if (pRx >= 0 && pRz >= 0) { pAzimuth = pPhi; }
678  if (pRx < 0 && pRz >= 0) { pAzimuth = TMath::Pi() + pPhi; }
679  if (pRx < 0 && pRz < 0) { pAzimuth = TMath::Pi() + pPhi; }
680  if (pRx >= 0 && pRz < 0) { pAzimuth = 2*TMath::Pi() + pPhi; }
681 
682  /*
683  // Loop over planes, cells
684  for (unsigned int i_plane = 0, n_plane = hitsByPlaneCell.size(); i_plane < n_plane; ++i_plane) {
685 
686  // Count the number of hit cells in this plane
687 
688  for (unsigned int i_cell = 0, n_cell = hitsByPlaneCell[i_plane].size(); i_cell < n_cell; ++i_cell) {
689 
690  if(hitsByPlaneCell[i_plane][i_cell] > 0.) { hitCells += 1.; }
691  }
692 
693  if (planeView[i_plane] == geo::kX && hitCells > 0.) {
694  planesX += 1.;
695  cellsX += hitCells;
696  }
697  if (planeView[i_plane] == geo::kY && hitCells > 0.) {
698  planesY += 1.;
699  cellsY += hitCells;
700  }
701  } // end loop over hit planes/cells
702  */
703 
704  // Filling the tree
705  particlePDG.push_back( particle->PdgCode() );
706  particleMass.push_back( particle->Mass() );
707  particleFLSHits.push_back( nflshits );
708  particleNplanesXZ.push_back(planesX);
709  particleNplanesYZ.push_back(planesY);
710  particlePx.push_back( fourMomentum.Px() );
711  particlePy.push_back( fourMomentum.Py() );
712  particlePz.push_back( fourMomentum.Pz() );
713  particleE.push_back( fourMomentum.E() );
714  particleStartX.push_back( entryPos[0] );
715  particleStartY.push_back( entryPos[1] );
716  particleStartZ.push_back( entryPos[2] );
717  particleEndX.push_back( exitPos[0] );
718  particleEndY.push_back( exitPos[1] );
719  particleEndZ.push_back( exitPos[2] );
720  particleDirX.push_back(pRx/abspR);
721  particleDirY.push_back(pRy/abspR);
722  particleDirZ.push_back(pRz/abspR);
723  particleLen.push_back( particleLength );
724  particleT.push_back( fourPosition.T() );
725  particleTheta.push_back( pZenith * TMath::RadToDeg() );
726  particlePhi.push_back( pAzimuth * TMath::RadToDeg() );
727  particleCosTheta.push_back( TMath::Cos(pZenith) );
728  particleCosPhi.push_back( TMath::Cos(pAzimuth) );
729 
730  //std::cout << "PARTICLE LEN: " << pathLen << std::endl;
731 
732 
733  }
734 
735 
736 
737  // If there are no muons in the event
738  if (particlePDG.size() == 0) {
739 
740  particlePDG.push_back(0);
741  particleMass.push_back(0);
742  particleFLSHits.push_back(0);
743  particleNplanesXZ.push_back(0);
744  particleNplanesYZ.push_back(0);
745  particlePx.push_back(0);
746  particlePy.push_back(0);
747  particlePz.push_back(0);
748  particleE.push_back(0);
749  particleStartX.push_back(0);
750  particleStartY.push_back(0);
751  particleStartZ.push_back(0);
752  particleEndX.push_back(0);
753  particleEndY.push_back(0);
754  particleEndZ.push_back(0);
755  particleDirX.push_back(0);
756  particleDirY.push_back(0);
757  particleDirZ.push_back(0);
758  particleLen.push_back(0);
759  particleT.push_back(0);
760  particleTheta.push_back(0);
761  particlePhi.push_back(0);
762  particleCosTheta.push_back(0);
763  particleCosPhi.push_back(0);
764 
765  }
766 
767  // Filling the tree
768  tree_Truth->Fill();
769 
770 
771 
772  }
773 
774 
775 
776 
777 
778  //############################################################
779  //
780  // RUN, SUBRUN, TIME, UTC TIME OF THE EVENT
781  //
782  //############################################################
783 
784  // For printing out a specific event
785  unsigned int eventToCheck = 0;
786 
787  // Extract Objects from the Event
788  art::Handle<std::vector<rb::Cluster> > my_slices;
789  e.getByLabel(slice_label_, my_slices);
790  if(my_slices.failedToGet()){
791  std::cout << "No slices in event " << e.event() << ", so doing nothing\n";
792  return;
793  }
794  art::FindManyP<rb::HoughResult> find_hough_results(my_slices, e, hough_label_);
795 
796 
797  // Local time
798  art::Timestamp ts = e.time();
799  const time_t timeSec = ts.timeHigh();
800  struct tm* timeStruct = localtime(&timeSec);
801 
802  // UTC time (from /nova/app/home/novasoft/slf6/novasoft/releases/development/EventDisplay/HeaderDrawer.cxx)
803  unsigned int UTCyear_, UTCmonth_, UTCday_, UTChour_, UTCminute_, UTCsecond_;
804 
805  unsigned long long int tsval = e.time().value();
806  const unsigned long int mask32 = 0xFFFFFFFFUL;
807  unsigned long int lup = ( tsval >> 32 ) & mask32;
808  unsigned long int llo = tsval & mask32;
809  TTimeStamp UTCts(lup, (int)llo);
810 
811  UTCts.GetDate(kTRUE,0,&UTCyear_, &UTCmonth_, &UTCday_);
812  UTCts.GetTime(kTRUE,0,&UTChour_, &UTCminute_, &UTCsecond_);
813 
814 
815  // Get the trigger information from event
816  art::Handle< std::vector<rawdata::RawTrigger> > rawTriggers;
817  e.getByLabel("daq", rawTriggers);
818  const rawdata::RawTrigger trig = rawTriggers->at(0);
819 
820  triggerStartUnixTime = trig.fTriggerTimingMarker_TimeStart;
821 
822  //triggerLength = trig.fTriggerRange_TriggerLength;
823 
824  struct timespec trigger_ts;
825  novadaq::timeutils::convertNovaTimeToUnixTime(triggerStartUnixTime, trigger_ts);
826 
827 
828 
829 
830 
831  ////////////////////////////////////////////
832 /*
833 
834 
835 
836 
837 
838 
839  static const int TDC_PER_US = 64;
840  static const int US_PER_MICROSLICE = 50; // I hope this is always true
841 
842  int64_t event_length_tdc = 0;
843  int64_t delta_tdc = 0;
844  art::Handle< std::vector<rawdata::FlatDAQData> > flatdaq;
845  art::Handle< std::vector<rawdata::RawTrigger> > rawtrigger;
846 
847  e.getByLabel("daq", rawtrigger);
848  e.getByLabel("daq", flatdaq);
849 
850 
851  const rawdata::FlatDAQData thisflatdaq = flatdaq->at(0);
852  const rawdata::RawTrigger thisrawtrigger = rawtrigger->at(0);
853 
854  // std::cerr << flatdaq->size() << std::endl;
855 
856 
857  // daqdataformats::RawEvent *raw;
858 
859 
860 // raw->readData(thisflatdaq.getRawBufferPointer());
861 
862 
863  //raw->setFloatingDataBlock(0);
864 
865 
866  // daqdataformats::RawDataBlock datablock = raw->getFloatingDataBlock();
867 
868  uint64_t event_start_time = 0xffffffffffffffff;
869  uint64_t event_end_time = 0x0000000000000000;
870 
871 
872  for(unsigned int di = 0; di < raw->getDataBlockNumber(); di++){
873  raw->setFloatingDataBlock(di);
874  datablock = (raw->getFloatingDataBlock());
875 
876 
877  if(datablock.getHeader()->getMarker() ==
878  daqdataformats::datablockheader::SummaryBlock_Marker ||
879  !datablock.getHeader()->checkMarker()) continue;
880 
881 
882  for(unsigned int mi = 0; mi < datablock.getNumMicroBlocks(); mi++){
883  datablock.setFloatingMicroBlock(mi);
884  daqdataformats::RawMicroBlock * ub = datablock.getFloatingMicroBlock();
885 
886  // The time is always in the second and third words of the
887  // microslice, which follows two words of microblock header, so
888  // just get it. Justin says you can also get it from getTime(),
889  // but this already works and I'm not touching it.
890  const uint32_t t_marker_low = ((uint32_t *)(ub->getBuffer()))[3];
891  const uint32_t t_marker_high = ((uint32_t *)(ub->getBuffer()))[4];
892 
893  uint64_t time_marker = t_marker_low;
894  time_marker |= (uint64_t)t_marker_high << 32;
895  if(time_marker < event_start_time) event_start_time = time_marker;
896  if(time_marker > event_end_time ) event_end_time = time_marker;
897  }
898  }
899 
900  delta_tdc = (int64_t)((thisrawtrigger.fTriggerTimingMarker_TimeStart - event_start_time));
901 
902  // Assume that microblocks are always 50us. I hope that's true for all
903  // relevant data.
904  event_length_tdc = ((int64_t)(event_end_time - event_start_time)) + US_PER_MICROSLICE*TDC_PER_US;
905  */
906 
907 
908 
909  ////////////////////////////////////////////
910 
911 
912 
913 
914 
915 
916 
917 
918 
919 
920  // Filling the Events Tree
921  run = e.run();
922  subrun = e.subRun();
923  event = e.event();
924 
925  year = timeStruct->tm_year + 1900;
926  month = timeStruct->tm_mon + 1;
927  day = timeStruct->tm_mday;
928  doy = timeStruct->tm_yday + 1;
929  hour = timeStruct->tm_hour + 1;
930  minute = timeStruct->tm_min;
931  second = timeStruct->tm_sec;
932  nanoSecond = trigger_ts.tv_nsec;
933 
934  UTCyear = UTCyear_;
935  UTCmonth = UTCmonth_;
936  UTCday = UTCday_;
937  UTCdoy = UTCts.GetDayOfYear();
938  UTChour = UTChour_;
939  UTCminute = UTCminute_;
940  UTCsecond = UTCsecond_;
941  UTCnanoSecond = trigger_ts.tv_nsec;
942 
943 
944 
945  //std::cerr << e.event() << ": " << UTChour << ":" << UTCminute << ":" << UTCsecond << "." << UTCnanoSecond << " | " << event_start_time << " " << event_length_tdc << std::endl;
946 
947 
948 
949 
950  //############################################################
951  //
952  // HOUGH TRACK RECONSTRUCTION FOR BOTH DATA AND MC FILES
953  //
954  //############################################################
955 
956  //starting a fresh vector (tracks tree)
957  nhits.clear();
958  nplanesXZ.clear();
959  nplanesYZ.clear();
960  nplanes.clear();
961  startPlaneXZ.clear();
962  endPlaneXZ.clear();
963  startPlaneYZ.clear();
964  endPlaneYZ.clear();
965  len.clear();
966  startTime.clear();
967  meanTime.clear();
968  calE.clear();
969 
970  startX.clear();
971  startY.clear();
972  startZ.clear();
973 
974  endX.clear();
975  endY.clear();
976  endZ.clear();
977 
978  dirX.clear();
979  dirY.clear();
980  dirZ.clear();
981 
982  thetaXZ.clear();
983  thetaYZ.clear();
984 
985  theta.clear();
986  phi.clear();
987  cosTheta.clear();
988  cosPhi.clear();
989 
990  zenith.clear();
991  azimuth.clear();
992  cosZenith.clear();
993  cosAzimuth.clear();
994 
995  //starting a fresh vector (hough tree)
996  houghView.clear();
997  houghRho.clear();
998  houghTheta.clear();
999  houghSlope.clear();
1000  houghIntercept.clear();
1001  houghHeight.clear();
1002 
1003 
1004  ntracksXZ = 0;
1005  ntracksYZ = 0;
1006 
1007 
1008  if (event == eventToCheck){
1009  std::cerr << "EVENT " << event << std::endl;
1010  std::cerr << "------------" << std::endl;
1011  std::cerr << std::endl;
1012  }
1013 
1014 
1015 
1016  //....................[ Creating 2D tracks ]....................
1017 
1018  std::map<geo::View_t, std::vector<Track2D> > hough_tracks;
1019 
1020  // Looping over slices
1021  for (size_t n_slice = 0; n_slice != my_slices->size(); ++n_slice)
1022  {
1023  std::vector<art::Ptr<rb::HoughResult> > hough_results;
1024  find_hough_results.get(n_slice, hough_results);
1025 
1026 
1027 
1028  for (auto const& hough_result : hough_results)
1029  {
1030 
1031  unsigned peakInSlice = 0;
1032 
1033  // Looping over Hough peaks
1034  for (auto const& peak : hough_result->fPeak)
1035  {
1036 
1037  // Hough peak cut
1038  if (peak.fH > hough_peak_cut_)
1039  {
1040  geo::View_t view = hough_result->fView;
1041 
1042  double slope = -9999;
1043  double intercept = -9999;
1044  hough_result->SlopeIcept(peakInSlice++, &slope, &intercept);
1045 
1046  houghView.push_back(view);
1047  houghRho.push_back(peak.fRho);
1048  houghTheta.push_back(peak.fTheta);
1049  houghSlope.push_back(slope);
1050  houghIntercept.push_back(intercept);
1051  houghHeight.push_back(peak.fH);
1052 
1053  Track2D hough_track(view, slope, intercept);
1054 
1055  if (event == eventToCheck) {
1056 
1057  std::cout << "2D: view, slope, theta: " << view << ", " << hough_track.slope << ", " << TMath::Abs(TMath::ATan(hough_track.slope)) << std::endl;
1058 
1059  }
1060 
1061  // Sort hits by plane
1062  art::PtrVector<rb::CellHit> hits_planeSorted = my_slices->at(n_slice).AllCells();
1063  SortByPlane(hits_planeSorted);
1064 
1065 
1066  // Looping over hits
1067  for (auto const& hit : hits_planeSorted)
1068  {
1069 
1070  if (hit->View() != view) continue;
1071 
1072  // Calculating the DOCA for every hit
1073  // doca = |v x w| / |v|
1074  // v is a vector along the Hough track and w a vector connecting the Hough track to the hit
1075  double hit_XY = 0.;
1076  double hit_Z = 0.;
1077 
1078  double hit_xyz[3];
1079  geometry_->CellInfo(hit->Plane(), hit->Cell(), &view, hit_xyz);
1080  TVector3 hit_vector(hit_xyz);
1081 
1082  if (view == geo::View_t::kX) hit_XY = hit_vector.x();
1083  if (view == geo::View_t::kY) hit_XY = hit_vector.y();
1084 
1085  hit_Z = hit_vector.z();
1086 
1087  double wX = hit_Z - ( (hit_XY - intercept) / slope );
1088  double doca = TMath::Abs( wX * TMath::Sin(TMath::ATan(slope)) );
1089 
1090  double houghTheta = TMath::Abs(TMath::ATan(hough_track.slope));
1091  double sinTheta = TMath::Sin(houghTheta);
1092  double max_hit_distance = hit_distance_cut_ + hit_distance_cut_ * sinTheta;
1093 
1094 
1095  // Adding hit to the cluster
1096  if (doca <= max_hit_distance) {
1097 
1098  hough_track.cluster.Add(hit);
1099 
1100  /*
1101  if (event == eventToCheck) {
1102  if (view == geo::View_t::kX) std::cerr << "[XZ ]" << std::endl;
1103  if (view == geo::View_t::kY) std::cerr << "[ YZ]" << std::endl;
1104  std::cerr << "HIT: " << hit_XY << ", " << hit_Z << "; " << hit->ADC() << std::endl;
1105  std::cerr << "TRK: " << slope * hit_vector.z() + intercept << ", " << hit_Z << std::endl;
1106  std::cerr << "DCA: " << doca << " < " << max_hit_distance << std::endl;
1107  std::cerr << std::endl;
1108  }
1109  */
1110 
1111  }
1112  }
1113 
1114  // 2D Hough track criteria
1115  if (hough_track.cluster.NCell() > track2D_min_hits_) {
1116 
1117  hough_tracks[view].push_back(hough_track);
1118 
1119  if (view == geo::View_t::kX) ntracksXZ++;
1120  if (view == geo::View_t::kY) ntracksYZ++;
1121 
1122  }
1123  }
1124  }
1125  }
1126  }
1127 
1128 
1129  tree_Hough->Fill();
1130 
1131 
1132  //....................[ If no 3D tracks are possible ]....................
1133 
1134  if (ntracksXZ == 0 || ntracksYZ == 0) {
1135 
1136  ntracks3D = 0;
1137  //std::cerr << "++ Failed 3D reconstruction in Event " << event << std::endl;
1138 
1139  // Tracks Tree
1140  nhits.push_back(0);
1141  nplanesXZ.push_back(0);
1142  nplanesYZ.push_back(0);
1143  nplanes.push_back(0);
1144  startPlaneXZ.push_back(0);
1145  endPlaneXZ.push_back(0);
1146  startPlaneYZ.push_back(0);
1147  endPlaneYZ.push_back(0);
1148  len.push_back(0);
1149  startTime.push_back(0);
1150  meanTime.push_back(0);
1151  calE.push_back(0);
1152 
1153  startX.push_back(0);
1154  startY.push_back(0);
1155  startZ.push_back(0);
1156  endX.push_back(0);
1157  endY.push_back(0);
1158  endZ.push_back(0);
1159 
1160  dirX.push_back(0);
1161  dirY.push_back(0);
1162  dirZ.push_back(0);
1163 
1164  thetaXZ.push_back(0);
1165  thetaYZ.push_back(0);
1166 
1167  theta.push_back(0);
1168  phi.push_back(0);
1169  cosTheta.push_back(0);
1170  cosPhi.push_back(0);
1171 
1172  zenith.push_back(0);
1173  azimuth.push_back(0);
1174  cosZenith.push_back(0);
1175  cosAzimuth.push_back(0);
1176 
1177  tree_Tracks->Fill();
1178 
1179  // Events Tree
1180  totalntracksX = totalntracksX + ntracksXZ;
1181  totalntracksY = totalntracksY + ntracksYZ;
1182 
1183  // Getting the elapsed time (for the rate)
1184  if (event_counter_ == 0) {
1185 
1186  startYear = year;
1187  startMonth = month;
1188  startDay = day;
1189  startDoy = doy;
1190  startHour = hour;
1191  startMinute = minute;
1192  startSecond = second;
1193 
1194  UTCstartYear = UTCyear;
1195  UTCstartMonth = UTCmonth;
1196  UTCstartDay = UTCday;
1197  UTCstartDoy = UTCdoy;
1198  UTCstartHour = UTChour;
1199  UTCstartMinute = UTCminute;
1200  UTCstartSecond = UTCsecond;
1201 
1202  startTimeInSeconds = (startHour * 3600) + (startMinute * 60) + startSecond;
1203  }
1204 
1205  endYear = year;
1206  endMonth = month;
1207  endDay = day;
1208  endDoy = doy;
1209  endHour = hour;
1210  endMinute = minute;
1211  endSecond = second;
1212 
1213  UTCendYear = UTCyear;
1214  UTCendMonth = UTCmonth;
1215  UTCendDay = UTCday;
1216  UTCendDoy = UTCdoy;
1217  UTCendHour = UTChour;
1218  UTCendMinute = UTCminute;
1219  UTCendSecond = UTCsecond;
1220 
1221  tree_Events->Fill();
1222 
1223  event_counter_++;
1224  failed3Devents++;
1225 
1226  return;
1227 
1228  }
1229  //........................................................................
1230 
1231 
1232 
1233 
1234  // Sort Hough Tracks by NCell
1235  std::sort(hough_tracks[geo::kX].begin(), hough_tracks[geo::kX].end(), htk::sortByNCell);
1236  std::sort(hough_tracks[geo::kY].begin(), hough_tracks[geo::kY].end(), htk::sortByNCell);
1237 
1238  // Order Hough Tracks by descending order of NCell
1239  std::reverse(hough_tracks[geo::kX].begin(), hough_tracks[geo::kX].end());
1240  std::reverse(hough_tracks[geo::kY].begin(), hough_tracks[geo::kY].end());
1241 
1242 
1243 
1244  //....................[ Removing far way hits from 2D tracks ]....................
1245 
1246  // We loop over every 2D track created and try to remove hits that are far appart from the real track.
1247  // To do this, for every track, we start in the middle of its hits vector.
1248  // It is very likely to fall in the bulk of the track's hits and, therefore,
1249  // we will be very likely to be within the real track.
1250  // Then, we move from the center of the vector toward both ends of the track,
1251  // adding every hit that has a close distance to the next neighbor.
1252  // If a big gap is found, it is most likely that the track has ended.
1253  // Every hit after that point must be non related to the track and is ignored.
1254 
1255 
1256 
1257  std::vector<Track2D> x_clusteredTracks;
1258  std::vector<Track2D> y_clusteredTracks;
1259 
1260  // Looping over XZ tracks
1261  for (auto x_track = hough_tracks[geo::kX].begin(); x_track != hough_tracks[geo::kX].end(); x_track++) {
1262 
1263  // Sorting hits in by increasing value of Z
1264  art::PtrVector<rb::CellHit> hits_planeSorted = x_track->cluster.AllCells();
1265  SortByPlane(hits_planeSorted);
1266 
1267  int trackSize = x_track->cluster.NCell();
1268  int halfTrackSize = trackSize/2;
1269 
1270  double x_slope = x_track->slope;
1271  double x_intercept = x_track->intercept;
1272 
1273  double x_theta = TMath::Abs(TMath::ATan(x_track->slope));
1274  double x_sinTheta = TMath::Sin(x_theta);
1275  double maxGap = max_distance_between_hits_ + max_distance_between_hits_ * x_sinTheta;
1276 
1277 
1278  if (event == eventToCheck) {
1279  std::cout << "XZ: NCell, theta, halftrack, maxGap: " << trackSize << ", " << x_theta * TMath::RadToDeg() << ", " << halfTrackSize << ", " << maxGap << std::endl;
1280  }
1281 
1282  Track2D x_clusteredTrack(geo::kX, x_slope, x_intercept);
1283 
1284  // From Z = halfTrackSize to Zmin
1285  for (int hit_i = halfTrackSize; hit_i > 1; hit_i--) {
1286 
1287  art::Ptr<rb::CellHit> currentHit = x_track->cluster.Cell(hit_i);
1288  art::Ptr<rb::CellHit> previousHit = x_track->cluster.Cell(hit_i - 1);
1289 
1290  double currentHit_xyz[3];
1291  geometry_->CellInfo(currentHit->Plane(), currentHit->Cell(), NULL, currentHit_xyz);
1292  TVector3 currentHit_vector(currentHit_xyz);
1293 
1294  //if (x_theta == 0 && event == eventToCheck) std::cout << "hitXZ: " << currentHit_vector.x() << ", " << currentHit_vector.z() << std::endl;
1295 
1296  double currentHit_X = currentHit_vector.x();
1297  double currentHit_Z = currentHit_vector.z();
1298 
1299  double previousHit_xyz[3];
1300  geometry_->CellInfo(previousHit->Plane(), previousHit->Cell(), NULL, previousHit_xyz);
1301  TVector3 previousHit_vector(previousHit_xyz);
1302 
1303  double previousHit_X = previousHit_vector.x();
1304  double previousHit_Z = previousHit_vector.z();
1305 
1306  double deltaXsquare = (previousHit_X - currentHit_X)*(previousHit_X - currentHit_X);
1307  double deltaYsquare = (previousHit_Z - currentHit_Z)*(previousHit_Z - currentHit_Z);
1308 
1309  double distance_between_hits = sqrt(deltaXsquare + deltaYsquare);
1310 
1311  x_clusteredTrack.cluster.Add(currentHit);
1312 
1313  // If the track's continuity is broken, any hit afterwards is ignored
1314  if (distance_between_hits > maxGap) { break; }
1315 
1316 
1317  }
1318 
1319  // From Z = halfTrackSize to Zmax
1320  for (int hit_i = halfTrackSize; hit_i < trackSize - 1; hit_i++) {
1321 
1322  art::Ptr<rb::CellHit> currentHit = x_track->cluster.Cell(hit_i);
1323  art::Ptr<rb::CellHit> nextHit = x_track->cluster.Cell(hit_i + 1);
1324 
1325  double currentHit_xyz[3];
1326  geometry_->CellInfo(currentHit->Plane(), currentHit->Cell(), NULL, currentHit_xyz);
1327  TVector3 currentHit_vector(currentHit_xyz);
1328 
1329  double currentHit_X = currentHit_vector.x();
1330  double currentHit_Z = currentHit_vector.z();
1331 
1332  double nextHit_xyz[3];
1333  geometry_->CellInfo(nextHit->Plane(), nextHit->Cell(), NULL, nextHit_xyz);
1334  TVector3 nextHit_vector(nextHit_xyz);
1335 
1336  double nextHit_X = nextHit_vector.x();
1337  double nextHit_Z = nextHit_vector.z();
1338 
1339  double deltaXsquare = (nextHit_X - currentHit_X)*(nextHit_X - currentHit_X);
1340  double deltaZsquare = (nextHit_Z - currentHit_Z)*(nextHit_Z - currentHit_Z);
1341 
1342  double distance_between_hits = sqrt(deltaXsquare + deltaZsquare);
1343 
1344  x_clusteredTrack.cluster.Add(currentHit);
1345 
1346  // If the track's continuity is broken, any hit afterwards is ignored
1347  if (distance_between_hits > maxGap) { break; }
1348 
1349 
1350  }
1351 
1352  // Adding track to the vector x_goodTracks
1353  if (x_clusteredTrack.cluster.NCell() > track2D_min_hits_) x_clusteredTracks.push_back(x_clusteredTrack);
1354 
1355  }
1356 
1357  // Looping over YZ tracks
1358  for (auto y_track = hough_tracks[geo::kY].begin(); y_track != hough_tracks[geo::kY].end(); y_track++) {
1359 
1360  // Sorting hits in by increasing value of Z
1361  art::PtrVector<rb::CellHit> hits_planeSorted = y_track->cluster.AllCells();
1362  SortByPlane(hits_planeSorted);
1363 
1364  int trackSize = y_track->cluster.NCell();
1365  int halfTrackSize = trackSize/2;
1366 
1367  double y_slope = y_track->slope;
1368  double y_intercept = y_track->intercept;
1369 
1370  double y_theta = TMath::Abs(TMath::ATan(y_track->slope));
1371  double y_sinTheta = TMath::Sin(y_theta);
1372  double maxGap = max_distance_between_hits_ + max_distance_between_hits_ * y_sinTheta;
1373 
1374  if (event == eventToCheck) {
1375  std::cout << "YZ: NCell, theta, halftrack, maxGap: " << trackSize << ", " << y_theta * TMath::RadToDeg() << ", " << halfTrackSize << ", " << maxGap << std::endl;
1376  }
1377 
1378  Track2D y_clusteredTrack(geo::kY, y_slope, y_intercept);
1379 
1380  // From Z = halfTrackSize to Zmin
1381  for (int hit_i = halfTrackSize; hit_i > 1; hit_i--) {
1382 
1383  art::Ptr<rb::CellHit> currentHit = y_track->cluster.Cell(hit_i);
1384  art::Ptr<rb::CellHit> previousHit = y_track->cluster.Cell(hit_i - 1);
1385 
1386  double currentHit_xyz[3];
1387  geometry_->CellInfo(currentHit->Plane(), currentHit->Cell(), NULL, currentHit_xyz);
1388  TVector3 currentHit_vector(currentHit_xyz);
1389 
1390  double currentHit_Y = currentHit_vector.y();
1391  double currentHit_Z = currentHit_vector.z();
1392 
1393  double previousHit_xyz[3];
1394  geometry_->CellInfo(previousHit->Plane(), previousHit->Cell(), NULL, previousHit_xyz);
1395  TVector3 previousHit_vector(previousHit_xyz);
1396 
1397  double previousHit_Y = previousHit_vector.y();
1398  double previousHit_Z = previousHit_vector.z();
1399 
1400  double deltaYsquare = (previousHit_Y - currentHit_Y)*(previousHit_Y - currentHit_Y);
1401  double deltaZsquare = (previousHit_Z - currentHit_Z)*(previousHit_Z - currentHit_Z);
1402 
1403  double distance_between_hits = sqrt(deltaYsquare + deltaZsquare);
1404 
1405  y_clusteredTrack.cluster.Add(currentHit);
1406 
1407  // If the track's continuity is broken, any hit afterwards is ignored
1408  if (distance_between_hits > maxGap) { break; }
1409 
1410 
1411  }
1412 
1413  // From Z = halfTrackSize to Zmax
1414  for (int hit_i = halfTrackSize; hit_i < trackSize - 1; hit_i++) {
1415 
1416  art::Ptr<rb::CellHit> currentHit = y_track->cluster.Cell(hit_i);
1417  art::Ptr<rb::CellHit> nextHit = y_track->cluster.Cell(hit_i + 1);
1418 
1419  double currentHit_xyz[3];
1420  geometry_->CellInfo(currentHit->Plane(), currentHit->Cell(), NULL, currentHit_xyz);
1421  TVector3 currentHit_vector(currentHit_xyz);
1422 
1423  double currentHit_Y = currentHit_vector.y();
1424  double currentHit_Z = currentHit_vector.z();
1425 
1426  double nextHit_xyz[3];
1427  geometry_->CellInfo(nextHit->Plane(), nextHit->Cell(), NULL, nextHit_xyz);
1428  TVector3 nextHit_vector(nextHit_xyz);
1429 
1430  double nextHit_Y = nextHit_vector.y();
1431  double nextHit_Z = nextHit_vector.z();
1432 
1433  double deltaYsquare = (nextHit_Y - currentHit_Y)*(nextHit_Y - currentHit_Y);
1434  double deltaZsquare = (nextHit_Z - currentHit_Z)*(nextHit_Z - currentHit_Z);
1435 
1436  double distance_between_hits = sqrt(deltaYsquare + deltaZsquare);
1437 
1438  y_clusteredTrack.cluster.Add(currentHit);
1439 
1440  // If the track's continuity is broken, any hit afterwards is ignored
1441  if (distance_between_hits > maxGap) { break; }
1442 
1443 
1444  }
1445 
1446  // Adding track to the vector x_goodTracks
1447  if (y_clusteredTrack.cluster.NCell() > track2D_min_hits_) y_clusteredTracks.push_back(y_clusteredTrack);
1448 
1449  }
1450 
1451 
1452 
1453 
1454 
1455  //....................[ Checking if 3D tracks are still possible after removing far away hits ]..............
1456 
1457  if (x_clusteredTracks.size() == 0 || y_clusteredTracks.size() == 0) {
1458 
1459  ntracks3D = 0;
1460  //std::cerr << "++ Failed 3D reconstruction in Event " << event << std::endl;
1461 
1462  // Tracks Tree
1463  nhits.push_back(0);
1464  nplanesXZ.push_back(0);
1465  nplanesYZ.push_back(0);
1466  nplanes.push_back(0);
1467  startPlaneXZ.push_back(0);
1468  endPlaneXZ.push_back(0);
1469  startPlaneYZ.push_back(0);
1470  endPlaneYZ.push_back(0);
1471  len.push_back(0);
1472  startTime.push_back(0);
1473  meanTime.push_back(0);
1474  calE.push_back(0);
1475 
1476  startX.push_back(0);
1477  startY.push_back(0);
1478  startZ.push_back(0);
1479  endX.push_back(0);
1480  endY.push_back(0);
1481  endZ.push_back(0);
1482 
1483  dirX.push_back(0);
1484  dirY.push_back(0);
1485  dirZ.push_back(0);
1486 
1487  thetaXZ.push_back(0);
1488  thetaYZ.push_back(0);
1489 
1490  theta.push_back(0);
1491  phi.push_back(0);
1492  cosTheta.push_back(0);
1493  cosPhi.push_back(0);
1494 
1495  zenith.push_back(0);
1496  azimuth.push_back(0);
1497  cosZenith.push_back(0);
1498  cosAzimuth.push_back(0);
1499 
1500  tree_Tracks->Fill();
1501 
1502  // Events Tree
1503  totalntracksX = totalntracksX + ntracksXZ;
1504  totalntracksY = totalntracksY + ntracksYZ;
1505 
1506  // Getting the elapsed time (for the rate)
1507  if (event_counter_ == 0) {
1508 
1509  startYear = year;
1510  startMonth = month;
1511  startDay = day;
1512  startDoy = doy;
1513  startHour = hour;
1514  startMinute = minute;
1515  startSecond = second;
1516 
1517  UTCstartYear = UTCyear;
1518  UTCstartMonth = UTCmonth;
1519  UTCstartDay = UTCday;
1520  UTCstartDoy = UTCdoy;
1521  UTCstartHour = UTChour;
1522  UTCstartMinute = UTCminute;
1523  UTCstartSecond = UTCsecond;
1524 
1525  startTimeInSeconds = (startHour * 3600) + (startMinute * 60) + startSecond;
1526  }
1527 
1528  endYear = year;
1529  endMonth = month;
1530  endDay = day;
1531  endDoy = doy;
1532  endHour = hour;
1533  endMinute = minute;
1534  endSecond = second;
1535 
1536  UTCendYear = UTCyear;
1537  UTCendMonth = UTCmonth;
1538  UTCendDay = UTCday;
1539  UTCendDoy = UTCdoy;
1540  UTCendHour = UTChour;
1541  UTCendMinute = UTCminute;
1542  UTCendSecond = UTCsecond;
1543 
1544  tree_Events->Fill();
1545 
1546  event_counter_++;
1547  failed3Devents++;
1548 
1549  return;
1550 
1551  }
1552 
1553 
1554 
1555 
1556 
1557  //....................[ Removing rogue 2D tracks created by nearby hits around a track ]....................
1558 
1559  // We loop over tracks and hits, comparing each hit of a track with each hit of a previous track. If they are
1560  // a plane or less appart, tracks are crossing or VERY close to one another and, most likely, the track
1561  // with lesser hits is a rogue one
1562 
1563 
1564  // Sort Hough Tracks by NCell
1565  std::sort(x_clusteredTracks.begin(), x_clusteredTracks.end(), htk::sortByNCell);
1566  std::sort(y_clusteredTracks.begin(), y_clusteredTracks.end(), htk::sortByNCell);
1567 
1568  // Order Hough Tracks by descending order of NCell
1569  std::reverse(x_clusteredTracks.begin(), x_clusteredTracks.end());
1570  std::reverse(y_clusteredTracks.begin(), y_clusteredTracks.end());
1571 
1572 
1573  // Creating the vectors of good tracks
1574  std::vector<Track2D> x_goodTracks;
1575  std::vector<Track2D> y_goodTracks;
1576 
1577  if (event == eventToCheck){
1578 
1579  std::cerr << "Removing rogue tracks" << std::endl;
1580  std::cerr << "x_clusteredTrack, y_clusteredTrack sizes: " << x_clusteredTracks.size() << ", " << y_clusteredTracks.size() << std::endl;
1581 
1582 
1583  }
1584 
1585 
1586  bool isGood2DTrack;
1587 
1588 
1589  // Creating a matrix of track comparisons to check if the comparison at i-th and j-th
1590  // position returns passed or not the criteria
1591  bool **isGood2DxTrackMatrix;
1592  isGood2DxTrackMatrix = new bool*[x_clusteredTracks.size()];
1593 
1594  for (unsigned int p = 0; p < x_clusteredTracks.size(); p++) {
1595 
1596  isGood2DxTrackMatrix[p] = new bool[x_clusteredTracks.size()];
1597  }
1598 
1599 
1600 
1601  for (unsigned int n = 0; n < x_clusteredTracks.size(); n++) {
1602 
1603  for (unsigned int m = 0; m < x_clusteredTracks.size(); m++) {
1604 
1605  isGood2DxTrackMatrix[n][m] = true;
1606  }
1607  }
1608 
1609 
1610  // Looping over XZ tracks
1611  for (unsigned int i = 0; i < x_clusteredTracks.size() - 1; i++)
1612  {
1613 
1614  Track2D x_track_i = x_clusteredTracks[i];
1615  art::PtrVector<rb::CellHit> x_hits_i = x_track_i.cluster.AllCells();
1616 
1617  for (unsigned int j = i + 1; j < x_clusteredTracks.size(); j++)
1618  {
1619 
1620  Track2D x_track_j = x_clusteredTracks[j];
1621  art::PtrVector<rb::CellHit> x_hits_j = x_track_j.cluster.AllCells();
1622 
1623  for (auto const &hit_i: x_hits_i)
1624  {
1625  // Fetching hit_i Z position
1626  double hit_i_xyz[3];
1627  geometry_->CellInfo(hit_i->Plane(), hit_i->Cell(), NULL, hit_i_xyz);
1628  TVector3 hit_i_vector(hit_i_xyz);
1629 
1630  double hit_i_Z = hit_i_vector.z();
1631 
1632  for (auto const &hit_j: x_hits_j)
1633  {
1634  // Fetching hit_j Z position
1635  double hit_j_xyz[3];
1636  geometry_->CellInfo(hit_j->Plane(), hit_j->Cell(), NULL, hit_j_xyz);
1637  TVector3 hit_j_vector(hit_j_xyz);
1638 
1639  double hit_j_Z = hit_j_vector.z();
1640 
1641  if (TMath::Abs(hit_i_Z - hit_j_Z) < fake_tracks_threshold_) {
1642 
1643  isGood2DxTrackMatrix[i][j] = false;
1644  continue;
1645  }
1646  }
1647 
1648  // No need to verify more hits if the false combination was already triggered
1649  if (isGood2DxTrackMatrix[i][j] == false) { continue; }
1650  }
1651  }
1652  }
1653 
1654 
1655  // Filling up the x_goodTracks vector
1656  for (unsigned int n = 0; n < x_clusteredTracks.size(); n++) {
1657 
1658  isGood2DTrack = true;
1659 
1660  for (unsigned int m = 0; m < x_clusteredTracks.size(); m++) {
1661 
1662  if (isGood2DxTrackMatrix[m][n] == false) { isGood2DTrack = false; }
1663 
1664  }
1665 
1666  if (isGood2DTrack == true) { x_goodTracks.push_back(x_clusteredTracks[n]); }
1667 
1668  }
1669 
1670 
1671  //..........................
1672 
1673 
1674  // Matrix for YZ tracks. Starts as a matrix filled with true
1675  bool **isGood2DyTrackMatrix;
1676  isGood2DyTrackMatrix = new bool*[y_clusteredTracks.size()];
1677 
1678  for (unsigned int p = 0; p < y_clusteredTracks.size(); p++) {
1679 
1680  isGood2DyTrackMatrix[p] = new bool[y_clusteredTracks.size()];
1681  }
1682 
1683  for (unsigned int n = 0; n < y_clusteredTracks.size(); n++) {
1684 
1685  for (unsigned int m = 0; m < y_clusteredTracks.size(); m++) {
1686 
1687  isGood2DyTrackMatrix[n][m] = true;
1688  }
1689  }
1690 
1691 
1692  // Looping over YZ tracks
1693  for (unsigned int i = 0; i < y_clusteredTracks.size() - 1; i++)
1694  {
1695 
1696  Track2D y_track_i = y_clusteredTracks[i];
1697  art::PtrVector<rb::CellHit> y_hits_i = y_track_i.cluster.AllCells();
1698 
1699  for (unsigned int j = i + 1; j < y_clusteredTracks.size(); j++)
1700  {
1701 
1702  Track2D y_track_j = y_clusteredTracks[j];
1703  art::PtrVector<rb::CellHit> y_hits_j = y_track_j.cluster.AllCells();
1704 
1705  for (auto const &hit_i: y_hits_i)
1706  {
1707  // Fetching hit_i Z position
1708  double hit_i_xyz[3];
1709  geometry_->CellInfo(hit_i->Plane(), hit_i->Cell(), NULL, hit_i_xyz);
1710  TVector3 hit_i_vector(hit_i_xyz);
1711 
1712  double hit_i_Z = hit_i_vector.z();
1713 
1714  for (auto const &hit_j: y_hits_j)
1715  {
1716  // Fetching hit_j Z position
1717  double hit_j_xyz[3];
1718  geometry_->CellInfo(hit_j->Plane(), hit_j->Cell(), NULL, hit_j_xyz);
1719  TVector3 hit_j_vector(hit_j_xyz);
1720 
1721  double hit_j_Z = hit_j_vector.z();
1722 
1723  if (TMath::Abs(hit_i_Z - hit_j_Z) < fake_tracks_threshold_) {
1724 
1725  isGood2DyTrackMatrix[i][j] = false;
1726  continue;
1727  }
1728  }
1729 
1730  // No need to verify more hits if the false combination was already triggered
1731  if (isGood2DyTrackMatrix[i][j] == false) { continue; }
1732  }
1733  }
1734  }
1735 
1736 
1737  // Filling up the y_goodTracks vector
1738  for (unsigned int n = 0; n < y_clusteredTracks.size(); n++) {
1739 
1740  isGood2DTrack = true;
1741 
1742  for (unsigned int m = 0; m < y_clusteredTracks.size(); m++) {
1743 
1744  if (isGood2DyTrackMatrix[m][n] == false) { isGood2DTrack = false; }
1745 
1746  }
1747 
1748  if (isGood2DTrack == true) { y_goodTracks.push_back(y_clusteredTracks[n]); }
1749 
1750  }
1751 
1752 
1753 
1754 
1755 
1756  // PRINTING THE MATRICES AND TRACK INFO
1757  if (event == eventToCheck) {
1758 
1759  std::cerr << std::endl;
1760  std::cerr << std::endl;
1761  std::cerr << "XZ MATRIX" << std::endl;
1762  std::cerr << std::endl;
1763  std::cerr << "TRACKS ";
1764 
1765  for (unsigned int n = 0; n < x_clusteredTracks.size(); n++) {
1766 
1767  std::cerr << n << " ";
1768 
1769  }
1770 
1771  std::cerr << std::endl;
1772 
1773  for (unsigned int n = 0; n < x_clusteredTracks.size(); n++) {
1774 
1775  std::cerr << " " << n << " | ";
1776 
1777  for (unsigned int m = 0; m < x_clusteredTracks.size(); m++) {
1778 
1779  if (m == 0) { std::cerr << isGood2DxTrackMatrix[n][m]; }
1780  if (m > 0) { std::cerr << " " << isGood2DxTrackMatrix[n][m]; }
1781 
1782  }
1783 
1784  std::cerr << " |" << std::endl;
1785 
1786  }
1787 
1788  //....................
1789 
1790 
1791  std::cerr << std::endl;
1792  std::cerr << std::endl;
1793  std::cerr << "YZ MATRIX" << std::endl;
1794  std::cerr << std::endl;
1795  std::cerr << "TRACKS ";
1796 
1797  for (unsigned int n = 0; n < y_clusteredTracks.size(); n++) {
1798 
1799  std::cerr << n << " ";
1800 
1801  }
1802 
1803  std::cerr << std::endl;
1804 
1805  for (unsigned int n = 0; n < y_clusteredTracks.size(); n++) {
1806 
1807  std::cerr << " " << n << " | ";
1808 
1809  for (unsigned int m = 0; m < y_clusteredTracks.size(); m++) {
1810 
1811  if (m == 0) { std::cerr << isGood2DyTrackMatrix[n][m]; }
1812  if (m > 0) { std::cerr << " " << isGood2DyTrackMatrix[n][m]; }
1813 
1814  }
1815 
1816  std::cerr << " |" << std::endl;
1817  }
1818 
1819  }
1820 
1821 
1822 
1823 
1824  //....................[ Creating 3D tracks ]....................
1825 
1826  ntracks3D = 0;
1827 
1828  double x_hit_startX = 9999;
1829  double x_hit_startZ = 9999;
1830 
1831  int xz_nplanes = 0;
1832  int xz_startPlane = 0;
1833  int xz_endPlane = 0;
1834  double xz_hit_start_time = 9999;
1835 
1836 
1837  double y_hit_startY = 9999;
1838  double y_hit_startZ = 9999;
1839 
1840  int yz_nplanes = 0;
1841  int yz_startPlane = 0;
1842  int yz_endPlane = 0;
1843  double yz_hit_start_time = 9999;
1844 
1845 
1846  double x_hit_endX = -9999;
1847  double x_hit_endZ = -9999;
1848 
1849  double xz_hit_end_time = 9999;
1850 
1851  double y_hit_endY = -9999;
1852  double y_hit_endZ = -9999;
1853 
1854  double yz_hit_end_time = 9999;
1855 
1856  // Sort Hough Tracks by NCell
1857  std::sort(x_goodTracks.begin(), x_goodTracks.end(), htk::sortByNCell);
1858  std::sort(y_goodTracks.begin(), y_goodTracks.end(), htk::sortByNCell);
1859 
1860  // Order Hough Tracks by descending order of NCell
1861  std::reverse(x_goodTracks.begin(), x_goodTracks.end());
1862  std::reverse(y_goodTracks.begin(), y_goodTracks.end());
1863 
1864 
1865 
1866  std::vector<Track_Match> track_matches;
1867  std::unique_ptr<std::vector<Track3D> > tracks_3D(new std::vector<Track3D>);
1868 
1869 
1870 
1871  // Loop over XZ tracks
1872  for (auto x_track = x_goodTracks.begin(); x_track != x_goodTracks.end(); x_track++)
1873  {
1874 
1875  art::PtrVector<rb::CellHit> x_hits_planeSorted = x_track->cluster.AllCells();
1876  SortByPlane(x_hits_planeSorted);
1877 
1878 
1879  double x_track_meanTime = 0;
1880  unsigned int nhitsXZ = 0;
1881  int this_xz_plane = -99;
1882 
1883  // Loop over XZ hits
1884  for (auto const& x_hit : x_hits_planeSorted)
1885  {
1886 
1887  // Fetching track's start.z/end.z positions
1888  double x_hit_xyz[3];
1889  geometry_->CellInfo(x_hit->Plane(), x_hit->Cell(), NULL, x_hit_xyz);
1890  TVector3 x_hit_vector(x_hit_xyz);
1891 
1892  if (nhitsXZ == 0) {
1893 
1894  x_hit_startX = x_hit_vector.x();
1895  x_hit_startZ = x_hit_vector.z();
1896 
1897  xz_startPlane = x_hit->Plane();
1898  xz_hit_start_time = x_hit->TNS();
1899 
1900  }
1901 
1902  if (nhitsXZ == x_track->cluster.NCell() - 1) {
1903 
1904  x_hit_endX = x_hit_vector.x();
1905  x_hit_endZ = x_hit_vector.z();
1906 
1907  xz_endPlane = x_hit->Plane();
1908  xz_hit_end_time = x_hit->TNS();
1909 
1910  }
1911 
1912  if (x_hit->Plane() != this_xz_plane) {
1913 
1914  this_xz_plane = x_hit->Plane();
1915  xz_nplanes++;
1916 
1917  }
1918 
1919  x_track_meanTime = x_track_meanTime + x_hit->TNS();
1920 
1921  nhitsXZ++;
1922  }
1923 
1924 
1925 
1926  // Loop over YZ tracks
1927  for (auto y_track = y_goodTracks.begin(); y_track != y_goodTracks.end(); y_track++)
1928  {
1929 
1930  art::PtrVector<rb::CellHit> y_hits_planeSorted = y_track->cluster.AllCells();
1931  SortByPlane(y_hits_planeSorted);
1932 
1933  double y_track_meanTime = 0;
1934  unsigned int nhitsYZ = 0;
1935  int this_yz_plane = -99;
1936 
1937  // Loop over YZ hits
1938  for (auto const& y_hit : y_hits_planeSorted)
1939  {
1940 
1941  // Fetching track's start.z/end.z positions
1942  double y_hit_xyz[3];
1943  geometry_->CellInfo(y_hit->Plane(), y_hit->Cell(), NULL, y_hit_xyz);
1944  TVector3 y_hit_vector(y_hit_xyz);
1945 
1946  if (nhitsYZ == 0) {
1947 
1948  y_hit_startY = y_hit_vector.y();
1949  y_hit_startZ = y_hit_vector.z();
1950 
1951  yz_startPlane = y_hit->Plane();
1952  yz_hit_start_time = y_hit->TNS();
1953 
1954 
1955  }
1956 
1957  if (nhitsYZ == y_track->cluster.NCell() - 1) {
1958 
1959  y_hit_endY = y_hit_vector.y();
1960  y_hit_endZ = y_hit_vector.z();
1961 
1962  yz_endPlane = y_hit->Plane();
1963  yz_hit_end_time = y_hit->TNS();
1964 
1965 
1966  }
1967 
1968  if (y_hit->Plane() != this_yz_plane) {
1969 
1970  this_yz_plane = y_hit->Plane();
1971  yz_nplanes++;
1972 
1973  }
1974 
1975  y_track_meanTime = y_track_meanTime + y_hit->TNS();
1976 
1977  nhitsYZ++;
1978  }
1979 
1980  x_track_meanTime = x_track_meanTime / x_track->cluster.NCell();
1981  y_track_meanTime = y_track_meanTime / y_track->cluster.NCell();
1982 
1983 
1984  // Track 3D selection criteria
1985  if (nhitsXZ + nhitsYZ < track3D_min_hits_) continue;
1986 
1987 
1988 
1989  if (event == eventToCheck){
1990  std::cerr << std::endl;
1991  std::cerr << "Track " << ntracks3D << " (XZ) (YZ)" << std::endl;
1992  std::cerr << "Time : " << x_track_meanTime << ", " << y_track_meanTime << std::endl;
1993  std::cerr << "NCell : " << x_track->cluster.NCell() << ", " << y_track->cluster.NCell() << std::endl;
1994  std::cerr << "Start : (" << x_hit_startX << ", " << x_hit_startZ << "), (" << y_hit_startY << ", " << y_hit_startZ << ")" << std::endl;
1995  std::cerr << "End : (" << x_hit_endX << ", " << x_hit_endZ << "), (" << y_hit_endY << ", " << y_hit_endZ << ")" << std::endl;
1996  std::cerr << "Theta : " << TMath::ATan(x_track->slope) * TMath::RadToDeg() << ", " << TMath::ATan(y_track->slope) * TMath::RadToDeg() << std::endl;
1997  std::cerr << std::endl;
1998  }
1999 
2000 
2001  // i) The start.z/end.z of the 2D tracks should coincide in a good 3D tracks
2002  if (TMath::Abs(x_hit_startZ - y_hit_startZ) > vertex_max_distance_ || TMath::Abs(x_hit_endZ - y_hit_endZ) > vertex_max_distance_) continue;
2003 
2004 
2005 
2006  // Verifying if it is indeed a good track match by doing a few checks
2007  bool isGood3DTrack = true;
2008 
2009  if (ntracks3D > 0) {
2010 
2011  for (unsigned int i = 0; i < startZ.size(); i++) {
2012 
2013  // i) XZ and YZ slopes cannot be repeated
2014  if (thetaXZ.at(i) == (TMath::ATan(x_track->slope) * TMath::RadToDeg()) || thetaYZ.at(i) == (TMath::ATan(y_track->slope) * TMath::RadToDeg()) ) {
2015 
2016  isGood3DTrack = false;
2017  continue;
2018  }
2019  }
2020  }
2021 
2022 
2023  if (isGood3DTrack == true) {
2024 
2025  // Since hits were only sorted by plane, we need to guarantee that we are not looking at the track backward in time
2026  // Amount of upward muons is neglectible and every track will be considered to be going downward
2027  // The only information stored in the ntuple are the start/end points, thus, I will only flip them if startY < endY
2028  if (y_hit_startY < y_hit_endY) {
2029 
2030  double hit_temp;
2031 
2032  //YZ view, Y
2033  hit_temp = y_hit_startY;
2034  y_hit_startY = y_hit_endY;
2035  y_hit_endY = hit_temp;
2036 
2037  //YZ view, Z
2038  hit_temp = y_hit_startZ;
2039  y_hit_startZ = y_hit_endZ;
2040  y_hit_endZ = hit_temp;
2041 
2042  //YZ view, time
2043  hit_temp = yz_hit_start_time;
2044  yz_hit_start_time = yz_hit_end_time;
2045  yz_hit_end_time = hit_temp;
2046 
2047 
2048  if ( (y_hit_startZ - y_hit_endZ) > 0 && (x_hit_startZ - x_hit_endZ) < 0 ) {
2049 
2050  //XZ view, X
2051  hit_temp = x_hit_startX;
2052  x_hit_startX = x_hit_endX;
2053  x_hit_endX = hit_temp;
2054 
2055  //XZ view, Z
2056  hit_temp = x_hit_startZ;
2057  x_hit_startZ = x_hit_endZ;
2058  x_hit_endZ = hit_temp;
2059 
2060  //XZ view, time and planes
2061  hit_temp = xz_hit_start_time;
2062  xz_hit_start_time = xz_hit_end_time;
2063  xz_hit_end_time = hit_temp;
2064 
2065  hit_temp = xz_startPlane;
2066  xz_startPlane = xz_endPlane;
2067  xz_endPlane = hit_temp;
2068 
2069  }
2070 
2071  if ( (y_hit_startZ - y_hit_endZ) < 0 && (x_hit_startZ - x_hit_endZ) > 0 ) {
2072 
2073  //XZ view, X
2074  hit_temp = x_hit_startX;
2075  x_hit_startX = x_hit_endX;
2076  x_hit_endX = hit_temp;
2077 
2078  //XZ view, Z
2079  hit_temp = x_hit_startZ;
2080  x_hit_startZ = x_hit_endZ;
2081  x_hit_endZ = hit_temp;
2082 
2083  //XZ view, time and planes
2084  hit_temp = xz_hit_start_time;
2085  xz_hit_start_time = xz_hit_end_time;
2086  xz_hit_end_time = hit_temp;
2087 
2088  hit_temp = xz_startPlane;
2089  xz_startPlane = xz_endPlane;
2090  xz_endPlane = hit_temp;
2091  }
2092 
2093  }
2094 
2095 
2096 
2097  // Building Zenith and Azimuth angles
2098  double Rx = x_hit_endX - x_hit_startX;
2099  double Ry = y_hit_endY - y_hit_startY;
2100  double Rz = y_hit_endZ - y_hit_startZ;
2101  double absR = TMath::Sqrt( Rx*Rx + Ry*Ry + Rz*Rz );
2102 
2103  double trackTheta = TMath::ACos( Ry / absR );
2104  double zen = TMath::Pi() - trackTheta;
2105  double trackPhi = TMath::ATan(Rz / Rx); // Starts at the X+ axis
2106  double azi = 0;
2107 
2108  if (Rx > 0 && Rz > 0) { azi = trackPhi; }
2109  if (Rx < 0 && Rz > 0) { azi = TMath::Pi() + trackPhi; }
2110  if (Rx < 0 && Rz < 0) { azi = TMath::Pi() + trackPhi; }
2111  if (Rx > 0 && Rz < 0) { azi = 2*TMath::Pi() + trackPhi; }
2112  if (Rx == 0 && Rz > 0) { azi = TMath::Pi()/2; }
2113  if (Rx == 0 && Rz < 0) { azi = 3*TMath::Pi()/2; }
2114  if (Rx > 0 && Rz == 0) { azi = 0; }
2115  if (Rx < 0 && Rz == 0) { azi = TMath::Pi(); }
2116 
2117  Track3D track3D(*x_track, *y_track);
2118 
2119  double track3DmeanTime = (x_track_meanTime + y_track_meanTime) / 2;
2120 
2121  // Filling the Tracks Tree
2122  nhits.push_back(track3D.cluster().NCell());
2123  nplanesXZ.push_back(xz_nplanes);
2124  nplanesYZ.push_back(yz_nplanes);
2125  nplanes.push_back(xz_nplanes + yz_nplanes);
2126  startPlaneXZ.push_back(xz_startPlane);
2127  endPlaneXZ.push_back(xz_endPlane);
2128  startPlaneYZ.push_back(yz_startPlane);
2129  endPlaneYZ.push_back(yz_endPlane);
2130 
2131  len.push_back(absR);
2132 
2133  //std::cout << "TRACK LEN: " << absR << std::endl;
2134 
2135 
2136  if (xz_hit_start_time < yz_hit_start_time) { startTime.push_back(xz_hit_start_time); }
2137  else { startTime.push_back(yz_hit_start_time); }
2138  meanTime.push_back(track3DmeanTime);
2139 
2140  calE.push_back(track3D.cluster().TotalGeV());
2141 
2142  startX.push_back(x_hit_startX);
2143  startY.push_back(y_hit_startY);
2144  if (x_hit_startZ > y_hit_startZ) startZ.push_back(x_hit_startZ);
2145  else startZ.push_back(y_hit_startZ);
2146 
2147 
2148  endX.push_back(x_hit_endX);
2149  endY.push_back(y_hit_endY);
2150  if (x_hit_endZ > y_hit_endZ) endZ.push_back(x_hit_endZ);
2151  else endZ.push_back(y_hit_endZ);
2152 
2153  dirX.push_back( Rx/absR );
2154  dirY.push_back( Ry/absR );
2155  dirZ.push_back( Rz/absR );
2156 
2157  thetaXZ.push_back( TMath::ATan(x_track->slope) * TMath::RadToDeg() ); // starting from Z towards X
2158  thetaYZ.push_back( TMath::ATan(y_track->slope) * TMath::RadToDeg() ); // starting from Z towards Y
2159 
2160  theta.push_back( zen * TMath::RadToDeg() );
2161  phi.push_back( azi * TMath::RadToDeg() );
2162  cosTheta.push_back( TMath::Cos(zen) );
2163  cosPhi.push_back( TMath::Cos(azi) );
2164 
2165  // ND angle from true north (http://ppd.fnal.gov/align/beams-doc-1148.html)
2166  double trueNorthCorrection = 0.668112; // in rad (equal to 38.28 deg)
2167 
2168  zenith.push_back( zen * TMath::RadToDeg() );
2169  if (azi - trueNorthCorrection > 0) { azimuth.push_back( (azi - trueNorthCorrection) * TMath::RadToDeg() ); }
2170  if (azi - trueNorthCorrection < 0) { azimuth.push_back( (TMath::Pi()*2 - TMath::Abs(azi - trueNorthCorrection)) * TMath::RadToDeg() ); }
2171  cosZenith.push_back( TMath::Cos(zen));
2172  cosAzimuth.push_back( TMath::Cos(azi - trueNorthCorrection));
2173 
2174  // Adding objects to vectors
2175  track_matches.emplace_back(x_track, y_track);
2176  tracks_3D->push_back(track3D);
2177 
2178  if (event == eventToCheck){
2179  std::cerr << std::endl;
2180  std::cerr << "PASSED TRACK" << std::endl;
2181  std::cerr << "Track " << ntracks3D << " (XZ) (YZ)" << std::endl;
2182  std::cerr << "Time : " << x_track_meanTime << ", " << y_track_meanTime << std::endl;
2183  std::cerr << "NCell : " << x_track->cluster.NCell() << ", " << y_track->cluster.NCell() << std::endl;
2184  std::cerr << "Start : (" << x_hit_startX << ", " << x_hit_startZ << "), (" << y_hit_startY << ", " << y_hit_startZ << ")" << std::endl;
2185  std::cerr << "End : (" << x_hit_endX << ", " << x_hit_endZ << "), (" << y_hit_endY << ", " << y_hit_endZ << ")" << std::endl;
2186  std::cerr << "Theta : " << TMath::ATan(x_track->slope) * TMath::RadToDeg() << ", " << TMath::ATan(y_track->slope) * TMath::RadToDeg() << std::endl;
2187  std::cerr << std::endl;
2188  }
2189 
2190  ntracks3D++;
2191 
2192  }
2193  }
2194 
2195  if (event == eventToCheck) std::cerr << std::endl;
2196  }
2197 
2198  //....................[ End of 3D tracks loop ]....................
2199 
2200 
2201  // Fill the Tracks Tree branches with zeros when there are no good 3D tracks
2202  if (ntracks3D == 0) {
2203 
2204  //std::cerr << "++ Failed 3D reconstruction in Event " << event << std::endl;
2205 
2206 
2207  nhits.push_back(0);
2208  nplanesXZ.push_back(0);
2209  nplanesYZ.push_back(0);
2210  nplanes.push_back(0);
2211  len.push_back(0);
2212  startTime.push_back(0);
2213  meanTime.push_back(0);
2214  calE.push_back(0);
2215 
2216  startX.push_back(0);
2217  startY.push_back(0);
2218  startZ.push_back(0);
2219 
2220  endX.push_back(0);
2221  endY.push_back(0);
2222  endZ.push_back(0);
2223 
2224  dirX.push_back(0);
2225  dirY.push_back(0);
2226  dirZ.push_back(0);
2227 
2228  thetaXZ.push_back(0);
2229  thetaYZ.push_back(0);
2230 
2231  theta.push_back(0);
2232  phi.push_back(0);
2233  cosTheta.push_back(0);
2234  cosPhi.push_back(0);
2235 
2236  zenith.push_back(0);
2237  azimuth.push_back(0);
2238  cosZenith.push_back(0);
2239  cosAzimuth.push_back(0);
2240  }
2241 
2242 
2243  tree_Tracks->Fill();
2244 
2245  totalntracksX = totalntracksX + ntracksXZ;
2246  totalntracksY = totalntracksY + ntracksYZ;
2247  total3DTracks = total3DTracks + ntracks3D;
2248 
2249 
2250 
2251  // Getting the elapsed time (for the rate)
2252  if (event_counter_ == 0) {
2253 
2254  startYear = year;
2255  startMonth = month;
2256  startDay = day;
2257  startDoy = doy;
2258  startHour = hour;
2259  startMinute = minute;
2260  startSecond = second;
2261 
2262  UTCstartYear = UTCyear;
2263  UTCstartMonth = UTCmonth;
2264  UTCstartDay = UTCday;
2265  UTCstartDoy = UTCdoy;
2266  UTCstartHour = UTChour;
2267  UTCstartMinute = UTCminute;
2268  UTCstartSecond = UTCsecond;
2269 
2270  startTimeInSeconds = (startHour * 3600) + (startMinute * 60) + startSecond;
2271  }
2272 
2273  endYear = year;
2274  endMonth = month;
2275  endDay = day;
2276  endDoy = doy;
2277  endHour = hour;
2278  endMinute = minute;
2279  endSecond = second;
2280 
2281  UTCendYear = UTCyear;
2282  UTCendMonth = UTCmonth;
2283  UTCendDay = UTCday;
2284  UTCendDoy = UTCdoy;
2285  UTCendHour = UTChour;
2286  UTCendMinute = UTCminute;
2287  UTCendSecond = UTCsecond;
2288 
2289  // Getting the the number of single/multi muon events
2290  if (ntracks3D == 1) nSinglemuEvents++;
2291  if (ntracks3D > 1) nMultimuEvents++;
2292 
2293 
2294  tree_Events->Fill();
2295 
2296  event_counter_++;
2297 
2298 
2299  // Deleting variables
2300  for (unsigned int k = 0; k < x_clusteredTracks.size(); k++) {
2301 
2302  delete[] isGood2DxTrackMatrix[k];
2303 
2304  }
2305 
2306  delete[] isGood2DxTrackMatrix;
2307  isGood2DxTrackMatrix = 0;
2308 
2309  for (unsigned int k = 0; k < y_clusteredTracks.size(); k++) {
2310 
2311  delete[] isGood2DyTrackMatrix[k];
2312 
2313  }
2314 
2315  delete[] isGood2DyTrackMatrix;
2316  isGood2DyTrackMatrix = 0;
2317 
2318 
2319 }
2320 
2321 
2322 
2323 //........................................[ SORT TRACKS BY NCELL ]........................................
2324 
2325 bool htk::sortByNCell(htk::Track2D &track_i, htk::Track2D &track_j) {
2326 
2327  return (track_i.cluster.NCell() < track_j.cluster.NCell());
2328 }
2329 
2330 
2331 
2332 //........................................[ SORT TRACKS BY FLSHITS ]........................................
2333 
2334 /*
2335  bool htk::sortByFLSHits(sim::Particle *particle_i, sim::Particle *particle_j) {
2336 
2337  return (art::ServiceHandle<cheat::BackTracker> bt->ParticleNavigator()->ParticleToFLSHit(particle_i->TrackId()).size()
2338  <
2339  art::ServiceHandle<cheat::BackTracker> bt->ParticleNavigator()->ParticleToFLSHit(particle_j->TrackId()).size() );
2340  }
2341  */
2342 
2343 
2344 
2345 
2346 //........................................[ END JOB ]........................................
2347 
2348 void htk::HoughTrack::endJob() {
2349 
2350 
2351 
2352 
2353 
2354  // Calculating the rate
2355  endTimeInSeconds = (endHour * 3600) + (endMinute * 60) + endSecond;
2356 
2357  if (startDoy == endDoy) { exposureTime = endTimeInSeconds - startTimeInSeconds; }
2358  if (startDoy != endDoy) { exposureTime = (86400 - startTimeInSeconds) + endTimeInSeconds; }
2359 
2360  singlemuRate = nSinglemuEvents / exposureTime;
2361  multimuRate = nMultimuEvents / exposureTime;
2362 
2363 
2364  // Filling the Rates tree
2365  tree_Rates->Fill();
2366 
2367 
2368  std::cout << std::endl;
2369  std::cout << std::endl;
2370  std::cout << "End Job" << std::endl;
2371  std::cout << std::endl;
2372  std::cout << std::endl;
2373  std::cout << std::endl;
2374  std::cout << "------ Summary ------" << std::endl;
2375  std::cout << "Total events : " << event_counter_ << std::endl;
2376  std::cout << "Total 2D Hough tracks : " << totalntracksX << " (X) + " << totalntracksY << " (Y) = "<< totalntracksX + totalntracksY << std::endl;
2377  std::cout << "Total 3D Hough tracks : " << total3DTracks << std::endl;
2378  std::cout << "Failed 3D merging : " << failed3Devents << std::endl;
2379  std::cout << std::endl;
2380  std::cout << "Total multimuon events: " << nMultimuEvents << std::endl;
2381  std::cout << "Start Time : " << startHour << ":" << startMinute << ":" << startSecond << " (" << startTimeInSeconds << " s)" << std::endl;
2382  std::cout << "End Time : " << endHour << ":" << endMinute << ":" << endSecond << " (" << endTimeInSeconds << " s)" << std::endl;
2383  std::cout << "Exposure Time : " << exposureTime << " s" << std::endl;
2384  std::cout << std::endl;
2385  std::cout << std::endl;
2386 
2387 
2388 }
2389 
2390 
2391 
2392 DEFINE_ART_MODULE(htk::HoughTrack)
len
const XML_Char int len
Definition: expat.h:262
htk::HoughTrack::triggerStartUnixTime
unsigned long long triggerStartUnixTime
Definition: HoughTrack_module.cc:162
htk::HoughTrack::particleTheta
std::vector< double > particleTheta
Definition: HoughTrack_module.cc:242
art::Event::subRun
SubRunNumber_t subRun() const
Definition: Event.h:72
htk::HoughTrack::nplanesYZ
std::vector< double > nplanesYZ
Definition: HoughTrack_module.cc:188
BackTracker.h
back track the reconstruction to the simulation
simb::MCParticle::Position
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:211
htk::HoughTrack::calE
std::vector< double > calE
Definition: HoughTrack_module.cc:197
htk::HoughTrack::particleEndX
std::vector< double > particleEndX
Definition: HoughTrack_module.cc:234
rb::Cluster::NCell
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
htk::HoughTrack::cosAzimuth
std::vector< double > cosAzimuth
Definition: HoughTrack_module.cc:216
it
set< int >::iterator it
Definition: ProduceBadRunsByFailureMode.C:101
htk::HoughTrack::particleDirX
std::vector< double > particleDirX
Definition: HoughTrack_module.cc:237
htk::HoughTrack::particleNplanesXZ
std::vector< double > particleNplanesXZ
Definition: HoughTrack_module.cc:225
htk::HoughTrack::particleStartY
std::vector< double > particleStartY
Definition: HoughTrack_module.cc:232
febshutoff_auto.end
end
Definition: febshutoff_auto.py:54
cheat::BackTracker::ParticleNavigator
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
htk::HoughTrack::endZ
std::vector< double > endZ
Definition: HoughTrack_module.cc:203
htk::Track_Match::x_track
std::vector< Track2D >::const_iterator x_track
Definition: HoughTrack_module.cc:101
htk::HoughTrack::cosZenith
std::vector< double > cosZenith
Definition: HoughTrack_module.cc:215
htk::HoughTrack::startZ
std::vector< double > startZ
Definition: HoughTrack_module.cc:200
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
geo::CellGeo::LocalToWorld
void LocalToWorld(const double *local, double *world) const
Definition: CellGeo.cxx:80
rb::CellHit::Plane
unsigned short Plane() const
Definition: CellHit.h:39
htk::Track2D::slope
double slope
Definition: Track2D.h:26
htk::HoughTrack::endPlaneYZ
std::vector< double > endPlaneYZ
Definition: HoughTrack_module.cc:193
htk::HoughTrack::dirZ
std::vector< double > dirZ
Definition: HoughTrack_module.cc:206
htk::Track3D::cluster
rb::Cluster cluster() const
Definition: Track3D.cxx:113
geo::PlaneGeo::Cell
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
p
const char * p
Definition: xmltok.h:285
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
cheat::BackTracker::ParticleToFLSHit
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
Definition: BackTracker_service.cc:1514
art::Timestamp::timeHigh
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:35
simb::MCParticle::Mass
double Mass() const
Definition: MCParticle.h:238
sim::ParticleNavigator::const_iterator
list_type::const_iterator const_iterator
Definition: ParticleNavigator.h:86
htk::Track2D::cluster
rb::Cluster cluster
Definition: Track2D.h:25
htk::HoughTrack::particleCosPhi
std::vector< double > particleCosPhi
Definition: HoughTrack_module.cc:245
geo::kX
Vertical planes which measure X.
Definition: PlaneGeo.h:28
htk::HoughTrack::endY
std::vector< double > endY
Definition: HoughTrack_module.cc:202
geo::PlaneGeo::Ncells
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
htk
Definition: Constants.h:4
om::cerr
OStream cerr
Definition: OStream.cxx:7
htk::HoughTrack::particlePx
std::vector< double > particlePx
Definition: HoughTrack_module.cc:227
htk::HoughTrack::cosTheta
std::vector< double > cosTheta
Definition: HoughTrack_module.cc:211
site_stats_from_log.get
get
Definition: site_stats_from_log.py:189
geo::GeometryBase::Plane
const PlaneGeo * Plane(unsigned int i) const
Definition: GeometryBase.cxx:749
htk::HoughTrack::particleEndZ
std::vector< double > particleEndZ
Definition: HoughTrack_module.cc:236
htk::HoughTrack::particleE
std::vector< double > particleE
Definition: HoughTrack_module.cc:230
DEFINE_ART_MODULE
DEFINE_ART_MODULE(TestTMapFile)
htk::HoughTrack::meanTime
std::vector< double > meanTime
Definition: HoughTrack_module.cc:196
novadaq::timeutils::convertNovaTimeToUnixTime
bool convertNovaTimeToUnixTime(uint64_t const &inputNovaTime, struct timespec &outputUnixTime)
Definition: TimingUtilities.cpp:163
htk::HoughTrack::endX
std::vector< double > endX
Definition: HoughTrack_module.cc:201
art::Timestamp::value
constexpr TimeValue_t value() const
Definition: Timestamp.h:24
geo::kY
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
htk::HoughTrack::startX
std::vector< double > startX
Definition: HoughTrack_module.cc:198
htk::sortByFLSHits
bool sortByFLSHits(sim::Particle particle_i, sim::Particle particle_j)
site_stats_from_log.reverse
reverse
Definition: site_stats_from_log.py:189
kY
const Var kY([](const caf::SRProxy *sr){float tmp=0.f;if(sr->mc.nu.empty()) return tmp;tmp=sr->mc.nu[0].y;return tmp;})
htk::HoughTrack::particleNplanesYZ
std::vector< double > particleNplanesYZ
Definition: HoughTrack_module.cc:226
geo::GeometryBase::CellInfo
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
Definition: GeometryBase.cxx:813
htk::HoughTrack::event
unsigned long event
Definition: HoughTrack_module.cc:160
rb::Cluster::AllCells
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
submit_syst.y
y
Definition: submit_syst.py:249
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:209
rb::Cluster::Add
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
htk::HoughTrack::particleEndY
std::vector< double > particleEndY
Definition: HoughTrack_module.cc:235
htk::HoughTrack::startPlaneXZ
std::vector< double > startPlaneXZ
Definition: HoughTrack_module.cc:190
m
#define m
Definition: MakeGENIEReweightTable_module.cc:407
htk::HoughTrack::startTime
std::vector< double > startTime
Definition: HoughTrack_module.cc:195
rb::CellHit::Cell
unsigned short Cell() const
Definition: CellHit.h:40
htk::Track_Match
Definition: HoughTrack_module.cc:96
submit_syst.x
x
Definition: submit_syst.py:249
htk::HoughTrack::failed3Devents
unsigned int failed3Devents
Definition: HoughTrack_module.cc:154
breakpoints::beginJob
void beginJob()
htk::HoughTrack::totalntracksX
unsigned int totalntracksX
Definition: HoughTrack_module.cc:151
htk::HoughTrack::totalntracksY
unsigned int totalntracksY
Definition: HoughTrack_module.cc:152
htk::HoughTrack::particleMass
std::vector< double > particleMass
Definition: HoughTrack_module.cc:223
htk::HoughTrack::particleFLSHits
std::vector< double > particleFLSHits
Definition: HoughTrack_module.cc:224
htk::HoughTrack::nplanes
std::vector< double > nplanes
Definition: HoughTrack_module.cc:189
htk::HoughTrack::slice_label_
std::string slice_label_
Definition: HoughTrack_module.cc:132
htk::HoughTrack::bt
art::ServiceHandle< cheat::BackTracker > bt
Definition: HoughTrack_module.cc:146
htk::HoughTrack::houghIntercept
std::vector< double > houghIntercept
Definition: HoughTrack_module.cc:289
htk::HoughTrack::hough_peak_cut_
unsigned int hough_peak_cut_
Definition: HoughTrack_module.cc:133
htk::HoughTrack::houghView
std::vector< double > houghView
Definition: HoughTrack_module.cc:284
Geo.h
Collect Geo headers and supply basic geometry functions.
htk::HoughTrack::particleDirY
std::vector< double > particleDirY
Definition: HoughTrack_module.cc:238
htk::HoughTrack::houghSlope
std::vector< double > houghSlope
Definition: HoughTrack_module.cc:288
htk::HoughTrack::hough_label_
std::string hough_label_
Definition: HoughTrack_module.cc:131
calib::j
const double j
Definition: BetheBloch.cxx:29
htk::sortByNCell
bool sortByNCell(Track2D &track_i, Track2D &track_j)
Definition: HoughTrack_module.cc:2325
htk::Track_Match::Track_Match
Track_Match(std::vector< Track2D >::const_iterator x, std::vector< Track2D >::const_iterator y)
Definition: HoughTrack_module.cc:98
htk::HoughTrack::thetaYZ
std::vector< double > thetaYZ
Definition: HoughTrack_module.cc:208
htk::HoughTrack::particleT
std::vector< double > particleT
Definition: HoughTrack_module.cc:241
art::Event::event
EventNumber_t event() const
Definition: Event.h:67
htk::HoughTrack::beginJob
void beginJob() override
Definition: HoughTrack_module.cc:323
htk::HoughTrack::timestamp
unsigned long long timestamp
Definition: HoughTrack_module.cc:161
htk::HoughTrack::len
std::vector< double > len
Definition: HoughTrack_module.cc:194
htk::HoughTrack::subrun
unsigned long subrun
Definition: HoughTrack_module.cc:159
View
Definition: View.py:1
htk::HoughTrack::analyze
void analyze(art::Event const &e) override
Definition: HoughTrack_module.cc:471
run
Definition: run.py:1
om::cout
OStream cout
Definition: OStream.cxx:6
Event.h
htk::HoughTrack::particleDirZ
std::vector< double > particleDirZ
Definition: HoughTrack_module.cc:239
htk::HoughTrack::houghHeight
std::vector< double > houghHeight
Definition: HoughTrack_module.cc:287
art::TFileDirectory::make
T * make(ARGS...args) const
Definition: TFileDirectory.h:82
HoughResult.h
Data resulting from a Hough transform on the cell hit positions.
htk::HoughTrack::nplanesXZ
std::vector< double > nplanesXZ
Definition: HoughTrack_module.cc:187
htk::HoughTrack::startY
std::vector< double > startY
Definition: HoughTrack_module.cc:199
htk::HoughTrack::particlePhi
std::vector< double > particlePhi
Definition: HoughTrack_module.cc:243
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
ana::Sqrt
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
rawdata::RawTrigger::fTriggerTimingMarker_TimeStart
unsigned long long fTriggerTimingMarker_TimeStart
Definition: RawTrigger.h:38
hit
Definition: structs.h:12
htk::HoughTrack::thetaXZ
std::vector< double > thetaXZ
Definition: HoughTrack_module.cc:207
htk::HoughTrack::particleStartZ
std::vector< double > particleStartZ
Definition: HoughTrack_module.cc:233
htk::HoughTrack::phi
std::vector< double > phi
Definition: HoughTrack_module.cc:210
geo::GeometryBase::GetPlanesByView
const std::set< unsigned int > & GetPlanesByView(View_t v=kXorY) const
Definition: GeometryBase.cxx:918
simb::MCParticle::Momentum
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
cheat::BackTracker::HaveTruthInfo
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
kX
Definition: geo.h:1
htk::HoughTrack::particlePy
std::vector< double > particlePy
Definition: HoughTrack_module.cc:228
htk::HoughTrack::HoughTrack
HoughTrack(fhicl::ParameterSet const &p)
Definition: HoughTrack_module.cc:303
htk::HoughTrack::theta
std::vector< double > theta
Definition: HoughTrack_module.cc:209
htk::HoughTrack::zenith
std::vector< double > zenith
Definition: HoughTrack_module.cc:213
art::Event::time
Timestamp time() const
Definition: Event.h:61
htk::HoughTrack::nhits
std::vector< double > nhits
Definition: HoughTrack_module.cc:186
rb::Cluster::TotalGeV
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
htk::HoughTrack::particleLen
std::vector< double > particleLen
Definition: HoughTrack_module.cc:240
htk::HoughTrack::cosPhi
std::vector< double > cosPhi
Definition: HoughTrack_module.cc:212
htk::HoughTrack::total3DTracks
unsigned int total3DTracks
Definition: HoughTrack_module.cc:153
art::FindManyP::get
size_type get(size_type i, reference item, data_reference data) const
Definition: FindManyP.h:469
htk::HoughTrack::geometry_
art::ServiceHandle< geo::Geometry > geometry_
Definition: HoughTrack_module.cc:149
e
Float_t e
Definition: plot.C:35
art::Event::run
RunNumber_t run() const
Definition: Event.h:77
htk::HoughTrack::max_distance_between_hits_
double max_distance_between_hits_
Definition: HoughTrack_module.cc:135
htk::HoughTrack::particleStartX
std::vector< double > particleStartX
Definition: HoughTrack_module.cc:231
htk::HoughTrack::dirY
std::vector< double > dirY
Definition: HoughTrack_module.cc:205
htk::HoughTrack::startPlaneYZ
std::vector< double > startPlaneYZ
Definition: HoughTrack_module.cc:192
htk::HoughTrack::particlePDG
std::vector< double > particlePDG
Definition: HoughTrack_module.cc:222
htk::HoughTrack::houghTheta
std::vector< double > houghTheta
Definition: HoughTrack_module.cc:286
htk::HoughTrack::endJob
void endJob() override
Definition: HoughTrack_module.cc:2348
Geometry.h
Encapsulate the geometry of one entire detector (near, far, ndos)
htk::HoughTrack::azimuth
std::vector< double > azimuth
Definition: HoughTrack_module.cc:214
htk::HoughTrack::houghRho
std::vector< double > houghRho
Definition: HoughTrack_module.cc:285
rb::SortByPlane
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).
Definition: CellHit.cxx:124
art::Handle::failedToGet
bool failedToGet() const
Definition: Handle.h:196
htk::Track_Match::y_track
std::vector< Track2D >::const_iterator y_track
Definition: HoughTrack_module.cc:101
htk::HoughTrack::particleCosTheta
std::vector< double > particleCosTheta
Definition: HoughTrack_module.cc:244
htk::HoughTrack::endPlaneXZ
std::vector< double > endPlaneXZ
Definition: HoughTrack_module.cc:191
htk::HoughTrack::particlePz
std::vector< double > particlePz
Definition: HoughTrack_module.cc:229
htk::HoughTrack::dirX
std::vector< double > dirX
Definition: HoughTrack_module.cc:204
msf_helper.timespec
timespec
Definition: msf_helper.py:31
string
enum BeamMode string
Generated on Fri Jan 15 2021 14:12:24 for NOvA by   doxygen 1.8.11