AirKalmanAna_module.cc
Go to the documentation of this file.
1 
2 ////////////////////////////////////////////////////////////////////////
3 // Class: AirShower
4 // Module Type: analyzer
5 // File: AirKalmanAna_module.cc
6 //
7 // Updated at February 17, 2015 by Mehreen Sultana using artmod
8 // from cetpkgsupport v1_04_02.
9 ////////////////////////////////////////////////////////////////////////
10 
11 // C++ includes
12 
13 
14 #include <boost/format.hpp>
15 #include <math.h>
16 #include <string>
17 #include <vector>
18 #include <memory>
19 #include <fstream>
20 #include <iostream>
21 
22 
23 // ROOT includes
24 #include "TH1D.h"
25 #include "TH2D.h"
26 #include "TH1F.h"
27 #include "TH2F.h"
28 #include "TH3D.h"
29 #include "TVector3.h"
30 #include "TVector2.h"
31 #include "TTree.h"
32 #include "TTimeStamp.h"
33 
34 // Framework includes
46 #include "fhiclcpp/ParameterSet.h"
47 
48 
49 // NOvA includes
50 #include "Calibrator/Calibrator.h"
51 #include "Geometry/Geometry.h"
53 #include "RecoBase/CellHit.h"
54 #include "RecoBase/Cluster.h"
55 #include "RecoBase/Track.h"
56 #include "RecoBase/RecoHit.h"
57 #include "VertexFinder/DOCAInfo.h"
58 #include "Geometry/Geometry.h"
59 #include "Geometry/LiveGeometry.h"
60 #include "GeometryObjects/Geo.h"
61 #include "MCCheater/BackTracker.h"
62 #include "Utilities/AssociationUtil.h"
63 #include "Utilities/func/MathUtil.h"
64 #include "NovaDAQConventions/DAQConventions.h"
65 #include "RecoBase/FilterList.h"
67 
68 namespace air {
69  class AirKalmanAna;
70 }
71 
73 public:
74  explicit AirKalmanAna(fhicl::ParameterSet const & p);
75  virtual ~AirKalmanAna();
76 
77  void beginJob();
78  void analyze(art::Event const & e) override;
79 
80  virtual void beginSubRun(const art::SubRun & sr);
81 
82  void endJob();
83 private:
85 
86  void OppSignTrkCount(art::Ptr<rb::Track> &ftrack, std::vector<art::Ptr<rb::Track> > &tracks, unsigned int &oppsigntrks, unsigned int &samesigntrks, unsigned int &i);
87 
88  bool CheckContainer(std::vector<art::Ptr<rb::Track> > &trkcollection, art::Ptr<rb::Track> &ftrack) const;
89 
90  bool DoTracksIntersectYZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2) const;
91 
92  bool DoTracksIntersectXZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2) const;
93 
94 
95 
96  std::map<std::string, TH1*> h_;
97 
98  TTree* TrackTree;
99  TTree* SliceTree;
100  TTree* EventTree;
101  TTree* AnaTree;
102  TTree* ResultsTree;
103 
104  void z2north (std::vector<TVector3> track_i,
105  std::vector<TVector3> track_o);
106 
107  void Unitize (art::Ptr<rb::Track> &track, TVector3 &unitvec);
108  TVector3 unitvec;
109 
110  void z2north (TVector3 &unitvec, TVector3 &unitvec_newxz);
111  TVector3 unitvec_newxz; // Rotated about zenith by the clockwise orientation of the detector relative to the true north (Phi)
112 
113  void AvgTrkVec (std::vector<TVector3> &trkcol, TVector3 &avgtrk);
114  std::vector<TVector3> trkcol; // any collection of track unit vectors (Start - End)
115  TVector3 avgtrk; // Averaging all the unit track vectors
116  TVector3 avgrottrk;
117  TVector3 avgunittrk;
118 
119  void StdDev (std::vector<double> &anglecol, double &stddev);
120  std::vector<double> anglecol;
121  double stddev;
122 
123  void GetThetaXYZ (art::Ptr<rb::Track> &track, double &thetaX, double &thetaY, double &thetaZ);
124 
125  void GetThetaXYZ (TVector3 &track, double &thetaX, double &thetaY, double &thetaZ);
126 
127  void AvgValue(std::vector<double> &pramcol, double &avgpram);
128  std::vector<double> pramcol;
129  double avgpram;
130 
131  void AngleBetween (TVector3 &avgtrk, TVector3 &trkvec, double &sigma);
132  TVector3 trkvec;
133 
134  void AngleBetweentrk (std::vector<TVector3> &trkcol, TVector3 &trkvec, double &sigma, std::vector<double> &sigmacol);
135 
136  void CheckAssociations(art::Ptr<rb::Cluster> slice,std::vector<art::Ptr<rb::Track> > tracks, bool filtered);
137 
138 
139 
140 
141 
142  /////////////////////////////////// DEFINING HISTOS ////////////////////////////////////////
143 
144  //Individual Track Level Histograms
145 
146 
147  // TH1D* fthetaY; // Theta Y of each track in degrees
148  // TH1D* fthetaX; // Theta X of each track in degrees
149  // TH1D* fthetaZ; // Theta Z of each track in degrees
150 
151  // TH1D* fthetaXX;
152  // TH1D* fthetaYY;
153  // TH1D* fthetaZZ;
154 
155  // TH1D* fStdDev;
156  // TH1D* fStdDev2;
157 
158 
159 
160 
161  TH1D* fDOCA;
162 
164 
165  TH1D* fAllDeltaTheta; // Angle between average track and original individual track
167 
171 
172  TH1D* fDeltaAngle;
173 
174  TH1D* tottracks;
175  TH1D* conttracks;
176  TH1D* inttracks;
177  TH1D* opptracks;
178  TH1D* lentracks;
179 
180  TH1D* fAzimuth;
181  TH1D* fAltitude;
182 
183 
184  TH1D* fTracksPerSlice; // Total reconstructed tracks within a slice
185  TH1D* fTracksPerEvent; // Total reconstructed tracks within an event
186  // Total containment test passed tracks in slice
187 
189 
190  TH1D* fNSlices;
191 
192  TH1D* fStdDev;
193 
194  TH1D* fStdDev2;
195 
197 
199 
201 
203 
204  TH1D* fExposure;
205 
206  TH1D* fCosmicRate;
207 
208 
209  //////////////////////////////// DEFINING VARIABLES ///////////////////
210 
211  unsigned HighEnergyThresh; // High Energy ADC Threshold defined in the trigger
212 
213  unsigned int fNevents; //< total number of events counter
214 
215  unsigned ntracksinslice; // all tracks reconstructed in slice
216  unsigned ntracksinevent; // all tracks reconstructed in an event
217 
218  unsigned n3Dtracksinslice; // only 3D tracks reconstructed in slice
219 
220  unsigned n3Dtracksinevent; // Total 3D tracks in an event
221 
222  unsigned n3Dtrackspassedinevent; // Total 3D tracks passing containment cuts in event
223 
224  unsigned n3Dtrackspassedinslice; // Total 3D tracks passing containment cuts in slice
225 
226  unsigned nEmptyEvents; // To keep a count of empty events found in a set of data
227 
228  unsigned nslices; // Number of slices in the event with tracks that passed delta theta thing
229 
230  int fRun;
231 
232  int fSubRun;
233 
234  int fNTracks;
235 
236 
237 
238 
239  // Computational Variables;
240 
241  double cosY; // cosY = (end - start)/total_length // end and start of only Y components
242  double thetaY; // theta relative to the Y axis only
243  double thetaYY;
244  double cosX; // (end - start)/total_length // end and start of only X components
245  double thetaX; // theta relative to the X axis only
246  double thetaXX;
247  double cosZ; // (end - start)/total_length // end and start of only Z components
248  double thetaZ; // theta relative to the X axis only
249  double thetaZZ;
250 
251  double sigma; // Angle between individual tracks
252 
253 
254 
255 
256  double savgadc; // Average ADC per slice
257  double stotadc; // Total ADC per slice
258  unsigned sncell; // N Cells per Slice
259  double avgadc; // Average ADC per track
260  double totallength; // total length of track
261  double totaladc; // total adc of track
262  unsigned ncells; // N Cells per track
263 
264  // TRACK CONTAINTMENT: VECTORS
265 
266  TVector3 recoStart; // Starting points of all KalmanMerge reconstructed 3D tracks
267  TVector3 recoEnd; // Ending points of all KalmanMerge reconstructed 3D tracks
268 
269  TVector3 iXYZ;
270 
271  // TRACK CONTAINMENT: VARIABLES
272 
273  double detHalfWidth; // Detector width (X) center from the center
274  double detHalfHeight; // Detector length (Y) from the center
275  double detLength; // Dectector length (Z)
276 
277  double dist2edge;
278  double recostartX;
279  double recostartY;
280  double recostartZ;
281  double recoendX;
282  double recoendY;
283  double recoendZ;
284 
285  // Theta Histogram Analysis
286  double deltathetaY;
287  double deltathetaX;
288  double deltathetaZ;
289 
290  double thetaYavg;
291  double thetaXavg;
292  double thetaZavg;
293 
294 
295 
296  // DEFINE CONSTANTS
297  double phideg = 332.0 + (03.0/60.0) + (58.071769/3600); // Angle to rotate detector relative to true north
298  double phi = phideg*M_PI/180; // Above rotation Angle in radians
299  double SMALL_NUM = 0.0001; // Small number parameter for intersect checking
300 
301 
302  double fLiveTime;
303 
304 
305 
306  // Making an Event List so I don't have to keep typing it within the code and altering if loops
307 
308  std::vector<unsigned> Geventcol = {6930,6956, 6983, 6990, 7006, 7022, 7083, 7090, 7089, 7103,
309  7111, 7133, 7151, 7177, 7182, 7216, 7241, 7337, 7392,
310  7485, 7500, 7536, 7555, 7567, 7714, 7728, 7746, 7748,
311  7760, 7766, 7796, 7864, 7934, 7962, 7980, 7981, 8002,
312  8050, 8072, 8093, 8094, 8156, 8158, 8235, 8277, 8308,
313  8366, 8411, 8645, 8489, 8427, 8443, 8456, 8527, 8545,
314  8560, 8564, 8633, 8643, 8695, 8699, 8727, 8748, 8953,
315  9051, 9081, 9086, 9106, 9137, 9167, 9168};
316 
317  std::vector<unsigned> BadEventcol = {7392, 7702, 7672, 7839, 7629, 7587, 8149, 8065};
318 
319 
322 
323  unsigned eventpassed;
324  unsigned goldenevent;
325 
326 
327 
328 
329 
330  // GET GEOMETRY OPTIONS
331 
333 
335 
336 
337 
338  TTree* fTree;
339 
340  unsigned int fnAll;
341  unsigned int fnCont;
342  unsigned int fnInt;
343  unsigned int fnOpp;
344  unsigned int fnLen;
345 
346  TTree *subrun_tree_;
347  std::map<std::string, double> srt_;
348 
349  TTree *event_tree_;
350  std::map<std::string, double> et_;
351 
352  TTimeStamp startevent;
353 
354 };
355 
356 
357 
359  : EDAnalyzer(p)
360  , hit_label_(p.get<std::string>("hit_label"))
361  , slice_label_(p.get<std::string>("slice_label"))
362  , track_label_(p.get<std::string>("track_label"))
363  , fRun (-999)
364  , fSubRun (-999)
365  , fNTracks (0)
366  , fnAll (0)
367  , fnCont (0)
368  , fnInt (0)
369  , fnOpp (0)
370  , fnLen (0)
371 {
372 }
373 
374 
375 
377 {
378 }
379 
380 
381 
383 {
385 
386  subrun_tree_ = tfs->make<TTree>("SubRunTree", "SubRun Details");
387  std::vector<std::string> subrun_names {"exposure", "totevents", "tottracks", "passedtracks", "passedevents", "cosmicrate"};
388  for (auto const& name : subrun_names)
389  subrun_tree_->Branch(name.c_str(), &srt_[name], (name + "/D").c_str());
390 
391  event_tree_ = tfs->make<TTree>("EventTree", "Event Details");
392  std::vector<std::string> event_names {"time", "timestamp"};
393  for (auto const& name : event_names)
394  event_tree_->Branch(name.c_str(), &et_[name], (name + "/D").c_str());
395 
396  // Containment Tree
397 
398  fDOCA = tfs->make<TH1D>("DOCA",
399  "Closest distance between two tracks; Distance (m) ;N 3DTracks / 1 m",
400  100,-0.5, 99.5);
401 
402  fStartEndDist = tfs->make<TH1D>("StartEndDist",
403  "enddist - startdist (between 2 tracks) ; Distance (cm) ;N 3DTracks / 2.0 cm",
404  500,-1000.5, 999.5);
405 
406 
407  fAllDeltaTheta = tfs->make<TH1D>("All",
408  "#Delta #theta ;#Delta #theta (degrees) ;N 3DTracks / 0.01 degrees",
409  9000,-0.005, 89.995);
410 
411  fContDeltaTheta = tfs->make<TH1D>("Containment",
412  "#Delta #theta ;#Delta #theta (degrees) ;N 3DTracks / 0.01 degrees",
413  9000,-0.005, 89.995);
414 
415  fIntDeltaTheta = tfs->make<TH1D>("Intersect",
416  "#Delta #theta;#Delta #theta (degrees);N 3DTracks / 0.01 degrees",
417  9000,-0.005, 89.995);
418 
419 
420  fSlopeDeltaTheta = tfs->make<TH1D>("OppositeSlope",
421  "#Delta #theta;#Delta #theta (degrees) ;N 3DTracks / 0.01 degrees",
422  9000,-0.005, 89.995);
423 
424  fLenDeltaTheta = tfs->make<TH1D>("Length",
425  "#Delta #theta;#Delta #theta (degrees) ;N 3DTracks / 0.01 degrees",
426  9000,-0.005, 89.995);
427 
428  fAzimuth = tfs->make<TH1D>("Azimuth",
429  "Azimutal Angle;Azimuth (degrees); N Tracks / 1.0^{o}",
430  180,-89.5, 89.5);
431 
432  fAltitude = tfs->make<TH1D>("Altitude",
433  "Altitude Angle;Altitude (degrees);N Tracks / 1.0^{o}",
434  90, -0.50 , 89.5);
435 
436 
437  fDeltaAngle = tfs->make<TH1D>("DeltaAngle",
438  "(#Delta#theta)^{2} /2 ; (#Delta#theta)^{2} /2 (deg^{2}); N Muon Pairs/ 0.01 (deg^{2})",
439  9000, -0.5, 89.5);
440 
441  fAltMultiplicity = tfs->make<TH2D>("DeltaThetaAverage",
442  "#Delta #theta_{average} vs Tracks per Slice; #Delta #theta_{average} per slice (Degrees); N Tracks per slice",
443  90,-0.5, 89.5,
444  50, -0.5, 49.5);
445 
446 
447  fAzimuthMultiplicity = tfs->make<TH2D>("AzMulti",
448  "Azimuth_{average} vs Tracks per Slice; Azimuth_{average} per slice (Degrees); N Tracks per slice",
449  180,-89.5, 90.5,
450  50, -0.5, 49.5);
451 
452 
453  fSigmaAltitude = tfs->make<TH2D>("SigmaAlt",
454  "#sigma_{average} vs Altitude_{average} per slice; #sigma_{average} per slice (Degrees); Altitude_{average} per slice",
455  15,-0.5, 14.5,
456  90,-0.5, 89.5);
457 
458  // Event Level Histograms
459 
460  fTracksPerEvent = tfs->make<TH1D>("TracksPerEvent",
461  "Tracks Per Event;N 3D Tracks; NEvents",
462  201,-0.5,200.5);
463 
464  fTracksPerSlice = tfs->make<TH1D>("TracksPerSlice",
465  "Tracks Per Slice;N 3D Tracks; NSlices",
466  501,-0.5,500.5);
467 
468  fTracksPerSliceFinal = tfs->make<TH1D>("TracksPerSliceFinal",
469  "3D Tracks Passed Per Slice;N 3D Tracks; NSlices",
470  501,-0.5,500.5);
471 
472  fNSlices = tfs->make<TH1D>("NSlices",
473  "N Slices with PASSED Tracks per Event; N Slices per Event; N Events",
474  50, -0.5, 49.5);
475 
476 
477 
478  fExposure = tfs->make<TH1D>("exposure",";Exposure (s)", 1000, 0., 10);
479 
480 
481  fCosmicRate = tfs->make<TH1D>("cosmicrate",";Cosmic Rate (kHz)", 500, 0., 100);
482 
483 
484  tottracks = tfs->make<TH1D>("alltrkcount", "ntracks;adc/track", 100, -0.5, 999.5);
485  conttracks = tfs->make<TH1D>("conttrkcount", "ntracks;adc/track", 100, -0.5, 999.5);
486  inttracks = tfs->make<TH1D>("inttrkcount", "ntracks;adc/track", 100, -0.5, 999.5);
487  opptracks = tfs->make<TH1D>("opptrkcount", "ntracks;adc/track", 100, -0.5, 999.5);
488  lentracks = tfs->make<TH1D>("lentrkcount", "ntracks;adc/track", 100, -0.5, 999.5);
489 
490 }
491 
492 
493 
494 
496 {
497 
498  if(filtered && tracks.size() > 0){
499 
500  mf::LogVerbatim("KalmanTrackMerge")<<"Filtered slice has non-zero number of reconstructed tracks. This should never happen.";
501 
502  // Should never create tracks on slices that have been filtered
503  abort();
504 
505  }
506 
507  if(slice->IsNoise() && tracks.size() > 0){
508 
509  mf::LogVerbatim("KalmanTrackMerge")<<"Noise slice has non-zero number of reconstructed tracks. This should never happen.";
510 
511  // Should never create tracks on slices that have been filtered
512 
513  abort();
514 
515  }
516 
517  // check that reco'd track hits all belong to the slice hits
518 
519 
520  art::PtrVector<rb::CellHit> slicehits = slice->AllCells();
521 
522  for(unsigned int itrack = 0; itrack < tracks.size(); ++itrack){
523 
524  // get all of the track hits for this track
525 
526  art::PtrVector<rb::CellHit> trackhits = tracks[itrack]->AllCells();
527  // loop over the track hits
528  for(unsigned int itrhit = 0; itrhit < trackhits.size(); ++itrhit){
529 
530  bool foundhit = false;
531 
532  // loop over all the slice hits
533 
534  for(unsigned int islhit = 0; islhit < slicehits.size(); ++islhit){
535  if(trackhits[itrhit] == slicehits[islhit]){
536 
537  foundhit = true;
538 
539  break;
540 
541  }
542 
543  }
544 
545  // after looking at all the slice hits, did we find this track hit. If not, yell.
546 
547  if(!foundhit){
548 
549  mf::LogVerbatim("KalmanTrackAna")<<"Track made from hits on a slice other than the slice the track is associated to. This should never happen.";
550 
551  // Should never create tracks from hits on slices other than the slice its associated with
552  abort();
553  }
554  }
555  }
556 
557 }
558 
559 
561  // SETTING UP VECTORS CONTAINING START AND END POINTS OF EACH TRACK
562 
563  recoStart.SetXYZ(track->Start().X(),
564  track->Start().Y(),
565  track->Start().Z());
566 
567  recoEnd.SetXYZ(track->Stop().X(),
568  track->Stop().Y(),
569  track->Stop().Z());
570 
571  TVector3 highest, lowest;
572  if (recoStart.Y() > recoEnd.Y())
573  {
574  highest = recoStart;
575  lowest = recoEnd;
576  }
577  else
578  {
579  highest=recoEnd;
580  lowest=recoStart;
581  }
582 
583 
584  TVector3 length = highest - lowest;
585 
586  totallength = length.Mag();
587 
588 
589  // Subtracting End and start points of Track 1 and dividing by magnitude to get unit vector for track 1
590 
591  unitvec.SetXYZ(length.X()/totallength,
592  length.Y()/totallength,
593  length.Z()/totallength);
594 
595  return;
596 
597 }
598 
599 
600 // Rotating unit vectors of tracks to align with local true north. Inverse of original rotation vector does counter clock wise rotation to align z axis towards true north
602  unitvec_newxz.SetXYZ( (unitvec.X()*cos(phi) ) - (unitvec.Z()*sin(phi) ),
603  unitvec.Y(),
604  (unitvec.X()*sin(phi) ) + ( unitvec.Z()*cos(phi) ) );
605 
606  return;
607 }
608 
609 // How to get an average Trackvector
610 
611 void air::AirKalmanAna::AvgTrkVec (std::vector<TVector3> &trkcol, TVector3 &avgtrk){
612 
613  TVector3 sumtrk;
614 
615  for (unsigned int u = 0; u < trkcol.size(); u++){
616  sumtrk += trkcol[u];
617  }
618 
619  avgtrk.SetXYZ(sumtrk.X()/trkcol.size(),
620  sumtrk.Y()/trkcol.size(),
621  sumtrk.Z()/trkcol.size());
622 
623  return;
624 
625 }
626 
627 // How to get Standard Deviation of a collection
628 // void air::AirKalmanAna::StdDev (std::vector<double> &anglecol, double &stddev){
629 
630 // double sum = 0.0;
631 
632 // for (unsigned int u = 0; u < anglecol.size(); u++){
633 // sum += anglecol[u];
634 // }
635 
636 // double avgangle = sum/anglecol.size();
637 
638 // double diff;
639 
640 // for (unsigned int u = 0; u < anglecol.size(); u++){
641 
642 // diff += (avgangle - anglecol[u])*(avgangle - anglecol[u]) ;
643 
644 // }
645 
646 // double variance = (diff)/anglecol.size();
647 
648 // stddev = TMath::Sqrt(variance);
649 
650 // //std::cout << "stddev is " << stddev << std::endl;
651 
652 // return;
653 
654 // }
655 
656 // Getting Theta XYZ of any reco track
657 
659 
660  recoStart.SetXYZ(track->Start().X(),
661  track->Start().Y(),
662  track->Start().Z());
663 
664  recoEnd.SetXYZ(track->Stop().X(),
665  track->Stop().Y(),
666  track->Stop().Z());
667 
668  TVector3 highest, lowest;
669  if (recoStart.Y() > recoEnd.Y())
670  {
671  highest = recoStart;
672  lowest = recoEnd;
673  }
674  else
675  {
676  highest=recoEnd;
677  lowest=recoStart;
678  }
679 
680 
681  TVector3 length = highest - lowest;
682  totallength = length.Mag();
683 
684  cosY = TMath::Abs(length.Y()/totallength);
685  thetaY = (acos(cosY))*180/M_PI;
686 
687  cosX = TMath::Abs(length.X())/totallength;
688  thetaX = (acos(cosX))*180/M_PI;
689 
690  cosZ = TMath::Abs(length.Z()/totallength);
691  thetaZ = (acos(cosZ))*180/M_PI;
692 
693  return;
694 
695 }
696 
697 // Getting the Average of any parameter
698 
699 void air::AirKalmanAna::AvgValue(std::vector<double> &pramcol, double &avgpram){
700  double sum = 0.0;
701  for (unsigned int u = 0; u < pramcol.size(); u++){
702  sum += pramcol[u];
703  }
704  avgpram = sum/pramcol.size();
705  return;
706 }
707 
708 //angle between specified tracks
709 void air::AirKalmanAna::AngleBetween (TVector3 &ftrack, TVector3 &trkvec,
710  double &sigma) {
711  double scalarproduct;
712  scalarproduct = (avgtrk.X()*trkvec.X()) + (avgtrk.Y()*trkvec.Y()) + (avgtrk.Z()*trkvec.Z());
713  sigma = (acos(scalarproduct))*180/M_PI;
714  return;
715 
716 }
717 
718 // angle between a track and other tracks from that same collection
719 
720 void air::AirKalmanAna::AngleBetweentrk (std::vector<TVector3> &trkcol,TVector3 &trkvec, double &sigma, std::vector<double> &sigmacol) {
721 
722  double scalarproduct;
723  for (unsigned int n = 0; n < trkcol.size(); n++) {
724 
725  if (trkcol[n] == trkvec) continue;
726  scalarproduct = (trkcol[n].X()*trkvec.X()) + (trkcol[n].Y()*trkvec.Y()) + (trkcol[n].Z()*trkvec.Z());
727  sigma = (acos(scalarproduct))*180/M_PI;
728  sigmacol.push_back(sigma);
729  }
730  return;
731 }
732 
733 
734 bool air::AirKalmanAna::CheckContainer(std::vector<art::Ptr<rb::Track> > &trkcollection, art::Ptr<rb::Track> &ftrack) const
735 {
736 
737  unsigned int sametrk = 0;
738  for (unsigned int i = 0; i < trkcollection.size() ; ++i)
739 
740  {
741  art::Ptr<rb::Track> trackin = trkcollection[i];
742 
743  if (ftrack == trackin )
744  {
745  ++sametrk;
746  return false;
747  }
748 
749  }
750 
751  if (sametrk > 0 )
752  {
753  return false;
754  }
755 
756  return true;
757 
758 }
759 
760 bool air::AirKalmanAna::DoTracksIntersectXZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2) const
761 {
762  double uxz[2];
763  uxz[0] = dirvec1.X();
764  uxz[1] = dirvec1.Z();
765 
766  double vxz[2];
767  vxz[0] = dirvec2.X();
768  vxz[1] = dirvec2.Z();
769 
770  double wxz[2];
771  wxz[0] = highest1.X() - highest2.X();
772  wxz[1] = highest1.Z() - highest2.Z();
773 
774  double I0xz[2] = {0.};//, I1xz[2]; // Define array for two intersect points
775 
776  double Dxz = uxz[0]*vxz[1] - uxz[1]*vxz[0];
777 
778  double highxz1[2] = {highest1.X(), highest1.Z()};
779 
780  // the segments are skew and may intersect in a point
781 
782  // get the intersect parameter for S1
783  double sIxz = (vxz[0]*wxz[1] - vxz[1]*wxz[0]) / Dxz;
784 
785  // get the intersect parameter for S2
786  double tIxz = (uxz[0]*wxz[1] - uxz[1]*wxz[0]) / Dxz;
787 
788  if ( (sIxz < 0 || sIxz > 1) &&
789  (tIxz < 0 || tIxz > 1) ) // no intersect with S1
790  {
791  //std::cout << "Track " << i << " does not intersect at sIxz < 0 || sIxz > 1" << std::endl;
792 
793  //std::cout << "Track " << j << " does not intersect at tIxz < 0 || tIxz > 1" << std::endl;
794 
795  return false;
796  }
797  else if ( (sIxz >= 0 && sIxz <= 1) &&
798  (tIxz >= 0 && tIxz <= 1) ){
799 
800  I0xz[0] = highxz1[0] + sIxz * uxz[0];
801  I0xz[1] = highxz1[1] + sIxz * uxz[1];// compute S1 intersect point
802 
803  //I1xz[0] = highxz2[0] + tIxz * vxz[0];
804  //I1xz[1] = highxz2[1] + tIxz * vxz[1];// compute S2 intersect point
805 
806  if (TMath::Abs(I0xz[0]) < 800) {
807  std::cout << "There is an XZ intersect on track AT point (" << I0xz[0] << ", 0, " << I0xz[1] << ")." << std::endl;
808 
809 
810  // std::cout << "There is an XZ intersect on track AT point (" << I1xz[0] << ", 0, " << I1xz[1] << ")." << std::endl;
811 
812  return true;
813  }
814 
815  return false;
816  } // End of ELSE loop computing the Intersect point on ith and jth track
817 
818  return false;
819 }
820 
821 
822 bool air::AirKalmanAna::DoTracksIntersectYZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2) const
823 {
824  double uyz[2];
825  uyz[0] = dirvec1.Y();
826  uyz[1] = dirvec1.Z();
827 
828  double vyz[2];
829  vyz[0] = dirvec2.Y();
830  vyz[1] = dirvec2.Z();
831 
832  double wyz[2];
833  wyz[0] = highest1.Y() - highest2.Y();
834  wyz[1] = highest1.Z() - highest2.Z();
835 
836  double highyz1[2] = {highest1.Y(), highest1.Z()};
837 
838  double I0yz[2] = {0.};//, I1yz[2];
839 
840  double Dyz = uyz[0]*vyz[1] - uyz[1]*vyz[0];
841 
842  // the segments are skew and may intersect in a point
843 
844  // get the intersect parameter for S1
845  double sIyz = (vyz[0]*wyz[1] - vyz[1]*wyz[0]) / Dyz;
846 
847  // get the intersect parameter for S2
848  double tIyz = (uyz[0]*wyz[1] - uyz[1]*wyz[0]) / Dyz;
849 
850  if ( (sIyz < 0 || sIyz > 1) &&
851  (tIyz < 0 || tIyz > 1) ) // no intersect with S1 or S2
852  {
853  //std::cout << "Track " << i << " does not intersect at tIxz < 0 || tIxz > 1" << std::endl;
854 
855  //std::cout << "Track " << j << " does not intersect at tIxz < 0 || tIxz > 1" << std::endl;
856 
857  return false;
858  }
859  else if ( (sIyz >= 0 && sIyz <= 1) &&
860  (tIyz >= 0 && tIyz <= 1) ) {
861 
862  I0yz[0] = highyz1[0] + sIyz * uyz[0];
863  I0yz[1] = highyz1[1] + sIyz * uyz[1];// compute S1 intersect point
864 
865 
866  //I1yz[0] = highyz2[0] + tIyz * vyz[0];
867  //I1yz[1] = highyz2[1] + tIyz * vyz[1];// compute S2 intersect point
868 
869 
870  if (TMath::Abs(I0yz[0]) < 800 ) {
871  std::cout << "There is a YZ intersect on track AT point (0," << I0yz[0] << ", " << I0yz[1] << ")." << std::endl;
872 
873 
874  //std::cout << "There is a YZ intersect on track AT point (0," << I1yz[0] << ", " << I1yz[1] << ")." << std::endl;
875 
876 
877  return true;
878  }
879  // USE Track.ID() -> Make a boolean function
880  // Fine me Track ID ->
881  } // End of Else statement for checking intersects, computing i'th and j'th track intersect
882 
883  return false;
884 }
885 
886 
887 void air::AirKalmanAna::OppSignTrkCount(art::Ptr<rb::Track> &ftrack, std::vector<art::Ptr<rb::Track> > &tracks, unsigned int &oppsigntrks, unsigned int &samesigntrks, unsigned int &i)
888 {
889 
890  for (unsigned int j = i+1; j < (tracks.size() - 1); ++j)
891  {
892  art::Ptr<rb::Track> track = tracks[j];
893 
894  TVector3 recoStart1, recoStart2, recoEnd1, recoEnd2;
895 
896 
897  recoStart1.SetXYZ(ftrack->Start().X(),
898  ftrack->Start().Y(),
899  ftrack->Start().Z());
900 
901  recoEnd1.SetXYZ(ftrack->Stop().X(),
902  ftrack->Stop().Y(),
903  ftrack->Stop().Z());
904 
905 
906 
907  TVector3 highest1, lowest1;
908  if (recoStart1.Y() > recoEnd1.Y())
909  {
910  highest1 = recoStart1;
911  lowest1 = recoEnd1;
912  }
913  else
914  {
915  highest1=recoEnd1;
916  lowest1=recoStart1;
917  }
918 
919  recoStart2.SetXYZ(track->Start().X(),
920  track->Start().Y(),
921  track->Start().Z());
922 
923  recoEnd2.SetXYZ(track->Stop().X(),
924  track->Stop().Y(),
925  track->Stop().Z());
926 
927  TVector3 highest2, lowest2;
928  if (recoStart2.Y() > recoEnd2.Y())
929  {
930  highest2 = recoStart2;
931  lowest2 = recoEnd2;
932  }
933  else
934  {
935  highest2=recoEnd2;
936  lowest2=recoStart2;
937  }
938 
939  TVector3 dirvec1, dirvec2, diffvec;
940 
941  dirvec1 = highest1-lowest1;
942  dirvec2 = highest2 - lowest2;
943 
944  double slope1Y = (highest1.Y() - lowest1.Y())/(highest1.Z() - lowest1.Z());
945 
946  double slope2Y = (highest2.Y() - lowest2.Y())/(highest2.Z() - lowest2.Z());
947 
948 
949  double slope1X = (highest1.X() - lowest1.X())/(highest1.Z() - lowest1.Z());
950 
951  double slope2X = (highest2.X() - lowest2.X())/(highest2.Z() - lowest2.Z());
952 
953  //std::cout << "The (YZ) SLOPE OF TRACK " << i << " = " << slope1Y << " and Track " << j << " = " << slope2Y << std::endl;
954 
955  double signtestX = slope1X*slope2X;
956 
957  double signtestY = slope1Y*slope2Y;
958 
959 
960  if (signtestX < 0 || signtestY < 0)
961  {
962 
963  std::cout << "Track " << i << " and Track " << j << " have opposite slopes. " << std::endl;
964 
965  ++oppsigntrks;
966 
967  }
968  else {
969 
970  ++samesigntrks;
971 
972  }
973 
974  } // end of j track loop
975 
976  return;
977 
978 
979 }
980 
981 
983 {
984 
986  e.getByLabel(hit_label_, hits);
987 
988  // Get Slices
989 
991  e.getByLabel(slice_label_,slicecol);
992 
993 
994  // Make Slice List
995  art::PtrVector<rb::Cluster> slicelist;
996  for(unsigned int i = 0; i<slicecol->size();++i)
997  {
998  art::Ptr<rb::Cluster>slice(slicecol,i);
999  slicelist.push_back(slice);
1000  }
1001 
1002 
1003 
1004  // setup the associations between tracks and slices
1005  art::FindManyP<rb::Track> fmtrack(slicecol,e,track_label_);
1006 
1007  if (fRun < 0) fRun = e.run();
1008  if (fSubRun < 0) fSubRun = e.subRun();
1009 
1010  std::cout << "BEGIN ANALYZING EVENT " << e.event() << " NOW." << std::endl;
1011 
1012  art::Timestamp ts = e.time();
1013  // make it into a TTimeStamp
1014  TTimeStamp tts(ts.timeHigh(), ts.timeLow());
1015 
1016  TTimeStamp nowevent = tts;
1017 
1018  et_["time"] = tts.AsDouble();
1019  et_["timestamp"] = tts;
1020  event_tree_->Fill();
1021  // std::cout << "Approximate live time so far is " << nowevent.AsDouble() << " seconds." << std::endl;
1022 
1023  ++fNevents; // Count the number of events in the subrun
1024 
1025 
1026  // const time_t timeSec = ts.timeHigh();
1027  // struct tm* timeStruct = localtime(&timeSec);
1028 
1029  // unsigned short year = timeStruct->tm_year + 1900;
1030  // unsigned short month = timeStruct->tm_mon + 1;
1031  // unsigned short day = timeStruct->tm_mday;
1032  // unsigned short doy = timeStruct->tm_yday + 1;
1033  // unsigned short hour = timeStruct->tm_hour + 1;
1034  // unsigned short minute = timeStruct->tm_min;
1035  // unsigned short second = timeStruct->tm_sec;
1036 
1037 
1038  std::ofstream outfile("OutFile_CosmicExposure.txt", std::ios::app);
1039  outfile << " *************************************************************** " << std::endl;
1040 
1041 
1042  ///////////////////////////////////////////////////////////////////////
1043 
1044  // Defining Detector Dimensions
1045 
1048  detLength = geom->DetLength();
1049 
1050  ///////////////////////////////////////////////////////////////////////
1051  // std::cout<<"WIDTH :"
1052  // <<detHalfWidth
1053  // <<" , HEIGHT :"
1054  // <<detHalfHeight
1055  // <<" , LENGTH :"
1056  // <<detLength<<std::endl;
1057  ///////////////////////////////////////////////////////////////////////
1058  // std::cout << "phi is " << phideg << " degrees." << std::endl;
1059 
1060  /////////////////////////////////////////////////////////////////////////////////////////////////
1061 
1062 
1063  // To Identify if this event is a golden event or not
1064  // To use a predefined list of specific events of interest
1065  // (hard coded above in GoldenEvents vector for DDenergy data from FD run 16820 taken on August 16th, 2014)
1066  // file found at /nova/ana/users/mfrank/data <-- But No longer available.
1067 
1068  // unsigned GEvent = 1;
1069 
1070  // for (unsigned i=0; i < Geventcol.size(); ++i){
1071  // if (e.event() == Geventcol[i]){
1072  // GEvent = 0;
1073  // std::cout << "! ! ! THIS IS A GOLDEN EVENT ! ! !" << std::endl;
1074  // }
1075  // }
1076 
1077  /////////////////////////////////////////////////////////////////////////////////////////////////
1078 
1079  nslices = 0;
1080 
1081 
1082  /////////////////////////////////////////////////////////////////////////////////////////////////
1083 
1084 
1085 
1086  for (unsigned int islice = 0; islice < slicelist.size(); ++islice)
1087  {
1088 
1089  std::vector<art::Ptr<rb::Track> > tracks; // A container for all tracks in the slice
1090 
1091 
1092  if(fmtrack.isValid())
1093  {
1094  tracks = fmtrack.at(islice);
1095  }
1096 
1097  //bool filtered = rb::IsFiltered(e, slicelist[islice]);
1098 
1099  //CheckAssociations(slicelist[islice],tracks, filtered);
1100 
1101  std::vector<art::Ptr<rb::Track> > n3Dtracks;
1102  std::vector<art::Ptr<rb::Track> > passedtracks; // A container for all tracks that passed containment in the slice
1103 
1104  std::vector<art::Ptr<rb::Track> > nocrosstracks; // A container for tracks that do not cross in a slice
1105 
1106  std::vector<art::Ptr<rb::Track> > sameslopetracks; // container for tracks that have the same signed slope
1107 
1108  std::vector<art::Ptr<rb::Track > > minlengthtracks;
1109 
1110  double scalarproduct;
1111 
1112  double deltatheta;
1113 
1114  std::vector<TVector3> alltracksvec;
1115  std::vector<TVector3> unittracks;
1116  std::vector<TVector3> passedveccol;
1117  std::vector<TVector3> samesignveccol;
1118  std::vector<TVector3> nocrossveccol;
1119  std::vector<TVector3> minlengthveccol;
1120  std::vector<TVector3> sigmapassedveccol;
1121  std::vector<art::Ptr<rb::Track> > sigmapassedtracks;
1122 
1123 
1124  if(slicelist[islice]->IsNoise()) continue; // skipping over noise slices
1125 
1126  ++nslices;
1127  // ADD condition "|| GEvent == 1)" to only get events of interest }
1128 
1129  //if(slicelist[islice]->TotalADC() < 375000) continue;
1130 
1131  //if(GEvent == 1) continue;
1132 
1133 
1134 
1135  // std::vector<double> thetaYcol;
1136  // std::vector<double> thetaXcol;
1137  // std::vector<double> thetaZcol;
1138  // std::vector<double> thetaYYcol;
1139  // std::vector<double> thetaYYYcol;
1140 
1141 
1142 
1143  // ===================================================================
1144  // First Cut: Check if Tracks Are "Contained." Want tracks that aren't.
1145  // ===================================================================
1146 
1147 
1148  // std::cout << "........ Now Testing Containment of Tracks." << std::endl;
1149 
1150 
1151  for (unsigned int i = 0; i < tracks.size(); ++i)
1152 
1153  {
1154  art::Ptr<rb::Track> track = tracks[i];
1155  if (track->Is3D())
1156  {
1157  TVector3 unitvec;
1158  TVector3 unitvec_newxz;
1159 
1160  Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
1161  z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
1162 
1164 
1165  alltracksvec.push_back(trkvec);
1166  n3Dtracks.push_back(track);
1167  tottracks->Fill(track->TotalADC());
1168  }
1169 
1170  }
1171 
1172 
1173  fNTracks += n3Dtracks.size();
1174  if(n3Dtracks.size() < 2) continue;
1175 
1176  fnAll += n3Dtracks.size();
1177  fTracksPerSlice->Fill(n3Dtracks.size());
1178 
1179  for (unsigned int i=0; i< n3Dtracks.size(); ++i)
1180  {
1181  art::Ptr<rb::Track> track = n3Dtracks[i];
1182 
1183  TVector3 unitvec;
1184  TVector3 unitvec_newxz;
1185 
1186  Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
1187  z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
1188 
1190  trkcol = alltracksvec;
1191 
1192 
1193  std::vector<double> sigmacol;
1194 
1195  for (unsigned int j = i+1; j < (trkcol.size() - 1); j++) {
1196  scalarproduct = (trkcol[j].X()*trkvec.X()) + (trkcol[j].Y()*trkvec.Y()) + (trkcol[j].Z()*trkvec.Z());
1197  deltatheta = (acos(scalarproduct))*180/M_PI;
1198 
1199  //std::cout << "TRACK " << i << " and TRACK " << j << " have deltatheta = " << deltatheta << std::endl;
1200 
1201  fAllDeltaTheta->Fill(deltatheta);
1202  } // end of for trkcol
1203 
1204  } // end of for n3D tracks
1205 
1206 
1207  for (unsigned int i=0;i<n3Dtracks.size();i++)
1208  {
1209  art::Ptr<rb::Track> track = n3Dtracks[i];
1210 
1211 
1212  if (track->Is3D() )
1213  {
1214  ///////////////////////////////////////////////////////////////////
1215  /// FROM TrackFit/CosmicTrackAna_module.cc
1216 
1217 
1218 
1219  // Get the end point of the track in a double array
1220  double xyz[3] = { track->Stop().X(), track->Stop().Y(), track->Stop().Z() };
1221  double xyzstart[3] = {track->Start().X(), track->Start().Y(), track->Start().Z() };
1222 
1223 
1224  double highest1[3], lowest1[3];
1225  if (xyzstart[1] > xyz[1])
1226  {
1227  highest1[0] = xyzstart[0];
1228  highest1[1] = xyzstart[1];
1229  highest1[2] = xyzstart[2];
1230  lowest1[0] = xyz[0];
1231  lowest1[1] = xyz[1];
1232  lowest1[2] = xyz[2];
1233  }
1234  else
1235  {
1236  highest1[0]=xyz[0];
1237  highest1[1]=xyz[1];
1238  highest1[2]=xyz[2];
1239  lowest1[0]=xyzstart[0];
1240  lowest1[1]=xyzstart[1];
1241  lowest1[2]=xyzstart[2];
1242  }
1243 
1244 
1245  // OLD definition of distance to edge of detector (NOW people use LiveGeometry)
1246  double dist2edgeOLD = geo::DistToEdge(lowest1,
1247  detHalfWidth,
1248  detHalfHeight,
1249  detLength);
1250 
1251 
1252  double dist2edgeStart = geo::DistToEdge(highest1,
1253  detHalfWidth,
1254  detHalfHeight,
1255  detLength);
1256 
1257  //////// LENGTH //////////////////////
1258 
1259  recoStart.SetXYZ(highest1[0],
1260  highest1[1],
1261  highest1[2]);
1262 
1263  recoEnd.SetXYZ(lowest1[0],
1264  lowest1[1],
1265  lowest1[2]);
1266 
1267 
1268  TVector3 length = recoEnd - recoStart;
1269  totallength = length.Mag()/100;
1270 
1271  //if (totallength < 5.0 ) continue; (might need to use length cut later. I don't want to be too restrictive earlier on)
1272 
1273  /////// CONTAINMENT CUTS ////////////////
1274 
1275  // Let's define a fiducial box
1276  // (i.e. a box where start point of track must be)
1277  if (dist2edgeStart > 100.0){
1278  // Let's define a containment box
1279  // (i.e. a box where end point of track must be)
1280  if (dist2edgeOLD > 100.0){
1281  // function returns a value greater
1282  // than 25.0 cm inside detector in order to be contained
1283  continue;
1284  } // end of not fiducial and no contained loop -> Tracks to get rid of
1285  } // end of Fiducial if Loop
1286  else {
1287  passedtracks.push_back(track);
1288 
1289  } // End of Fiducial Else loop
1290  } // end Is3D loop
1291  } // end For (Tracks) loop
1292 
1293 
1294  // ===================== END CONTAINMENT LOOP ========================
1295 
1296  // Only looking at slices with more than 1 containment passed track
1297  if (passedtracks.size() < 2) continue;
1298  fnCont += passedtracks.size();
1299  for (unsigned int i=0; i<passedtracks.size(); ++i)
1300  {
1301 
1302  art::Ptr<rb::Track> track = passedtracks[i];
1303 
1304  TVector3 unitvec;
1305  TVector3 unitvec_newxz;
1306 
1307  Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
1308  z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
1309 
1311 
1312  passedveccol.push_back(trkvec);
1313 
1314  }
1315 
1316  /// MAKE THE FIRST DELTATHETA PLOT FOR CONTAINMENT PASSED TRACKS
1317 
1318  std::cout << "............Now Calculating DeltaTheta of PassedTracks." << std::endl;
1319 
1320  for (unsigned int i=0; i<passedtracks.size(); ++i)
1321  {
1322  art::Ptr<rb::Track> track = passedtracks[i];
1323 
1324  TVector3 unitvec;
1325  TVector3 unitvec_newxz;
1326 
1327  Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
1328  z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
1329 
1331  trkcol = passedveccol;
1332 
1333  std::vector<double> sigmacol;
1334 
1335  for (unsigned int j = i+1; j < (trkcol.size() - 1); j++) {
1336  scalarproduct = (trkcol[j].X()*trkvec.X()) + (trkcol[j].Y()*trkvec.Y()) + (trkcol[j].Z()*trkvec.Z());
1337  deltatheta = (acos(scalarproduct))*180/M_PI;
1338  fContDeltaTheta->Fill(deltatheta);
1339  }
1340  conttracks->Fill(track->TotalADC());
1341  }
1342 
1343 
1344  // ===================================================================
1345  // Second Cut : Check if Tracks Intersect in XZ or YZ View.
1346  // ===================================================================
1347 
1348  std::cout << "...... Now Testing Intersection between track " << std::endl;
1349 
1350  for (unsigned int i = 0; i < passedtracks.size(); ++i)
1351 
1352  {
1353  art::Ptr<rb::Track> track1 = passedtracks[i];
1354 
1355  unsigned int XZcrosscount = 0;
1356  unsigned int YZcrosscount = 0;
1357 
1358  unsigned int noXZcross = 0;
1359  unsigned int noYZcross = 0;
1360 
1361  for (unsigned int j = i + 1; j < (passedtracks.size() - 1); ++j)
1362  {
1363  art::Ptr<rb::Track> track2 = passedtracks[j];
1364 
1365  std::vector<art::Ptr<rb::Track> > trkcollection = passedtracks;
1366 
1367  TVector3 recoStart1, recoStart2, recoEnd1, recoEnd2;
1368 
1369  //trackcollection = passedtracks;
1370 
1371  recoStart1.SetXYZ(track1->Start().X(),
1372  track1->Start().Y(),
1373  track1->Start().Z());
1374 
1375  recoEnd1.SetXYZ(track1->Stop().X(),
1376  track1->Stop().Y(),
1377  track1->Stop().Z());
1378 
1379 
1380 
1381  TVector3 highest1, lowest1;
1382  if (recoStart1.Y() > recoEnd1.Y())
1383  {
1384  highest1 = recoStart1;
1385  lowest1 = recoEnd1;
1386  }
1387  else
1388  {
1389  highest1=recoEnd1;
1390  lowest1=recoStart1;
1391  }
1392 
1393  recoStart2.SetXYZ(track2->Start().X(),
1394  track2->Start().Y(),
1395  track2->Start().Z());
1396 
1397  recoEnd2.SetXYZ(track2->Stop().X(),
1398  track2->Stop().Y(),
1399  track2->Stop().Z());
1400 
1401  TVector3 highest2, lowest2;
1402  if (recoStart2.Y() > recoEnd2.Y())
1403  {
1404  highest2 = recoStart2;
1405  lowest2 = recoEnd2;
1406  }
1407  else
1408  {
1409  highest2=recoEnd2;
1410  lowest2=recoStart2;
1411  }
1412 
1413  TVector3 dirvec1, dirvec2, diffvec;
1414 
1415  dirvec1 = highest1-lowest1;
1416  dirvec2 = highest2 - lowest2;
1417 
1418 
1419  if (DoTracksIntersectYZ(highest1, lowest1, highest2, lowest2, dirvec1, dirvec2) )
1420  {
1421 
1422  //std::cout << "There is a YZ intersect on track " << i << " with track " << j << " in Slice " << nslices << std::endl;
1423 
1424  ++YZcrosscount;
1425  //continue;
1426  }
1427  else {
1428 
1429  ++noYZcross;
1430 
1431  }
1432 
1433 
1434  if (DoTracksIntersectXZ(highest1, lowest1, highest2, lowest2, dirvec1, dirvec2) ) {
1435 
1436  //std::cout << "There is a XZ intersect on track " << i << " with track " << j << " in Slice " << nslices << std::endl;
1437 
1438  ++XZcrosscount;
1439  //continue;
1440  }
1441  else {
1442 
1443  //std::cout << "There are no XZ intersects on track " << i << " with track " << j << " in Slice " << nslices << std::endl;
1444 
1445  ++noXZcross;
1446 
1447  }
1448 
1449  } // End of j'th track loop
1450 
1451  unsigned int sumcross = YZcrosscount + XZcrosscount;
1452 
1453  unsigned int sumnocross = noXZcross + noYZcross;
1454 
1455  if (sumcross < sumnocross || sumcross == 0)
1456 
1457  {
1458  //std::cout << "TRACK " << i << " DOES NOT INTERSECT WITH ANY TRACKS ON SLICE " << nslices << std::endl;
1459  nocrosstracks.push_back(track1);
1460  }
1461 
1462  if (sumcross > sumnocross)
1463 
1464  {
1465  std::cout << "TRACK " << i << " INTERSECTS " << sumcross << " TIMES IN XZ OR YZ VIEW. " << std::endl;
1466  continue;
1467 
1468  }
1469 
1470 
1471 
1472  } // End of i'th tracks loop
1473 
1474  if (nocrosstracks.size() < 2) continue;
1475 
1476  fnInt += nocrosstracks.size();
1477 
1478  for (unsigned int i = 0; i < nocrosstracks.size(); ++i)
1479 
1480  {
1481  art::Ptr<rb::Track> track = nocrosstracks[i];
1482 
1483 
1484  TVector3 unitvec;
1485  TVector3 unitvec_newxz;
1486 
1487  Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
1488  z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
1489 
1491 
1492  nocrossveccol.push_back(trkvec);
1493  inttracks->Fill(track->TotalADC());
1494  }
1495 
1496  std::cout << "............Now Calculating DeltaTheta of NoCrossTracks." << std::endl;
1497 
1498  //std::cout << "There are " << nocrosstracks.size() << "number of tracks in slice " << nslices << std::endl;
1499 
1500  for (unsigned int i=0; i<nocrosstracks.size(); ++i)
1501  {
1502  art::Ptr<rb::Track> track = nocrosstracks[i];
1503 
1504  TVector3 unitvec;
1505  TVector3 unitvec_newxz;
1506 
1507  Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
1508  z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
1509 
1511  trkcol = nocrossveccol;
1512 
1513  std::vector<double> sigmacol;
1514  double deltatheta;
1515 
1516  for (unsigned int j = i+1; j < (trkcol.size() - 1); j++) {
1517  scalarproduct = (trkcol[j].X()*trkvec.X()) + (trkcol[j].Y()*trkvec.Y()) + (trkcol[j].Z()*trkvec.Z());
1518  deltatheta = (acos(scalarproduct))*180/M_PI;
1519 
1520  fIntDeltaTheta->Fill(deltatheta);
1521  }
1522 
1523  }
1524 
1525  //=====================================================================
1526 
1527  // FINDING DISTANCE BETWEEN TRACKS
1528 
1529  // See Links below for reference/skeleton code used. Pretty much the same. Just trying to learn how to code for now. Will fix it/optimize it later if it can be optimized.
1530 
1531  // http://geomalgorithms.com/a05-_intersect-1.html
1532 
1533  // http://geomalgorithms.com/a07-_distance.html#dist3D_Segment_to_Segment
1534 
1535 
1536  //======================================================================
1537 
1538 
1539  std::vector<art::Ptr<rb::Track > > samesigncol;
1540 
1541 
1542  std::cout << ".........Now Checking Slope Signs." << std::endl;
1543 
1544  for (unsigned int i = 0; i < nocrosstracks.size(); ++i)
1545 
1546  {
1547  art::Ptr<rb::Track> track1 = nocrosstracks[i];
1548 
1549  art::Ptr<rb::Track> ftrack = track1;
1550 
1551  std::vector<art::Ptr<rb::Track > > tracks = nocrosstracks;
1552 
1553 
1554  std::vector<art::Ptr<rb::Track> > trkcollection = nocrosstracks;
1555 
1556  unsigned int oppsigntrks = 0;
1557  unsigned int samesigntrks = 0;
1558 
1559 
1560  OppSignTrkCount(ftrack, tracks, oppsigntrks, samesigntrks, i);
1561 
1562  if (oppsigntrks > samesigntrks) {
1563 
1564  std::cout << "Track " << i << " on slice " << nslices << "has opposite slope to " << oppsigntrks << " other tracks. " << std::endl;
1565  std::cout << "TRACK " << i << " IS A BAD TRACK." << std::endl;
1566  continue;
1567  }
1568 
1569  else {
1570 
1571  samesigncol.push_back(ftrack);
1572 
1573  }
1574 
1575 
1576  } // end of i track loop
1577 
1578  // if (opptracks > 0) continue;
1579 
1580  if (samesigncol.size() < 2) continue;
1581 
1582  fnOpp += samesigncol.size();
1583 
1584  for (unsigned int i = 0; i < samesigncol.size(); ++i)
1585 
1586  {
1587  art::Ptr<rb::Track> track = samesigncol[i];
1588 
1589 
1590  TVector3 unitvec;
1591  TVector3 unitvec_newxz;
1592 
1593  Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
1594  z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
1595 
1597 
1598  samesignveccol.push_back(trkvec);
1599  opptracks->Fill(track->TotalADC());
1600  }
1601 
1602  std::cout << "............Now Calculating DeltaTheta of SameSignTracks." << std::endl;
1603 
1604  for (unsigned int i=0; i<samesigncol.size(); ++i)
1605  {
1606  art::Ptr<rb::Track> track = samesigncol[i];
1607 
1608  TVector3 unitvec;
1609  TVector3 unitvec_newxz;
1610 
1611  Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
1612  z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
1613 
1615  trkcol = samesignveccol;
1616 
1617  for (unsigned int j = i+1; j < (trkcol.size() - 1); j++) {
1618  scalarproduct = (trkcol[j].X()*trkvec.X()) + (trkcol[j].Y()*trkvec.Y()) + (trkcol[j].Z()*trkvec.Z());
1619  deltatheta = (acos(scalarproduct))*180/M_PI;
1620  fSlopeDeltaTheta->Fill(deltatheta);
1621  }
1622 
1623  if (track->TotalLength() > 200.0) {
1624  minlengthtracks.push_back(track);
1625  }
1626  }
1627 
1628  if (minlengthtracks.size() < 2) continue;
1629  fnLen += minlengthtracks.size();
1630 
1631  for (unsigned int i = 0; i < minlengthtracks.size(); ++i)
1632 
1633  {
1634  art::Ptr<rb::Track> track = minlengthtracks[i];
1635 
1636 
1637  TVector3 unitvec;
1638  TVector3 unitvec_newxz;
1639 
1640  Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
1641  z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
1642 
1644 
1645  minlengthveccol.push_back(trkvec);
1646 
1647  }
1648 
1649  std::cout << "............Now Calculating DeltaTheta of MinLengthTracks." << std::endl;
1650 
1651 
1652  for (unsigned int i=0; i < minlengthtracks.size(); i++)
1653  {
1654  art::Ptr<rb::Track> track = minlengthtracks[i];
1655 
1656  TVector3 unitvec;
1657  TVector3 unitvec_newxz;
1658 
1659  Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
1660  z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
1661 
1663  trkcol = minlengthveccol;
1664 
1665  for (unsigned int j = i+1; j < (minlengthtracks.size() - 1); j++) {
1666  scalarproduct = (trkcol[j].X()*trkvec.X()) + (trkcol[j].Y()*trkvec.Y()) + (trkcol[j].Z()*trkvec.Z());
1667  deltatheta = (acos(scalarproduct))*180/M_PI;
1668  //std::cout << "TRACK " << i << " and TRACK " << j << " have deltatheta = " << deltatheta << std::endl;
1669 
1670  fLenDeltaTheta->Fill(deltatheta);
1671  }
1672  lentracks->Fill(track->TotalADC());
1673  }
1674 
1675  std::cout << ".......... Now Calculating Closest Distance Between Tracks" << std::endl;
1676 
1677  std::vector< double > DOCAmincol;
1678 
1679  for (unsigned int i = 0 ; i < minlengthtracks.size() ; ++i)
1680  {
1681  art::Ptr<rb::Track> track1 = minlengthtracks[i];
1682 
1683  std::vector<double> DOCAcol;
1684  for (unsigned int j = i + 1; j < (minlengthtracks.size() - 1); j++)
1685  {
1686  art::Ptr<rb::Track> track2 = minlengthtracks[j];
1687 
1688  TVector3 recoStart1, recoStart2, recoEnd1, recoEnd2;
1689 
1690  recoStart1.SetXYZ(track1->Start().X(),
1691  track1->Start().Y(),
1692  track1->Start().Z());
1693 
1694  recoEnd1.SetXYZ(track1->Stop().X(),
1695  track1->Stop().Y(),
1696  track1->Stop().Z());
1697 
1698  TVector3 highest1, lowest1;
1699  if (recoStart1.Y() > recoEnd1.Y())
1700  {
1701  highest1 = recoStart1;
1702  lowest1 = recoEnd1;
1703  }
1704  else
1705  {
1706  highest1=recoEnd1;
1707  lowest1=recoStart1;
1708  }
1709 
1710  recoStart2.SetXYZ(track2->Start().X(),
1711  track2->Start().Y(),
1712  track2->Start().Z());
1713 
1714  recoEnd2.SetXYZ(track2->Stop().X(),
1715  track2->Stop().Y(),
1716  track2->Stop().Z());
1717 
1718  TVector3 highest2, lowest2;
1719  if (recoStart2.Y() > recoEnd2.Y())
1720  {
1721  highest2 = recoStart2;
1722  lowest2 = recoEnd2;
1723  }
1724  else
1725  {
1726  highest2=recoEnd2;
1727  lowest2=recoStart2;
1728  }
1729 
1730  TVector3 dirvec1, dirvec2, diffvec;
1731 
1732  dirvec1 = highest1-lowest1;
1733  dirvec2 = highest2 - lowest2;
1734 
1735  diffvec = highest1 - highest2;
1736 
1737 
1738 
1739  double a = dirvec1.X()*dirvec1.X() + dirvec1.Y()*dirvec1.Y() + dirvec1.Z()*dirvec1.Z(); // Always >= 0
1740  double b = dirvec1.X()*dirvec2.X() + dirvec1.Y()*dirvec2.Y() + dirvec1.Z()*dirvec2.Z();
1741 
1742  double c = dirvec2.X()*dirvec2.X() + dirvec2.Y()*dirvec2.Y() + dirvec2.Z()*dirvec2.Z(); // Always >= 0
1743 
1744  double d = dirvec1.X()*diffvec.X() + dirvec1.Y()*diffvec.Y() + dirvec1.Z()*diffvec.Z();
1745 
1746  double e = dirvec2.X()*diffvec.X() + dirvec2.Y()*diffvec.Y() + dirvec2.Z()*diffvec.Z();
1747 
1748  double D = a*c - b*b; // Always >=0
1749 
1750  double sc, sN, sD = D; // "steps" sc = sN / sD, default sD = D >= 0
1751  double tc, tN, tD = D; // tc = tN / tD, default tD = D >= 0
1752 
1753 
1754  if (D < SMALL_NUM) { // The lines are almost parallel
1755 
1756  sN = 0.0; // force using point Highest1 on linesegment (track1)
1757 
1758  sD = 1.0; // to prevent possible division by 0.0 later
1759 
1760  tN = e;
1761 
1762  tD = c;
1763 
1764  }
1765  else { // get the closest points on the infinite lines
1766  sN = (b*e - c*d);
1767  tN = (a*e - b*d);
1768  if (sN < 0.0) { // sc < 0 => the s=0 edge is visible
1769  sN = 0.0;
1770  tN = e;
1771  tD = c;
1772  }
1773  else if (sN > sD) { // sc > 1 => the s=1 edge is visible
1774  sN = sD;
1775  tN = e + b;
1776  tD = c;
1777  }
1778  }
1779 
1780  if (tN < 0.0) { // tc < 0 => the t=0 edge is visible
1781  tN = 0.0;
1782  // recompute sc for this edge
1783  if (-d < 0.0)
1784  sN = 0.0;
1785  else if (-d > a)
1786  sN = sD;
1787  else {
1788  sN = -d;
1789  sD = a;
1790  }
1791  }
1792  else if (tN > tD) { // tc > 1 => the t=1 edge is visible
1793  tN = tD;
1794  // recompute sc for this edge
1795  if ((-d + b) < 0.0)
1796  sN = 0;
1797  else if ((-d + b) > a)
1798  sN = sD;
1799  else {
1800  sN = (-d + b);
1801  sD = a;
1802  }
1803  }
1804  // finally do the division to get sc and tc
1805  sc = (TMath::Abs(sN) < SMALL_NUM ? 0.0 : sN / sD);
1806  tc = (TMath::Abs(tN) < SMALL_NUM ? 0.0 : tN / tD);
1807 
1808  // get the difference of the two closest points
1809  TVector3 dP = diffvec + (sc * dirvec1 ) - (tc * dirvec2); // = S1(sc) - S2(tc)
1810 
1811  double closestdist = dP.Mag(); // return the closest distance
1812 
1813 
1814  DOCAcol.push_back(closestdist);
1815 
1816 
1817  // Here: Which track to add to the next list? How to loop over without double counting? Do I need to define a function here or can I use a ternary statement/condition?
1818 
1819  // I want to use "parallelness" as a test. Which means their slopes should not be to difference from each other. But how different is different?
1820 
1821  } // end for J tracks loop
1822 
1823  // Get minimum DOCA value.
1824  double DOCAmin = 10000;
1825  for (unsigned int k = 0; k < DOCAcol.size(); k++)
1826  {
1827  if (DOCAcol[i] < DOCAmin)
1828  DOCAmin = DOCAcol[i];
1829  }
1830 
1831  DOCAmincol.push_back(DOCAmin);
1832 
1833  double DOCAmeter = DOCAmin/100.0;
1834  fDOCA->Fill(DOCAmeter);
1835 
1836  } // end For (Tracks i ) loop
1837 
1838 
1839 
1840  /////////////////////////////////////////////////////////////////////
1841 
1842  // HOW TO GET ANGLES BETWEEN INDIVIDUAL TRACKS
1843 
1844  /////////////////////////////////////////////////////////////////////
1845 
1846 
1847 
1848  std::vector<double> azimuthcol;
1849  std::vector<double> altcol;
1850  std::vector<double> sigmacol;
1851 
1852  for (unsigned int i = 0; i < minlengthtracks.size(); i++){
1853 
1854  art::Ptr<rb::Track> track = minlengthtracks[i];
1855 
1856  TVector3 unitvec;
1857  TVector3 unitvec_newxz;
1858 
1859  Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
1860  z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
1861 
1863 
1864 
1865  trkcol = minlengthveccol;
1866 
1867 
1868  for (unsigned int j = i+1; j < (trkcol.size() - 1); j++) {
1869  scalarproduct = (trkcol[j].X()*trkvec.X()) + (trkcol[j].Y()*trkvec.Y()) + (trkcol[j].Z()*trkvec.Z());
1870  deltatheta = (acos(scalarproduct))*180/M_PI;
1871  double deltaangle = (deltatheta*deltatheta)/2;
1872  fDeltaAngle->Fill(deltaangle);
1873  }
1874 
1875  //////// LENGTH //////////////////////
1876 
1877  recoStart.SetXYZ(track->Start().X(),
1878  track->Start().Y(),
1879  track->Start().Z());
1880 
1881  recoEnd.SetXYZ(track->Stop().X(),
1882  track->Stop().Y(),
1883  track->Stop().Z());
1884 
1885  TVector3 highest, lowest;
1886  if (recoStart.Y() > recoEnd.Y())
1887  {
1888  highest = recoStart;
1889  lowest = recoEnd;
1890  }
1891  else
1892  {
1893  highest=recoEnd;
1894  lowest=recoStart;
1895  }
1896 
1897  TVector3 length = highest - lowest;
1898  totallength = length.Mag()/100;
1899 
1900 
1901  double alt = 90 - acos(trkvec.Y())*180/M_PI; // Altitude should just be complimentary to the zenith angle (90-zenith)
1902 
1903  // The following angle is based on the XZ projection of the track.
1904  // Not sure if this is the right method.
1905  // Since legnth of the unitvector is just 1, r = sin(Az).
1906  // Where r = magnitude of XZ projection
1907 
1908  double azimuth = asin(trkvec.Z()/TMath::Sqrt( (trkvec.X()*trkvec.X()) + (trkvec.Z()*trkvec.Z()) ) )*180/M_PI;
1909  azimuthcol.push_back(azimuth);
1910  altcol.push_back(alt);
1911 
1912  fAzimuth->Fill(azimuth);
1913  fAltitude->Fill(alt);
1914 
1915  }
1916 
1917 
1918  fTracksPerSliceFinal->Fill(minlengthtracks.size());
1919  pramcol = sigmacol;
1921 
1922  double avgsigma = avgpram;
1923 
1924  std::cout << "Average sigma per slice is " << avgsigma << std::endl;
1925 
1926  std::cout << "...............Finished Processing SLICE " << nslices << " with " << minlengthtracks.size() << " passed tracks in event " << e.event() << std::endl;
1927 
1928  } // End Slice loop
1929 
1930 
1931  fNSlices->Fill(nslices);
1932 
1933 
1934  std::cout << "******FINISHED PROCESSING EVENT NUMBER " << e.event() << " ******" << std::endl;
1935 
1936  std::cout << "************************************************************************" << std::endl;
1937 
1938  std::cout << "There are: " << fnAll << " total tracks in this job so far. " << std::endl;
1939 
1940  std::cout << "There are: " << fnCont << " containment passed tracks in this job so far. " << std::endl;
1941 
1942  std::cout << "There are: " << fnInt << " intersect passed tracks in this job so far. " << std::endl;
1943 
1944  std::cout << "There are: " << fnOpp << " slope passed tracks in this job so far. " << std::endl;
1945 
1946  std::cout << "There are: " << fnLen << " length passed tracks in this job so far. " << std::endl;
1947 
1948  std::cout << "There are: " << fNevents << " events in this job so far. " << std::endl;
1949 
1950 
1951  std::cout << "There are: " << fNTracks << " total tracks in this job so far. " << std::endl;
1952 
1953 
1954 
1955 
1956 
1957 } // End of Analyze Event Loop
1958 
1959 
1960 //......................................................................
1961 
1963 {
1965  sr.getByLabel("exposure", ce);
1966  if(ce.failedToGet()){
1967  mf::LogError("AirKalmanAna") << "No Cosmic Exposure Info";
1968  return;
1969  }
1970  fLiveTime += ce->totlivetime;
1971  LOG_INFO("AirKalmanAna") << "The live time is now " << fLiveTime;
1972  return;
1973 }
1974 
1975 //......................................................................
1976 
1977 
1978 
1979 
1980 //......................................................................
1981 
1983 
1984 {
1985 
1986  fExposure->Fill(fLiveTime);
1987 
1988 
1989  if(fLiveTime != 0){
1990 
1991  double cosmicrate = fNTracks/(1000.0*fLiveTime);
1992 
1993  srt_["cosmicrate"]=cosmicrate;
1994  srt_["exposure"]=fLiveTime;
1995  srt_["tottracks"]=fNTracks;
1996  srt_["passedevets"]=fNevents;
1997  srt_["passedtracks"]=fnLen;
1998 
1999 
2000 
2001  LOG_VERBATIM("AirShower") << "There are: " << fNTracks
2002 
2003  << " cosmic tracks selected in this job.\n"
2004 
2005  << "Livetime is " << fLiveTime
2006 
2007  << " according to CosmicExposure.\n"
2008 
2009  << "Note, livetime is for FULL SUBRUN, "
2010 
2011  << "so rate is wrong if haven't run full subrun.\n"
2012 
2013  << "The cosmic rate is therefore "
2014 
2015  << cosmicrate << " Hz \n";
2016 
2017 
2018 
2019 
2020  fCosmicRate->Fill(cosmicrate);
2021  subrun_tree_->Fill();
2022  }
2023  return;
2024 }
2025 
2026 
2027 
2028 
2029  // //==============================================================
2030  // //Code for Checking Intersects -> Don't need it since
2031  // //I defined the function.
2032  // //===============================================================
2033 
2034  // std::cout << "CHECKING FUNCTION AGAINST CODE" << std::endl;
2035 
2036  // double uxz[2];
2037  // uxz[0] = dirvec1.X();
2038  // uxz[1] = dirvec1.Z();
2039 
2040  // double vxz[2];
2041  // vxz[0] = dirvec2.X();
2042  // vxz[1] = dirvec2.Z();
2043 
2044  // double wxz[2];
2045  // wxz[0] = highest1.X() - highest2.X();
2046  // wxz[1] = highest1.Z() - highest2.Z();
2047 
2048  // double I0xz[2], I1xz[2]; // Define array for two intersect points
2049 
2050  // double Dxz = uxz[0]*vxz[1] - uxz[1]*vxz[0];
2051 
2052  // double highxz1[2], highxz2[2];
2053  // highxz1[0] = highest1.X();
2054  // highxz1[1] = highest1.Z();
2055  // highxz2[0] = highest2.X();
2056  // highxz2[1] = highest2.Z();
2057 
2058  // // the segments are skew and may intersect in a point
2059  // // get the intersect parameter for S1
2060  // double sIxz = (vxz[0]*wxz[1] - vxz[1]*wxz[0]) / Dxz;
2061  // if (sIxz < 0 || sIxz > 1) // no intersect with S1
2062  // {
2063  // //std::cout << "Track " << i << " does not intersect at sIxz < 0 || sIxz > 1" << std::endl;
2064  // }
2065 
2066  // // get the intersect parameter for S2
2067  // double tIxz = (uxz[0]*wxz[1] - uxz[1]*wxz[0]) / Dxz;
2068  // if (tIxz < 0 || tIxz > 1) // no intersect with S2
2069  // {
2070  // //std::cout << "Track " << j << " does not intersect at tIxz < 0 || tIxz > 1" << std::endl;
2071  // }
2072  // else {
2073 
2074  // I0xz[0] = highxz1[0] + sIxz * uxz[0];
2075  // I0xz[1] = highxz1[1] + sIxz * uxz[1];// compute S1 intersect point
2076 
2077  // I1xz[0] = highxz2[0] + tIxz * vxz[0];
2078  // I1xz[1] = highxz2[1] + tIxz * vxz[1];// compute S2 intersect point
2079 
2080  // std::cout << "There is an XZ intersect on track " << i << " with Track " << j << " in Slice " << nslices << " AT point (" << I0xz[0] << ", 0, " << I0xz[1] << ")." << std::endl;
2081 
2082 
2083  // std::cout << "There is an XZ intersect on track " << j << " with track " << i << " in Slice " << nslices << " AT point (" << I1xz[0] << ", 0, " << I1xz[1] << ")." << std::endl;
2084 
2085  // ++ntrackscrossXZ; // Count number of XZ intersecting tracks
2086  // ++ntrackscrossXZslice;
2087 
2088  // continue;
2089 
2090  // } // End of ELSE loop computing the Intersect point on ith and jth track
2091 
2092 
2093 
2094 
2095 
2096  // //===================================================================
2097 
2098  // // COMPUTE INTERSECT FOR YZ track segments
2099 
2100  // //===================================================================
2101  // double uyz[2];
2102  // uyz[0] = dirvec1.Y();
2103  // uyz[1] = dirvec1.Z();
2104 
2105  // double vyz[2];
2106  // vyz[0] = dirvec2.Y();
2107  // vyz[1] = dirvec2.Z();
2108 
2109  // double wyz[2];
2110  // wyz[0] = highest1.Y() - highest2.Y();
2111  // wyz[1] = highest1.Z() - highest2.Z();
2112 
2113  // double highyz1[2], highyz2[2];
2114  // highyz1[0] = highest1.Y();
2115  // highyz1[1] = highest1.Z();
2116  // highyz2[0] = highest2.Y();
2117  // highyz2[1] = highest2.Z();
2118 
2119  // double I0yz[2], I1yz[2];
2120 
2121  // double Dyz = uyz[0]*vyz[1] - uyz[1]*vyz[0];
2122 
2123  // // the segments are skew and may intersect in a point
2124  // // get the intersect parameter for S1
2125  // double sIyz = (vyz[0]*wyz[1] - vyz[1]*wyz[0]) / Dyz;
2126  // if (sIyz < 0 || sIyz > 1) // no intersect with S1
2127  // {
2128  // //std::cout << "Track " << i << " does not intersect at tIxz < 0 || tIxz > 1" << std::endl;
2129  // }
2130 
2131  // // get the intersect parameter for S2
2132  // double tIyz = (uyz[0]*wyz[1] - uyz[1]*wyz[0]) / Dyz;
2133  // if (tIyz < 0 || tIyz > 1) // no intersect with S2
2134  // {
2135  // //std::cout << "Track " << j << " does not intersect at tIxz < 0 || tIxz > 1" << std::endl;
2136  // }
2137  // else {
2138 
2139  // I0yz[0] = highyz1[0] + sIyz * uyz[0];
2140  // I0yz[1] = highyz1[1] + sIyz * uyz[1];// compute S1 intersect point
2141 
2142  // std::cout << "There is a YZ intersect on track " << i << " with track " << j << " in Slice " << nslices << " AT point (0," << I0yz[0] << ", " << I0yz[1] << ")." << std::endl;
2143 
2144  // I1yz[0] = highyz2[0] + tIyz * vyz[0];
2145  // I1yz[1] = highyz2[1] + tIyz * vyz[1];// compute S2 intersect point
2146 
2147  // std::cout << "There is a YZ intersect on track " << j << " with track " << i << " AND in Slice " << nslices << " AT point (0," << I1yz[0] << ", " << I1yz[1] << ")." << std::endl;
2148 
2149  // ++ntrackscrossYZ; // count number of YZ intersecting tracks
2150  // ++ntrackscrossYZslice;
2151 
2152  // // USE Track.ID() -> Make a boolean function
2153  // // Fine me Track ID ->
2154  // } // End of Else statement for checking intersects, computing i'th and j'th track intersect
2155 
2156 
2157 
2158 
2159 
2160 
AirKalmanAna(fhicl::ParameterSet const &p)
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
const XML_Char * name
Definition: expat.h:151
void Unitize(art::Ptr< rb::Track > &track, TVector3 &unitvec)
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
art::ServiceHandle< geo::LiveGeometry > livegeom
std::map< std::string, TH1 * > h_
constexpr std::uint32_t timeLow() const
Definition: Timestamp.h:30
virtual void beginSubRun(const art::SubRun &sr)
std::vector< double > anglecol
std::vector< double > pramcol
std::vector< unsigned > BadEventcol
double DistToEdge(double *point, double detHalfWidth, double detHalfHeight, double detLength)
Find the distance from the given point to the edge of the detector.
Definition: Geo.cxx:626
void AvgTrkVec(std::vector< TVector3 > &trkcol, TVector3 &avgtrk)
const char * p
Definition: xmltok.h:285
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:35
void AngleBetween(TVector3 &avgtrk, TVector3 &trkvec, double &sigma)
Definition: event.h:19
double DetLength() const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
void CheckAssociations(art::Ptr< rb::Cluster > slice, std::vector< art::Ptr< rb::Track > > tracks, bool filtered)
std::map< std::string, double > et_
T acos(T number)
Definition: d0nt_math.hpp:54
#define M_PI
Definition: SbMath.h:34
void analyze(art::Event const &e) override
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
Float_t Y
Definition: plot.C:38
void AvgValue(std::vector< double > &pramcol, double &avgpram)
Float_t Z
Definition: plot.C:38
bool DoTracksIntersectXZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2) const
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
double TotalADC() const
Sum of the ADC of all the contained hits.
Definition: Cluster.cxx:360
void hits()
Definition: readHits.C:15
length
Definition: demo0.py:21
void OppSignTrkCount(art::Ptr< rb::Track > &ftrack, std::vector< art::Ptr< rb::Track > > &tracks, unsigned int &oppsigntrks, unsigned int &samesigntrks, unsigned int &i)
std::map< std::string, double > srt_
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
const double a
void StdDev(std::vector< double > &anglecol, double &stddev)
Collect Geo headers and supply basic geometry functions.
Float_t d
Definition: plot.C:236
void AngleBetweentrk(std::vector< TVector3 > &trkcol, TVector3 &trkvec, double &sigma, std::vector< double > &sigmacol)
caf::StandardRecord * sr
const double j
Definition: BetheBloch.cxx:29
art::ServiceHandle< geo::Geometry > geom
EventNumber_t event() const
Definition: Event.h:67
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
void GetThetaXYZ(art::Ptr< rb::Track > &track, double &thetaX, double &thetaY, double &thetaZ)
size_type size() const
Definition: PtrVector.h:308
OStream cout
Definition: OStream.cxx:6
std::vector< TVector3 > trkcol
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
T sin(T number)
Definition: d0nt_math.hpp:132
double DetHalfWidth() const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
app
Definition: demo.py:189
bool DoTracksIntersectYZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2) const
void z2north(std::vector< TVector3 > track_i, std::vector< TVector3 > track_o)
const hit & b
Definition: hits.cxx:21
T cos(T number)
Definition: d0nt_math.hpp:78
Timestamp time() const
Definition: Event.h:61
unsigned n3Dtracksdeltathetapassedinevent
#define LOG_VERBATIM(category)
unsigned n3Dtracksdeltathetapassedinslice
#define LOG_INFO(stream)
Definition: Messenger.h:144
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
std::vector< unsigned > Geventcol
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
bool CheckContainer(std::vector< art::Ptr< rb::Track > > &trkcollection, art::Ptr< rb::Track > &ftrack) const
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
Definition: Prong.h:71
bool failedToGet() const
Definition: Handle.h:196
FILE * outfile
Definition: dump_event.C:13
T asin(T number)
Definition: d0nt_math.hpp:60