AirKalmanAna_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       AirShower
// Module Type: analyzer
// File:        AirKalmanAna_module.cc
//
// Updated at February 17, 2015  by Mehreen Sultana using artmod
// from cetpkgsupport v1_04_02.
////////////////////////////////////////////////////////////////////////

// C++ includes


#include <boost/format.hpp>
#include <math.h>
#include <string>
#include <vector>
#include <memory>
#include <fstream>
#include <iostream>


// ROOT includes
#include "TH1D.h"
#include "TH2D.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3D.h"
#include "TVector3.h"
#include "TVector2.h"
#include "TTree.h"
#include "TTimeStamp.h"

// Framework includes
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Provenance/Timestamp.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"


// NOvA includes
#include "Calibrator/Calibrator.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "RecoBase/RecoHit.h"
#include "VertexFinder/DOCAInfo.h"
#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "GeometryObjects/Geo.h"
#include "MCCheater/BackTracker.h"
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/MathUtil.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/FilterList.h"
#include "SummaryData/CosmicExposure.h"

namespace air {
  class AirKalmanAna;
}

class air::AirKalmanAna : public art::EDAnalyzer {
public:
  explicit AirKalmanAna(fhicl::ParameterSet const & p);
  virtual ~AirKalmanAna();

  void beginJob();
  void analyze(art::Event const & e) override;

  virtual void beginSubRun(const art::SubRun & sr);

  void endJob();
private:
  std::string hit_label_, slice_label_, track_label_;

  void OppSignTrkCount(art::Ptr<rb::Track> &ftrack, std::vector<art::Ptr<rb::Track> > &tracks, unsigned int &oppsigntrks, unsigned int &samesigntrks, unsigned int &i);

  bool CheckContainer(std::vector<art::Ptr<rb::Track> > &trkcollection, art::Ptr<rb::Track> &ftrack) const;

  bool DoTracksIntersectYZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2) const;

  bool DoTracksIntersectXZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2) const;

  
  
  std::map<std::string, TH1*> h_;
    
  TTree* TrackTree;
  TTree* SliceTree;
  TTree* EventTree;
  TTree* AnaTree;
  TTree* ResultsTree;

  void z2north (std::vector<TVector3> track_i,
                std::vector<TVector3> track_o);

  void Unitize (art::Ptr<rb::Track> &track, TVector3 &unitvec);
  TVector3 unitvec;

  void z2north (TVector3 &unitvec, TVector3 &unitvec_newxz);
  TVector3 unitvec_newxz; // Rotated about zenith by the clockwise orientation of the detector relative to the true north (Phi)

  void AvgTrkVec (std::vector<TVector3> &trkcol, TVector3 &avgtrk);
  std::vector<TVector3> trkcol; // any collection of track unit vectors (Start - End)
  TVector3 avgtrk; // Averaging all the unit track vectors
  TVector3 avgrottrk;
  TVector3 avgunittrk;

  void StdDev (std::vector<double> &anglecol, double &stddev);
  std::vector<double> anglecol;
  double stddev;

  void GetThetaXYZ (art::Ptr<rb::Track> &track, double &thetaX, double &thetaY, double &thetaZ);

  void GetThetaXYZ (TVector3 &track, double &thetaX, double &thetaY, double &thetaZ);

  void AvgValue(std::vector<double> &pramcol, double &avgpram);
  std::vector<double> pramcol;
  double avgpram;

  void AngleBetween (TVector3 &avgtrk, TVector3 &trkvec, double &sigma);
  TVector3 trkvec;

  void AngleBetweentrk (std::vector<TVector3> &trkcol, TVector3 &trkvec, double &sigma, std::vector<double> &sigmacol);

  void CheckAssociations(art::Ptr<rb::Cluster> slice,std::vector<art::Ptr<rb::Track> > tracks, bool filtered);

  



  /////////////////////////////////// DEFINING HISTOS ////////////////////////////////////////

  //Individual Track Level Histograms


  // TH1D* fthetaY;                      // Theta Y of each track in degrees
  // TH1D* fthetaX;                     // Theta X of each track in degrees
  // TH1D* fthetaZ;                     // Theta Z of each track in degrees

  // TH1D* fthetaXX;
  // TH1D* fthetaYY;
  // TH1D* fthetaZZ;

  // TH1D* fStdDev;
  // TH1D* fStdDev2;




  TH1D* fDOCA;

  TH1D* fStartEndDist; 
  
  TH1D* fAllDeltaTheta;  // Angle between average track and original individual track
  TH1D* fContDeltaTheta;

  TH1D* fIntDeltaTheta;
  TH1D* fSlopeDeltaTheta;
  TH1D* fLenDeltaTheta;

  TH1D* fDeltaAngle;

  TH1D* tottracks;
  TH1D* conttracks;
  TH1D* inttracks;
  TH1D* opptracks;
  TH1D* lentracks;

  TH1D* fAzimuth;
  TH1D* fAltitude;


  TH1D* fTracksPerSlice;             // Total reconstructed tracks within a slice
  TH1D* fTracksPerEvent;             // Total reconstructed tracks within an event
      // Total containment test passed tracks in slice

  TH1D* fTracksPerSliceFinal;  

  TH1D*   fNSlices;

  TH1D*   fStdDev;

  TH1D*   fStdDev2;

  TH1D*   fDeltaAnglePerSlice;

  TH2D*   fAltMultiplicity;

  TH2D*   fAzimuthMultiplicity;

  TH2D*   fSigmaAltitude;

  TH1D*   fExposure;    

  TH1D*   fCosmicRate;


  ////////////////////////////////  DEFINING VARIABLES ///////////////////

  unsigned        HighEnergyThresh;  // High Energy ADC Threshold defined in the trigger 
  
  unsigned int     fNevents;         //< total number of events counter 

  unsigned ntracksinslice;          // all tracks reconstructed in slice
  unsigned ntracksinevent;          // all tracks reconstructed in an event
  
  unsigned n3Dtracksinslice;        // only 3D tracks reconstructed in slice

  unsigned n3Dtracksinevent;        // Total 3D tracks in an event 

  unsigned n3Dtrackspassedinevent;  // Total 3D tracks passing containment cuts in event
  
  unsigned n3Dtrackspassedinslice;  // Total 3D tracks passing containment cuts in slice 

  unsigned nEmptyEvents;            // To keep a count of empty events found in a set of data
  
  unsigned nslices;                 // Number of slices in the event with tracks that passed delta theta thing

  int fRun;

  int fSubRun;

  int fNTracks; 




  // Computational Variables; 

  double cosY;                      // cosY = (end - start)/total_length // end and start of only Y components 
  double thetaY;                    // theta relative to the Y axis only
  double thetaYY;
  double cosX;                      // (end - start)/total_length // end and start of only X components 
  double thetaX;                    // theta relative to the X axis only
  double thetaXX;
  double cosZ;                      // (end - start)/total_length // end and start of only Z components 
  double thetaZ;                    // theta relative to the X axis only
  double thetaZZ;
  
  double sigma;                       // Angle between individual tracks


 
  
  double savgadc;                   // Average ADC per slice
  double stotadc;                   // Total ADC per slice
  unsigned sncell;                 // N Cells per Slice
  double avgadc;                    // Average ADC per track
  double totallength;               // total length of track
  double totaladc;                  // total adc of track
  unsigned ncells;                  // N Cells per track

  // TRACK CONTAINTMENT: VECTORS

  TVector3 recoStart;              // Starting points of all KalmanMerge reconstructed 3D tracks
  TVector3 recoEnd;                // Ending points of all KalmanMerge reconstructed 3D tracks

  TVector3 iXYZ;

  // TRACK CONTAINMENT: VARIABLES
  
  double detHalfWidth;              // Detector width (X) center from the center
  double detHalfHeight;             // Detector length (Y) from the center
  double detLength;                 // Dectector length (Z)

  double dist2edge;                 
  double recostartX;
  double recostartY;
  double recostartZ;
  double recoendX;
  double recoendY;
  double recoendZ; 
  
  // Theta Histogram Analysis 
  double deltathetaY;
  double deltathetaX;
  double deltathetaZ;

  double thetaYavg;
  double thetaXavg;
  double thetaZavg;
  


  // DEFINE CONSTANTS
  double phideg = 332.0 + (03.0/60.0) + (58.071769/3600); // Angle to rotate detector relative to true north
  double phi = phideg*M_PI/180; // Above rotation Angle in radians
  double SMALL_NUM = 0.0001; // Small number parameter for intersect checking


  double fLiveTime;



  // Making an Event List so I don't have to keep typing it within the code and altering if loops

  std::vector<unsigned> Geventcol = {6930,6956, 6983, 6990, 7006, 7022, 7083, 7090, 7089, 7103, 
                                     7111, 7133, 7151, 7177, 7182, 7216, 7241, 7337, 7392, 
                                     7485, 7500, 7536, 7555, 7567, 7714, 7728, 7746, 7748, 
                                     7760, 7766, 7796, 7864, 7934, 7962, 7980, 7981, 8002, 
                                     8050, 8072, 8093, 8094, 8156, 8158, 8235, 8277, 8308, 
                                     8366, 8411, 8645, 8489, 8427, 8443, 8456, 8527, 8545,
                                     8560, 8564, 8633, 8643, 8695, 8699, 8727, 8748, 8953,
                                     9051, 9081, 9086, 9106, 9137, 9167, 9168};

  std::vector<unsigned> BadEventcol = {7392, 7702, 7672, 7839, 7629, 7587, 8149, 8065};


  unsigned n3Dtracksdeltathetapassedinevent;
  unsigned n3Dtracksdeltathetapassedinslice;

  unsigned eventpassed;
  unsigned goldenevent;





  // GET GEOMETRY OPTIONS

  art::ServiceHandle<geo::Geometry> geom; //OLD

  art::ServiceHandle<geo::LiveGeometry> livegeom; // NEW



  TTree*                               fTree;

  unsigned int fnAll;
  unsigned int fnCont;
  unsigned int fnInt;
  unsigned int fnOpp;
  unsigned int fnLen;

  TTree *subrun_tree_;
  std::map<std::string, double> srt_;

  TTree *event_tree_;
  std::map<std::string, double> et_;

  TTimeStamp startevent;

};



air::AirKalmanAna::AirKalmanAna(fhicl::ParameterSet const & p) 
  : EDAnalyzer(p)
  , hit_label_(p.get<std::string>("hit_label"))
  , slice_label_(p.get<std::string>("slice_label"))
  , track_label_(p.get<std::string>("track_label"))
  , fRun (-999)
  , fSubRun (-999)
  , fNTracks (0)
  , fnAll (0)
  , fnCont (0)
  , fnInt (0)
  , fnOpp (0)
  , fnLen (0)
{
}



air::AirKalmanAna::~AirKalmanAna()
{
}



void air::AirKalmanAna::beginJob()
{
  art::ServiceHandle<art::TFileService> tfs;

  subrun_tree_ = tfs->make<TTree>("SubRunTree", "SubRun Details");
  std::vector<std::string> subrun_names {"exposure", "totevents", "tottracks", "passedtracks", "passedevents", "cosmicrate"};
  for (auto const& name : subrun_names)
    subrun_tree_->Branch(name.c_str(), &srt_[name], (name + "/D").c_str());
  
  event_tree_ = tfs->make<TTree>("EventTree", "Event Details");
  std::vector<std::string> event_names {"time", "timestamp"};
  for (auto const& name : event_names)
    event_tree_->Branch(name.c_str(), &et_[name], (name + "/D").c_str());
  
  // Containment Tree

  fDOCA                          = tfs->make<TH1D>("DOCA",
                                                   "Closest distance between two tracks; Distance (m) ;N 3DTracks / 1 m", 
                                                   100,-0.5, 99.5);

  fStartEndDist                  = tfs->make<TH1D>("StartEndDist",
                                                           "enddist - startdist (between 2 tracks) ; Distance (cm) ;N 3DTracks / 2.0 cm", 
                                                           500,-1000.5, 999.5);

  
  fAllDeltaTheta                          = tfs->make<TH1D>("All",
                                                         "#Delta #theta ;#Delta #theta (degrees) ;N 3DTracks / 0.01 degrees", 
                                                         9000,-0.005, 89.995);

  fContDeltaTheta                          = tfs->make<TH1D>("Containment",
                                                           "#Delta #theta ;#Delta #theta (degrees) ;N 3DTracks / 0.01 degrees", 
                                                           9000,-0.005, 89.995);

  fIntDeltaTheta                         = tfs->make<TH1D>("Intersect",
                                                          "#Delta #theta;#Delta #theta (degrees);N 3DTracks / 0.01 degrees", 
                                                          9000,-0.005, 89.995);


  fSlopeDeltaTheta                         = tfs->make<TH1D>("OppositeSlope",
                                                            "#Delta #theta;#Delta #theta (degrees) ;N 3DTracks / 0.01 degrees", 
                                                            9000,-0.005, 89.995);

  fLenDeltaTheta                         = tfs->make<TH1D>("Length",
                                                            "#Delta #theta;#Delta #theta (degrees) ;N 3DTracks / 0.01 degrees", 
                                                            9000,-0.005, 89.995);

  fAzimuth                         = tfs->make<TH1D>("Azimuth",
                                                     "Azimutal Angle;Azimuth (degrees); N Tracks / 1.0^{o}", 
                                                     180,-89.5, 89.5);

  fAltitude                         = tfs->make<TH1D>("Altitude",
                                                      "Altitude Angle;Altitude (degrees);N Tracks / 1.0^{o}", 
                                                      90, -0.50 , 89.5);


  fDeltaAngle                = tfs->make<TH1D>("DeltaAngle",
                                                       "(#Delta#theta)^{2} /2 ; (#Delta#theta)^{2} /2 (deg^{2});  N Muon Pairs/ 0.01 (deg^{2})",
                                                       9000, -0.5, 89.5);

  fAltMultiplicity                = tfs->make<TH2D>("DeltaThetaAverage",
                                                    "#Delta #theta_{average} vs Tracks per Slice; #Delta #theta_{average} per slice (Degrees); N Tracks per slice",
                                                    90,-0.5, 89.5,
                                                    50, -0.5, 49.5);
                                                    

  fAzimuthMultiplicity                = tfs->make<TH2D>("AzMulti",
                                                        "Azimuth_{average} vs Tracks per Slice; Azimuth_{average} per slice (Degrees); N Tracks per slice",
                                                        180,-89.5, 90.5,
                                                        50, -0.5, 49.5);


  fSigmaAltitude                      = tfs->make<TH2D>("SigmaAlt",
                                                        "#sigma_{average} vs Altitude_{average} per slice; #sigma_{average} per slice (Degrees); Altitude_{average} per slice",
                                                        15,-0.5, 14.5,
                                                        90,-0.5, 89.5);

  // Event Level Histograms

  fTracksPerEvent               = tfs->make<TH1D>("TracksPerEvent",
                                                  "Tracks Per Event;N 3D Tracks; NEvents",
                                                  201,-0.5,200.5); 

  fTracksPerSlice               = tfs->make<TH1D>("TracksPerSlice",
                                                  "Tracks Per Slice;N 3D Tracks; NSlices",
                                                  501,-0.5,500.5);

  fTracksPerSliceFinal           = tfs->make<TH1D>("TracksPerSliceFinal",
                                                   "3D Tracks Passed Per Slice;N 3D Tracks; NSlices",
                                                   501,-0.5,500.5);

  fNSlices                       = tfs->make<TH1D>("NSlices",
                                                   "N Slices with PASSED Tracks per Event; N Slices per Event; N Events",
                                                   50, -0.5, 49.5);



  fExposure   = tfs->make<TH1D>("exposure",";Exposure (s)", 1000, 0., 10);


  fCosmicRate = tfs->make<TH1D>("cosmicrate",";Cosmic Rate (kHz)", 500, 0., 100);


  tottracks = tfs->make<TH1D>("alltrkcount", "ntracks;adc/track", 100, -0.5, 999.5);
  conttracks = tfs->make<TH1D>("conttrkcount", "ntracks;adc/track", 100, -0.5, 999.5);
  inttracks = tfs->make<TH1D>("inttrkcount", "ntracks;adc/track", 100, -0.5, 999.5);
  opptracks = tfs->make<TH1D>("opptrkcount", "ntracks;adc/track", 100, -0.5, 999.5);
  lentracks = tfs->make<TH1D>("lentrkcount", "ntracks;adc/track", 100, -0.5, 999.5);

}




void air::AirKalmanAna::CheckAssociations(art::Ptr<rb::Cluster> slice,std::vector<art::Ptr<rb::Track> > tracks, bool filtered)  
{
        
  if(filtered && tracks.size() > 0){
        
    mf::LogVerbatim("KalmanTrackMerge")<<"Filtered slice has non-zero number of reconstructed tracks. This should never happen.";
        
    // Should never create tracks on slices that have been filtered     
    abort();
        
  }
        
  if(slice->IsNoise() && tracks.size() > 0){
        
    mf::LogVerbatim("KalmanTrackMerge")<<"Noise slice has non-zero number of reconstructed tracks. This should never happen.";
        
    // Should never create tracks on slices that have been filtered
        
    abort();
        
  }
                
  // check that reco'd track hits all belong to the slice hits
        

  art::PtrVector<rb::CellHit> slicehits = slice->AllCells();
        
  for(unsigned int itrack = 0; itrack < tracks.size(); ++itrack){
        
    // get all of the track hits for this track
        
    art::PtrVector<rb::CellHit> trackhits = tracks[itrack]->AllCells(); 
    // loop over the track hits 
    for(unsigned int itrhit = 0; itrhit < trackhits.size(); ++itrhit){
        
      bool foundhit = false;
        
      // loop over all the slice hits
        
      for(unsigned int islhit = 0; islhit < slicehits.size(); ++islhit){        
        if(trackhits[itrhit] == slicehits[islhit]){ 
        
          foundhit = true;
        
          break;
        
        }
        
      }
        
      // after looking at all the slice hits, did we find this track hit. If not, yell.
        
      if(!foundhit){
        
        mf::LogVerbatim("KalmanTrackAna")<<"Track made from hits on a slice other than the slice the track is associated to. This should never happen.";
        
        // Should never create tracks from hits on slices other than the slice its associated with
        abort();
      } 
    }   
  }
        
}


void air::AirKalmanAna::Unitize (art::Ptr<rb::Track> &track, TVector3 &unitvec) {
  // SETTING UP VECTORS CONTAINING START AND END POINTS OF EACH TRACK
          
  recoStart.SetXYZ(track->Start().X(),
                   track->Start().Y(),
                   track->Start().Z());
            
  recoEnd.SetXYZ(track->Stop().X(),
                 track->Stop().Y(),
                 track->Stop().Z());

  TVector3 highest, lowest;
  if (recoStart.Y() > recoEnd.Y())
  {
    highest = recoStart;
    lowest = recoEnd;
  }
  else
  {
    highest=recoEnd;
    lowest=recoStart;
  }


  TVector3 length = highest - lowest;

  totallength = length.Mag();

          
  // Subtracting End and start points of Track 1 and dividing by magnitude to get unit vector for track 1

  unitvec.SetXYZ(length.X()/totallength,
                 length.Y()/totallength,
                 length.Z()/totallength);

  return;

}


// 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
void air::AirKalmanAna::z2north (TVector3 &unitvec, TVector3 &unitvec_newxz){
  unitvec_newxz.SetXYZ( (unitvec.X()*cos(phi) ) - (unitvec.Z()*sin(phi) ), 
                        unitvec.Y(),
                        (unitvec.X()*sin(phi) ) + ( unitvec.Z()*cos(phi) ) );
  
  return;
}

// How to get an average Trackvector

void air::AirKalmanAna::AvgTrkVec (std::vector<TVector3> &trkcol, TVector3 &avgtrk){

  TVector3 sumtrk;

  for (unsigned int u = 0; u < trkcol.size(); u++){
    sumtrk += trkcol[u];
  }
  
  avgtrk.SetXYZ(sumtrk.X()/trkcol.size(),
                sumtrk.Y()/trkcol.size(),
                sumtrk.Z()/trkcol.size());

  return;
        
}

// How to get Standard Deviation of a collection
// void air::AirKalmanAna::StdDev (std::vector<double> &anglecol, double &stddev){
  
//   double sum = 0.0;

//   for (unsigned int u = 0; u < anglecol.size(); u++){
//     sum += anglecol[u];
//   }

//   double avgangle = sum/anglecol.size(); 
  
//   double diff;

//   for (unsigned int u = 0; u < anglecol.size(); u++){

//     diff += (avgangle - anglecol[u])*(avgangle - anglecol[u]) ;

//   }

//   double variance = (diff)/anglecol.size();
  
//   stddev = TMath::Sqrt(variance);

//   //std::cout << "stddev is " << stddev << std::endl;

//   return;

// }

// Getting Theta XYZ of any reco track

void air::AirKalmanAna::GetThetaXYZ (art::Ptr<rb::Track> &track, double &thetaX, double &thetaY, double &thetaZ){
          
  recoStart.SetXYZ(track->Start().X(),
                   track->Start().Y(),
                   track->Start().Z());
            
  recoEnd.SetXYZ(track->Stop().X(),
                 track->Stop().Y(),
                 track->Stop().Z());

  TVector3 highest, lowest;
  if (recoStart.Y() > recoEnd.Y())
  {
    highest = recoStart;
    lowest = recoEnd;
  }
  else
  {
    highest=recoEnd;
    lowest=recoStart;
  }


  TVector3 length = highest - lowest;
  totallength = length.Mag();

  cosY = TMath::Abs(length.Y()/totallength);  
  thetaY = (acos(cosY))*180/M_PI;

  cosX = TMath::Abs(length.X())/totallength; 
  thetaX = (acos(cosX))*180/M_PI;

  cosZ = TMath::Abs(length.Z()/totallength);
  thetaZ = (acos(cosZ))*180/M_PI;

  return;

}

// Getting the Average of any parameter 

void air::AirKalmanAna::AvgValue(std::vector<double> &pramcol, double &avgpram){
  double sum = 0.0;
  for (unsigned int u = 0; u < pramcol.size(); u++){
    sum += pramcol[u];
  }
  avgpram = sum/pramcol.size();
  return;
}

//angle between specified tracks
void air::AirKalmanAna::AngleBetween (TVector3 &ftrack, TVector3 &trkvec, 
                                   double &sigma) {
  double scalarproduct;
  scalarproduct = (avgtrk.X()*trkvec.X()) + (avgtrk.Y()*trkvec.Y()) + (avgtrk.Z()*trkvec.Z());
  sigma = (acos(scalarproduct))*180/M_PI;
  return;

}

// angle between a track and other tracks from that same collection

void air::AirKalmanAna::AngleBetweentrk (std::vector<TVector3> &trkcol,TVector3 &trkvec, double &sigma, std::vector<double> &sigmacol) {

  double scalarproduct;
  for (unsigned int n = 0; n < trkcol.size(); n++) {

    if (trkcol[n] == trkvec) continue;
    scalarproduct = (trkcol[n].X()*trkvec.X()) + (trkcol[n].Y()*trkvec.Y()) + (trkcol[n].Z()*trkvec.Z());
    sigma = (acos(scalarproduct))*180/M_PI;
    sigmacol.push_back(sigma);
  }
  return;
}


bool air::AirKalmanAna::CheckContainer(std::vector<art::Ptr<rb::Track> > &trkcollection, art::Ptr<rb::Track> &ftrack) const
{

  unsigned int sametrk = 0;
  for (unsigned int i = 0; i < trkcollection.size() ; ++i) 

  {
    art::Ptr<rb::Track> trackin = trkcollection[i];

    if (ftrack  == trackin )
    {
      ++sametrk;
      return false;
    }

  }

  if (sametrk > 0 ) 
  {
    return false;
  }

  return true;
    
}

bool air::AirKalmanAna::DoTracksIntersectXZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2) const
{
  double    uxz[2];
  uxz[0] = dirvec1.X();
  uxz[1] = dirvec1.Z();
        
  double    vxz[2];
  vxz[0] = dirvec2.X();
  vxz[1] = dirvec2.Z();
        
  double    wxz[2];
  wxz[0] = highest1.X() - highest2.X();
  wxz[1] = highest1.Z() - highest2.Z();
        
  double I0xz[2] = {0.};//, I1xz[2]; // Define array for two intersect points

  double     Dxz = uxz[0]*vxz[1] - uxz[1]*vxz[0];

  double highxz1[2] = {highest1.X(), highest1.Z()};
        
  // the segments are skew and may intersect in a point

  // get the intersect parameter for S1
  double     sIxz = (vxz[0]*wxz[1] - vxz[1]*wxz[0]) / Dxz;

  // get the intersect parameter for S2
  double     tIxz = (uxz[0]*wxz[1] - uxz[1]*wxz[0]) / Dxz;

  if ( (sIxz < 0 || sIxz > 1) &&
       (tIxz < 0 || tIxz > 1) )                // no intersect with S1
  {
    //std::cout << "Track " << i << " does not intersect at sIxz < 0 || sIxz > 1" << std::endl; 

    //std::cout << "Track " << j << " does not intersect at tIxz < 0 || tIxz > 1" << std::endl;

    return false;
  }
  else if ( (sIxz >= 0 && sIxz <= 1) &&
            (tIxz >= 0 && tIxz <= 1) ){
    
    I0xz[0] = highxz1[0] + sIxz * uxz[0];
    I0xz[1] = highxz1[1] + sIxz * uxz[1];// compute S1 intersect point

    //I1xz[0] = highxz2[0] + tIxz * vxz[0];
    //I1xz[1] = highxz2[1] + tIxz * vxz[1];// compute S2 intersect point
         
    if (TMath::Abs(I0xz[0]) < 800) {
      std::cout << "There is an XZ intersect on track AT point (" << I0xz[0] << ", 0, " << I0xz[1] << ")." << std::endl;


      // std::cout << "There is an XZ intersect on track AT point (" << I1xz[0] << ", 0, " << I1xz[1] << ")." << std::endl;

      return true;
    }

    return false;
  } // End of ELSE loop computing the Intersect point on ith and jth track

  return false;
}


bool air::AirKalmanAna::DoTracksIntersectYZ(TVector3 &highest1, TVector3 &lowest1, TVector3 &highest2, TVector3 &lowest2, TVector3 &dirvec1, TVector3 &dirvec2) const
{
  double    uyz[2];
  uyz[0] = dirvec1.Y();
  uyz[1] = dirvec1.Z();
        
  double    vyz[2];
  vyz[0] = dirvec2.Y();
  vyz[1] = dirvec2.Z();
        
  double    wyz[2];
  wyz[0] = highest1.Y() - highest2.Y();
  wyz[1] = highest1.Z() - highest2.Z();

  double highyz1[2] = {highest1.Y(), highest1.Z()};
        
  double I0yz[2] = {0.};//, I1yz[2];

  double     Dyz = uyz[0]*vyz[1] - uyz[1]*vyz[0];
        
  // the segments are skew and may intersect in a point

  // get the intersect parameter for S1
  double     sIyz = (vyz[0]*wyz[1] - vyz[1]*wyz[0]) / Dyz;

  // get the intersect parameter for S2
  double     tIyz = (uyz[0]*wyz[1] - uyz[1]*wyz[0]) / Dyz;
  
  if ( (sIyz < 0 || sIyz > 1) &&
       (tIyz < 0 || tIyz > 1) )                // no intersect with S1 or S2
  {
    //std::cout << "Track " << i << " does not intersect at tIxz < 0 || tIxz > 1" << std::endl;

    //std::cout << "Track " << j << " does not intersect at tIxz < 0 || tIxz > 1" << std::endl;
    
    return false;
  }
  else if ( (sIyz >= 0 && sIyz <= 1) &&
            (tIyz >= 0 && tIyz <= 1) ) {
     
    I0yz[0] = highyz1[0] + sIyz * uyz[0];
    I0yz[1] = highyz1[1] + sIyz * uyz[1];// compute S1 intersect point

    
    //I1yz[0] = highyz2[0] + tIyz * vyz[0];
    //I1yz[1] = highyz2[1] + tIyz * vyz[1];// compute S2 intersect point
         

    if (TMath::Abs(I0yz[0]) < 800 ) {
      std::cout << "There is a YZ intersect on track AT point (0," << I0yz[0] << ", " << I0yz[1] << ")." << std::endl;


      //std::cout << "There is a YZ intersect on track AT point (0," << I1yz[0] << ", " << I1yz[1] << ")." << std::endl;
  
      
      return true;
    }
    // USE Track.ID() -> Make a boolean function 
    // Fine me Track ID -> 
  } // End of Else statement for checking intersects, computing i'th and j'th track intersect

  return false;
}


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)
{

  for (unsigned int j = i+1; j < (tracks.size() - 1); ++j)
  {
    art::Ptr<rb::Track> track = tracks[j];

    TVector3 recoStart1, recoStart2, recoEnd1, recoEnd2;

          
    recoStart1.SetXYZ(ftrack->Start().X(),
                      ftrack->Start().Y(),
                      ftrack->Start().Z());
            
    recoEnd1.SetXYZ(ftrack->Stop().X(),
                    ftrack->Stop().Y(),
                    ftrack->Stop().Z());



    TVector3 highest1, lowest1;
    if (recoStart1.Y() > recoEnd1.Y())
    {
      highest1 = recoStart1;
      lowest1 = recoEnd1;
    }
    else
    {
      highest1=recoEnd1;
      lowest1=recoStart1;
    }

    recoStart2.SetXYZ(track->Start().X(),
                      track->Start().Y(),
                      track->Start().Z());
            
    recoEnd2.SetXYZ(track->Stop().X(),
                    track->Stop().Y(),
                    track->Stop().Z());

    TVector3 highest2, lowest2;
    if (recoStart2.Y() > recoEnd2.Y())
    {
      highest2 = recoStart2;
      lowest2 = recoEnd2;
    }
    else
    {
      highest2=recoEnd2;
      lowest2=recoStart2;
    }

    TVector3 dirvec1, dirvec2, diffvec;

    dirvec1 = highest1-lowest1;
    dirvec2 = highest2 - lowest2;

    double slope1Y = (highest1.Y() - lowest1.Y())/(highest1.Z() - lowest1.Z());

    double slope2Y = (highest2.Y() - lowest2.Y())/(highest2.Z() - lowest2.Z());

        
    double slope1X = (highest1.X() - lowest1.X())/(highest1.Z() - lowest1.Z());

    double slope2X = (highest2.X() - lowest2.X())/(highest2.Z() - lowest2.Z());

    //std::cout << "The (YZ) SLOPE OF TRACK " << i << " = " << slope1Y << " and Track " << j << " = " << slope2Y << std::endl;

    double signtestX = slope1X*slope2X;

    double signtestY = slope1Y*slope2Y;

                
    if (signtestX < 0 || signtestY < 0)
    {

      std::cout << "Track " << i << " and Track " << j << " have opposite slopes. " << std::endl;

      ++oppsigntrks;
          
    }
    else {

      ++samesigntrks;
         
    } 

  } // end of j track loop

  return;


}


void air::AirKalmanAna::analyze(art::Event const & e)
{

  art::Handle<std::vector<rb::CellHit> > hits;
  e.getByLabel(hit_label_, hits);

  // Get Slices 

  art::Handle< std::vector<rb::Cluster> > slicecol;     
  e.getByLabel(slice_label_,slicecol);
    
    
  // Make Slice List
  art::PtrVector<rb::Cluster> slicelist;        
  for(unsigned int i = 0; i<slicecol->size();++i)
  {     
    art::Ptr<rb::Cluster>slice(slicecol,i);
    slicelist.push_back(slice);
  }



  // setup the associations between tracks and slices   
  art::FindManyP<rb::Track> fmtrack(slicecol,e,track_label_);

  if (fRun < 0) fRun = e.run();
  if (fSubRun < 0) fSubRun = e.subRun();

   std::cout << "BEGIN ANALYZING EVENT " << e.event() << " NOW." << std::endl;

  art::Timestamp ts = e.time();
  // make it into a TTimeStamp
  TTimeStamp tts(ts.timeHigh(), ts.timeLow());
  
  TTimeStamp nowevent = tts;
  
  et_["time"] = tts.AsDouble();
  et_["timestamp"] = tts;
  event_tree_->Fill();
  // std::cout << "Approximate live time so far is " << nowevent.AsDouble() << " seconds." << std::endl;

  ++fNevents; // Count the number of events in the subrun


  // const time_t timeSec = ts.timeHigh();
  // struct tm* timeStruct = localtime(&timeSec);

  // unsigned short year            = timeStruct->tm_year + 1900;
  // unsigned short month           = timeStruct->tm_mon + 1;
  // unsigned short day             = timeStruct->tm_mday;
  // unsigned short doy             = timeStruct->tm_yday + 1;
  // unsigned short hour            = timeStruct->tm_hour + 1;
  // unsigned short minute          = timeStruct->tm_min;
  // unsigned short second          = timeStruct->tm_sec;
  

  std::ofstream outfile("OutFile_CosmicExposure.txt", std::ios::app);
  outfile << " *************************************************************** " << std::endl;

  
  ///////////////////////////////////////////////////////////////////////

  // Defining Detector Dimensions

  detHalfWidth                 =  geom->DetHalfWidth();
  detHalfHeight                =  geom->DetHalfHeight();
  detLength                    =  geom->DetLength(); 

  ///////////////////////////////////////////////////////////////////////
  // std::cout<<"WIDTH  :"
  //       <<detHalfWidth
  //       <<"  ,   HEIGHT  :"
  //       <<detHalfHeight
  //       <<"  ,  LENGTH  :"
  //       <<detLength<<std::endl;
  ///////////////////////////////////////////////////////////////////////
  // std::cout << "phi is " << phideg << " degrees." << std::endl;

  /////////////////////////////////////////////////////////////////////////////////////////////////


  // To Identify if this event is a golden event or not
  // To use a predefined list of specific events of interest
  // (hard coded above in GoldenEvents vector for DDenergy data from FD run 16820 taken on August 16th, 2014) 
  // file found at /nova/ana/users/mfrank/data <-- But No longer available. 

  // unsigned GEvent = 1;

  // for (unsigned i=0; i < Geventcol.size(); ++i){
  //   if (e.event() == Geventcol[i]){
  //     GEvent = 0;
  //     std::cout << "! ! ! THIS IS A GOLDEN EVENT ! ! !" << std::endl;
  //   }
  // } 

  /////////////////////////////////////////////////////////////////////////////////////////////////
 
  nslices =  0;


  /////////////////////////////////////////////////////////////////////////////////////////////////

  
  
  for (unsigned int islice = 0; islice < slicelist.size(); ++islice)
  {
    
    std::vector<art::Ptr<rb::Track> > tracks;   // A container for all tracks in the slice


    if(fmtrack.isValid())
    {   
      tracks = fmtrack.at(islice);      
    }

    //bool filtered = rb::IsFiltered(e, slicelist[islice]);
    
    //CheckAssociations(slicelist[islice],tracks, filtered);

    std::vector<art::Ptr<rb::Track> > n3Dtracks;
    std::vector<art::Ptr<rb::Track> > passedtracks; // A container for all tracks that passed containment in the slice

    std::vector<art::Ptr<rb::Track> > nocrosstracks; // A container for tracks that do not cross in a slice

    std::vector<art::Ptr<rb::Track> > sameslopetracks; // container for tracks that have the same signed slope

    std::vector<art::Ptr<rb::Track > > minlengthtracks;

    double scalarproduct;

    double deltatheta;
    
    std::vector<TVector3> alltracksvec;
    std::vector<TVector3> unittracks;
    std::vector<TVector3> passedveccol;
    std::vector<TVector3> samesignveccol;
    std::vector<TVector3> nocrossveccol;
    std::vector<TVector3> minlengthveccol;    
    std::vector<TVector3> sigmapassedveccol;
    std::vector<art::Ptr<rb::Track> > sigmapassedtracks;
    

    if(slicelist[islice]->IsNoise()) continue; // skipping over noise slices 

    ++nslices;
    // ADD condition "|| GEvent == 1)" to only get events of interest } 

    //if(slicelist[islice]->TotalADC() < 375000) continue;
    
    //if(GEvent == 1) continue;
   


    // std::vector<double>   thetaYcol;
    // std::vector<double>   thetaXcol;
    // std::vector<double>   thetaZcol;
    // std::vector<double>   thetaYYcol;
    // std::vector<double>   thetaYYYcol;



    // ===================================================================
    // First Cut: Check if Tracks Are "Contained." Want tracks that aren't.
    // ===================================================================

   
    // std::cout << "........ Now Testing Containment of Tracks." << std::endl;
        
    
    for (unsigned int i = 0; i < tracks.size(); ++i)

    { 
      art::Ptr<rb::Track> track = tracks[i];
      if (track->Is3D())
      {
        TVector3 unitvec;
        TVector3 unitvec_newxz;
            
        Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
        z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
     
        trkvec = unitvec_newxz;

        alltracksvec.push_back(trkvec);
        n3Dtracks.push_back(track);
        tottracks->Fill(track->TotalADC());
      }

    }
   
    
    fNTracks += n3Dtracks.size();
    if(n3Dtracks.size() < 2) continue;

    fnAll += n3Dtracks.size();
    fTracksPerSlice->Fill(n3Dtracks.size());

    for (unsigned int i=0; i< n3Dtracks.size(); ++i)
    {
      art::Ptr<rb::Track> track = n3Dtracks[i];
      
      TVector3 unitvec;
      TVector3 unitvec_newxz;
            
      Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
      z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
     
      trkvec = unitvec_newxz;
      trkcol = alltracksvec;
      

      std::vector<double> sigmacol;
        
      for (unsigned int j = i+1; j < (trkcol.size() - 1); j++) {
        scalarproduct = (trkcol[j].X()*trkvec.X()) + (trkcol[j].Y()*trkvec.Y()) + (trkcol[j].Z()*trkvec.Z());
        deltatheta = (acos(scalarproduct))*180/M_PI;

        //std::cout << "TRACK " << i << " and TRACK " << j << " have deltatheta = " << deltatheta << std::endl;
          
        fAllDeltaTheta->Fill(deltatheta);
      } // end of for trkcol 
      
    } // end of for n3D tracks


    for (unsigned int i=0;i<n3Dtracks.size();i++)
    {
      art::Ptr<rb::Track> track = n3Dtracks[i];
       

      if (track->Is3D() )
      {
        ///////////////////////////////////////////////////////////////////
        /// FROM TrackFit/CosmicTrackAna_module.cc


              
        // Get the end point of the track in a double array
        double xyz[3] = { track->Stop().X(), track->Stop().Y(), track->Stop().Z() };
        double xyzstart[3] = {track->Start().X(), track->Start().Y(), track->Start().Z() };

        
        double highest1[3], lowest1[3];
        if (xyzstart[1] > xyz[1])
        {
          highest1[0] = xyzstart[0];
          highest1[1] = xyzstart[1];
          highest1[2] = xyzstart[2];
          lowest1[0] = xyz[0];
          lowest1[1] = xyz[1];
          lowest1[2] = xyz[2];
        }
        else
        {
          highest1[0]=xyz[0];
          highest1[1]=xyz[1];
          highest1[2]=xyz[2];
          lowest1[0]=xyzstart[0];
          lowest1[1]=xyzstart[1];
          lowest1[2]=xyzstart[2];
        }
        
              
        // OLD definition of distance to edge of detector (NOW people use LiveGeometry)
        double dist2edgeOLD = geo::DistToEdge(lowest1,                                      
                                              detHalfWidth,                                         
                                              detHalfHeight,                                        
                                              detLength); 
        
              
        double dist2edgeStart = geo::DistToEdge(highest1,                                           
                                                detHalfWidth,                                       
                                                detHalfHeight,                                      
                                                detLength); 

        //////// LENGTH  //////////////////////
      
        recoStart.SetXYZ(highest1[0],
                         highest1[1],
                         highest1[2]);
      
        recoEnd.SetXYZ(lowest1[0],
                       lowest1[1],
                       lowest1[2]);


        TVector3 length = recoEnd - recoStart;
        totallength = length.Mag()/100;

        //if (totallength < 5.0 ) continue; (might need to use length cut later. I don't want to be too restrictive earlier on)
                                            
        /////// CONTAINMENT CUTS ////////////////
                                              
        // Let's define a fiducial box 
        // (i.e. a box where start point of track must be)
        if (dist2edgeStart > 100.0){  
          // Let's define a containment box
          // (i.e. a box where end point of track must be)
          if (dist2edgeOLD > 100.0){         
            // function returns a value greater 
            // than 25.0 cm inside detector in order to be contained 
            continue;
          } // end of not fiducial and no contained loop -> Tracks to get rid of
        } // end of Fiducial if Loop
        else {
          passedtracks.push_back(track);
          
        } // End of Fiducial Else loop  
      } // end Is3D loop    
    } // end For (Tracks) loop


    // ===================== END CONTAINMENT LOOP ========================
      
    // Only looking at slices with more than 1 containment passed track
    if (passedtracks.size() < 2) continue;
    fnCont += passedtracks.size();
    for (unsigned int i=0; i<passedtracks.size(); ++i)
    {

      art::Ptr<rb::Track> track = passedtracks[i];

      TVector3 unitvec;
      TVector3 unitvec_newxz;
            
      Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
      z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
     
      trkvec = unitvec_newxz;

      passedveccol.push_back(trkvec);

    }

    /// MAKE THE FIRST DELTATHETA PLOT FOR CONTAINMENT PASSED TRACKS

    std::cout << "............Now Calculating DeltaTheta of PassedTracks." << std::endl;
    
    for (unsigned int i=0; i<passedtracks.size(); ++i)
    {
      art::Ptr<rb::Track> track = passedtracks[i];

      TVector3 unitvec;
      TVector3 unitvec_newxz;
            
      Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
      z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
     
      trkvec = unitvec_newxz;
      trkcol = passedveccol;

      std::vector<double> sigmacol;

      for (unsigned int j = i+1; j < (trkcol.size() - 1); j++) {
        scalarproduct = (trkcol[j].X()*trkvec.X()) + (trkcol[j].Y()*trkvec.Y()) + (trkcol[j].Z()*trkvec.Z());
        deltatheta = (acos(scalarproduct))*180/M_PI;
        fContDeltaTheta->Fill(deltatheta);
      }
      conttracks->Fill(track->TotalADC());
    }


    // ===================================================================
    // Second Cut : Check if Tracks Intersect in XZ or YZ View. 
    // ===================================================================

    std::cout << "...... Now Testing Intersection between track " << std::endl;

    for (unsigned int i = 0; i < passedtracks.size(); ++i)

    { 
      art::Ptr<rb::Track> track1 = passedtracks[i];

      unsigned int XZcrosscount = 0;
      unsigned int YZcrosscount = 0;

      unsigned int noXZcross = 0;
      unsigned int noYZcross = 0;

      for (unsigned int j = i + 1; j < (passedtracks.size() - 1); ++j)
      {
        art::Ptr<rb::Track> track2 = passedtracks[j];

        std::vector<art::Ptr<rb::Track> > trkcollection = passedtracks;

        TVector3 recoStart1, recoStart2, recoEnd1, recoEnd2;

        //trackcollection = passedtracks;
          
        recoStart1.SetXYZ(track1->Start().X(),
                          track1->Start().Y(),
                          track1->Start().Z());
            
        recoEnd1.SetXYZ(track1->Stop().X(),
                        track1->Stop().Y(),
                        track1->Stop().Z());



        TVector3 highest1, lowest1;
        if (recoStart1.Y() > recoEnd1.Y())
        {
          highest1 = recoStart1;
          lowest1 = recoEnd1;
        }
        else
        {
          highest1=recoEnd1;
          lowest1=recoStart1;
        }

        recoStart2.SetXYZ(track2->Start().X(),
                          track2->Start().Y(),
                          track2->Start().Z());
            
        recoEnd2.SetXYZ(track2->Stop().X(),
                        track2->Stop().Y(),
                        track2->Stop().Z());

        TVector3 highest2, lowest2;
        if (recoStart2.Y() > recoEnd2.Y())
        {
          highest2 = recoStart2;
          lowest2 = recoEnd2;
        }
        else
        {
          highest2=recoEnd2;
          lowest2=recoStart2;
        }

        TVector3 dirvec1, dirvec2, diffvec;

        dirvec1 = highest1-lowest1;
        dirvec2 = highest2 - lowest2;


        if (DoTracksIntersectYZ(highest1, lowest1, highest2, lowest2, dirvec1, dirvec2) )
        {

          //std::cout << "There is a YZ intersect on track " << i << " with track  " << j << " in Slice " << nslices << std::endl;

          ++YZcrosscount;
          //continue;
        }
        else {

          ++noYZcross;

        }

        
        if (DoTracksIntersectXZ(highest1, lowest1, highest2, lowest2, dirvec1, dirvec2) ) {

          //std::cout << "There is a XZ intersect on track " << i << " with track  " << j << " in Slice " << nslices << std::endl;
        
            ++XZcrosscount;
            //continue;
        }
        else {

          //std::cout << "There are no XZ intersects on track " << i << " with track  " << j << " in Slice " << nslices << std::endl;
            
            ++noXZcross;
            
        }

      } // End of j'th track loop 

      unsigned int sumcross = YZcrosscount + XZcrosscount;

      unsigned int sumnocross = noXZcross + noYZcross;

      if (sumcross < sumnocross || sumcross == 0)

      {
        //std::cout << "TRACK " << i << " DOES NOT INTERSECT WITH ANY TRACKS ON SLICE " << nslices << std::endl;
        nocrosstracks.push_back(track1);
      }

      if (sumcross > sumnocross)

      {
        std::cout << "TRACK " << i << " INTERSECTS " << sumcross << " TIMES IN XZ OR YZ VIEW. " << std::endl;
        continue;

      }



    } // End of i'th tracks loop

    if (nocrosstracks.size() < 2) continue;

    fnInt += nocrosstracks.size();

    for (unsigned int i = 0; i < nocrosstracks.size(); ++i)

    { 
      art::Ptr<rb::Track> track = nocrosstracks[i];


      TVector3 unitvec;
      TVector3 unitvec_newxz;
            
      Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
      z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
     
      trkvec = unitvec_newxz;

      nocrossveccol.push_back(trkvec);
      inttracks->Fill(track->TotalADC());
    }

    std::cout << "............Now Calculating DeltaTheta of NoCrossTracks." << std::endl;
    
    //std::cout << "There are " << nocrosstracks.size() << "number of tracks in slice " << nslices << std::endl;

    for (unsigned int i=0; i<nocrosstracks.size(); ++i)
    {
      art::Ptr<rb::Track> track = nocrosstracks[i];

      TVector3 unitvec;
      TVector3 unitvec_newxz;
            
      Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
      z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
     
      trkvec = unitvec_newxz;
      trkcol = nocrossveccol;

      std::vector<double> sigmacol;
      double deltatheta;

      for (unsigned int j = i+1; j < (trkcol.size() - 1); j++) {
        scalarproduct = (trkcol[j].X()*trkvec.X()) + (trkcol[j].Y()*trkvec.Y()) + (trkcol[j].Z()*trkvec.Z());
        deltatheta = (acos(scalarproduct))*180/M_PI;
        
        fIntDeltaTheta->Fill(deltatheta);
      }

    }
    
    //=====================================================================

    // FINDING DISTANCE BETWEEN TRACKS 

    // 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.

    // http://geomalgorithms.com/a05-_intersect-1.html 
    
    // http://geomalgorithms.com/a07-_distance.html#dist3D_Segment_to_Segment 

    
    //======================================================================


    std::vector<art::Ptr<rb::Track > > samesigncol;


    std::cout << ".........Now Checking Slope Signs." << std::endl;

    for (unsigned int i = 0; i < nocrosstracks.size(); ++i)

    { 
      art::Ptr<rb::Track> track1 = nocrosstracks[i];

      art::Ptr<rb::Track> ftrack = track1;

      std::vector<art::Ptr<rb::Track > > tracks = nocrosstracks;
      

      std::vector<art::Ptr<rb::Track> > trkcollection = nocrosstracks;

      unsigned int oppsigntrks = 0;
      unsigned int samesigntrks = 0;

      
      OppSignTrkCount(ftrack, tracks, oppsigntrks, samesigntrks, i);

      if (oppsigntrks > samesigntrks) {

        std::cout << "Track " << i << " on slice " << nslices << "has opposite slope to " << oppsigntrks << " other tracks. " << std::endl;
        std::cout << "TRACK " << i << " IS A BAD TRACK." << std::endl;
        continue;
      }

      else {

        samesigncol.push_back(ftrack);

      }
 

    } // end of i track loop

    // if (opptracks > 0) continue;

    if (samesigncol.size() < 2) continue;
          
    fnOpp += samesigncol.size();

    for (unsigned int i = 0; i < samesigncol.size(); ++i)

    { 
      art::Ptr<rb::Track> track = samesigncol[i];


      TVector3 unitvec;
      TVector3 unitvec_newxz;
            
      Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
      z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
     
      trkvec = unitvec_newxz;

      samesignveccol.push_back(trkvec);
      opptracks->Fill(track->TotalADC());
    }

    std::cout << "............Now Calculating DeltaTheta of SameSignTracks." << std::endl;
    
    for (unsigned int i=0; i<samesigncol.size(); ++i)
    {
      art::Ptr<rb::Track> track = samesigncol[i];

      TVector3 unitvec;
      TVector3 unitvec_newxz;
            
      Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
      z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
     
      trkvec = unitvec_newxz;
      trkcol = samesignveccol;

      for (unsigned int j = i+1; j < (trkcol.size() - 1); j++) {
        scalarproduct = (trkcol[j].X()*trkvec.X()) + (trkcol[j].Y()*trkvec.Y()) + (trkcol[j].Z()*trkvec.Z());
        deltatheta = (acos(scalarproduct))*180/M_PI;
        fSlopeDeltaTheta->Fill(deltatheta);
      }

      if (track->TotalLength() > 200.0) {
        minlengthtracks.push_back(track);
      }
    }

    if (minlengthtracks.size() < 2) continue;
    fnLen += minlengthtracks.size();

    for (unsigned int i = 0; i < minlengthtracks.size(); ++i)

    { 
      art::Ptr<rb::Track> track = minlengthtracks[i];


      TVector3 unitvec;
      TVector3 unitvec_newxz;
            
      Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
      z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
     
      trkvec = unitvec_newxz;

      minlengthveccol.push_back(trkvec);

    }

    std::cout << "............Now Calculating DeltaTheta of MinLengthTracks." << std::endl;
    

    for (unsigned int i=0; i < minlengthtracks.size(); i++)
    {
      art::Ptr<rb::Track> track = minlengthtracks[i];

      TVector3 unitvec;
      TVector3 unitvec_newxz;
            
      Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
      z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)
     
      trkvec = unitvec_newxz;
      trkcol = minlengthveccol;

      for (unsigned int j = i+1; j < (minlengthtracks.size() - 1); j++) {
        scalarproduct = (trkcol[j].X()*trkvec.X()) + (trkcol[j].Y()*trkvec.Y()) + (trkcol[j].Z()*trkvec.Z());
        deltatheta = (acos(scalarproduct))*180/M_PI;
        //std::cout << "TRACK " << i << " and TRACK " << j << " have deltatheta = " << deltatheta << std::endl;

        fLenDeltaTheta->Fill(deltatheta);
      }
      lentracks->Fill(track->TotalADC());
    }
    
    std::cout << ".......... Now Calculating Closest Distance Between Tracks" << std::endl;

    std::vector< double > DOCAmincol;

    for (unsigned int i = 0 ; i < minlengthtracks.size() ; ++i)
    {
      art::Ptr<rb::Track> track1 = minlengthtracks[i];
      
      std::vector<double> DOCAcol;
      for (unsigned int j = i + 1; j < (minlengthtracks.size() - 1); j++)
      {
        art::Ptr<rb::Track> track2 = minlengthtracks[j];

        TVector3 recoStart1, recoStart2, recoEnd1, recoEnd2;
          
        recoStart1.SetXYZ(track1->Start().X(),
                          track1->Start().Y(),
                          track1->Start().Z());
            
        recoEnd1.SetXYZ(track1->Stop().X(),
                        track1->Stop().Y(),
                        track1->Stop().Z());

        TVector3 highest1, lowest1;
        if (recoStart1.Y() > recoEnd1.Y())
        {
          highest1 = recoStart1;
          lowest1 = recoEnd1;
        }
        else
        {
          highest1=recoEnd1;
          lowest1=recoStart1;
        }

        recoStart2.SetXYZ(track2->Start().X(),
                          track2->Start().Y(),
                          track2->Start().Z());
            
        recoEnd2.SetXYZ(track2->Stop().X(),
                        track2->Stop().Y(),
                        track2->Stop().Z());

        TVector3 highest2, lowest2;
        if (recoStart2.Y() > recoEnd2.Y())
        {
          highest2 = recoStart2;
          lowest2 = recoEnd2;
        }
        else
        {
          highest2=recoEnd2;
          lowest2=recoStart2;
        }

        TVector3 dirvec1, dirvec2, diffvec;

        dirvec1 = highest1-lowest1;
        dirvec2 = highest2 - lowest2;

        diffvec = highest1 - highest2;

      

        double a = dirvec1.X()*dirvec1.X() + dirvec1.Y()*dirvec1.Y() + dirvec1.Z()*dirvec1.Z(); // Always >= 0
        double b = dirvec1.X()*dirvec2.X() + dirvec1.Y()*dirvec2.Y() + dirvec1.Z()*dirvec2.Z();

        double c = dirvec2.X()*dirvec2.X() + dirvec2.Y()*dirvec2.Y() + dirvec2.Z()*dirvec2.Z(); // Always >= 0
        
        double d = dirvec1.X()*diffvec.X() + dirvec1.Y()*diffvec.Y() + dirvec1.Z()*diffvec.Z();

        double e = dirvec2.X()*diffvec.X() + dirvec2.Y()*diffvec.Y() + dirvec2.Z()*diffvec.Z();
        
        double D = a*c - b*b; // Always >=0 

        double    sc, sN, sD = D;       // "steps" sc = sN / sD, default sD = D >= 0
        double    tc, tN, tD = D;       // tc = tN / tD, default tD = D >= 0


        if (D < SMALL_NUM) { // The lines are almost parallel

          sN = 0.0; // force using point Highest1 on linesegment (track1)
          
          sD = 1.0; // to prevent possible division by 0.0 later

          tN = e;

          tD = c;

        }
        else { // get the closest points on the infinite lines
          sN = (b*e - c*d);
          tN = (a*e - b*d);
          if (sN < 0.0) {        // sc < 0 => the s=0 edge is visible
            sN = 0.0;
            tN = e;
            tD = c;
          }
          else if (sN > sD) {  // sc > 1  => the s=1 edge is visible
            sN = sD;
            tN = e + b;
            tD = c;
          }
        }

        if (tN < 0.0) {            // tc < 0 => the t=0 edge is visible
          tN = 0.0;
          // recompute sc for this edge
          if (-d < 0.0)
            sN = 0.0;
          else if (-d > a)
            sN = sD;
          else {
            sN = -d;
            sD = a;
          }
        }
        else if (tN > tD) {      // tc > 1  => the t=1 edge is visible
          tN = tD;
          // recompute sc for this edge
          if ((-d + b) < 0.0)
            sN = 0;
          else if ((-d + b) > a)
            sN = sD;
          else {
            sN = (-d +  b);
            sD = a;
          }
        }
        // finally do the division to get sc and tc
        sc = (TMath::Abs(sN) < SMALL_NUM ? 0.0 : sN / sD);
        tc = (TMath::Abs(tN) < SMALL_NUM ? 0.0 : tN / tD);

        // get the difference of the two closest points
        TVector3   dP = diffvec + (sc * dirvec1 ) - (tc * dirvec2);  // =  S1(sc) - S2(tc)

        double   closestdist = dP.Mag();   // return the closest distance


        DOCAcol.push_back(closestdist);
        

        // 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? 

        // 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? 

      } // end for J tracks loop

      // Get minimum DOCA value. 
      double DOCAmin = 10000;
      for (unsigned int k = 0; k < DOCAcol.size(); k++)
      {
        if (DOCAcol[i] < DOCAmin)
          DOCAmin = DOCAcol[i];
      }

      DOCAmincol.push_back(DOCAmin);
      
      double DOCAmeter = DOCAmin/100.0;
      fDOCA->Fill(DOCAmeter);
  
    } // end For (Tracks i ) loop



    /////////////////////////////////////////////////////////////////////

    // HOW TO GET ANGLES BETWEEN INDIVIDUAL TRACKS

    /////////////////////////////////////////////////////////////////////

   
  
    std::vector<double> azimuthcol;
    std::vector<double> altcol;
    std::vector<double> sigmacol;

    for (unsigned int i = 0; i < minlengthtracks.size(); i++){
      
      art::Ptr<rb::Track> track = minlengthtracks[i];
      
      TVector3 unitvec;
      TVector3 unitvec_newxz;
            
      Unitize(track, unitvec); // Obtain unit vector of reco tracks based on start and end points
      z2north(unitvec, unitvec_newxz); // obtain new unitvectors with Z rotated relative to the true north (phi)

      trkvec = unitvec_newxz;


      trkcol = minlengthveccol;


      for (unsigned int j = i+1; j < (trkcol.size() - 1); j++) {
        scalarproduct = (trkcol[j].X()*trkvec.X()) + (trkcol[j].Y()*trkvec.Y()) + (trkcol[j].Z()*trkvec.Z());
        deltatheta = (acos(scalarproduct))*180/M_PI;
        double deltaangle = (deltatheta*deltatheta)/2;
        fDeltaAngle->Fill(deltaangle);
      }

      //////// LENGTH  //////////////////////
      
      recoStart.SetXYZ(track->Start().X(),
                       track->Start().Y(),
                       track->Start().Z());
      
      recoEnd.SetXYZ(track->Stop().X(),
                     track->Stop().Y(),
                     track->Stop().Z());

      TVector3 highest, lowest;
      if (recoStart.Y() > recoEnd.Y())
      {
        highest = recoStart;
        lowest = recoEnd;
      }
      else
      {
        highest=recoEnd;
        lowest=recoStart;
      }

      TVector3 length = highest - lowest;
      totallength = length.Mag()/100;
      

      double alt = 90 - acos(trkvec.Y())*180/M_PI; // Altitude should just be complimentary to the zenith angle (90-zenith)

      // The following angle is based on the XZ projection of the track. 
      // Not sure if this is the right method. 
      // Since legnth of the unitvector is just 1, r = sin(Az). 
      // Where r = magnitude of XZ projection

      double azimuth = asin(trkvec.Z()/TMath::Sqrt( (trkvec.X()*trkvec.X()) + (trkvec.Z()*trkvec.Z()) ) )*180/M_PI; 
      azimuthcol.push_back(azimuth);
      altcol.push_back(alt);

      fAzimuth->Fill(azimuth);
      fAltitude->Fill(alt);

    }
    

    fTracksPerSliceFinal->Fill(minlengthtracks.size());
    pramcol = sigmacol;
    AvgValue(pramcol, avgpram);
    
    double avgsigma = avgpram;
    
    std::cout << "Average sigma per slice is " << avgsigma << std::endl;

    std::cout << "...............Finished Processing SLICE " << nslices << " with " << minlengthtracks.size() << " passed tracks in event " << e.event() << std::endl;

  } // End Slice loop


  fNSlices->Fill(nslices);

    
  std::cout << "******FINISHED PROCESSING EVENT NUMBER  " << e.event() << " ******" << std::endl;

  std::cout << "************************************************************************" << std::endl;
  
  std::cout << "There are: " << fnAll << " total tracks in this job so far. " << std::endl;

  std::cout << "There are: " << fnCont << " containment passed tracks in this job so far. " << std::endl;

  std::cout << "There are: " << fnInt << " intersect passed tracks in this job so far. " << std::endl;

  std::cout << "There are: " << fnOpp << " slope passed tracks in this job so far. " << std::endl;

  std::cout << "There are: " << fnLen << " length passed tracks in this job so far. " << std::endl;
  
  std::cout << "There are: " << fNevents << " events in this job so far. " << std::endl;

    
  std::cout << "There are: " << fNTracks << " total tracks in this job so far. " << std::endl;
  


 
  
} // End of Analyze Event Loop


//......................................................................
        
void air::AirKalmanAna::beginSubRun(const art::SubRun & sr)     
{
  art::Handle<sumdata::CosmicExposure> ce;      
  sr.getByLabel("exposure", ce);
  if(ce.failedToGet()){
    mf::LogError("AirKalmanAna") << "No Cosmic Exposure Info";
    return;
    }
  fLiveTime += ce->totlivetime; 
  LOG_INFO("AirKalmanAna") << "The live time is now " << fLiveTime;
  return;
}
        
//......................................................................




//......................................................................  

void air::AirKalmanAna::endJob()

{ 

  fExposure->Fill(fLiveTime);


  if(fLiveTime != 0){

    double cosmicrate = fNTracks/(1000.0*fLiveTime);

    srt_["cosmicrate"]=cosmicrate;
    srt_["exposure"]=fLiveTime;
    srt_["tottracks"]=fNTracks;
    srt_["passedevets"]=fNevents;
    srt_["passedtracks"]=fnLen;



    LOG_VERBATIM("AirShower") << "There are: " << fNTracks 

                                   << " cosmic tracks selected in this job.\n"

                                   << "Livetime is " << fLiveTime 

                                   << " according to CosmicExposure.\n"

                                   << "Note, livetime is for FULL SUBRUN, "

                                   << "so rate is wrong if haven't run full subrun.\n"

                                   << "The cosmic rate is therefore " 

                                   << cosmicrate << " Hz \n";


      

    fCosmicRate->Fill(cosmicrate);
    subrun_tree_->Fill();
  }
  return;       
} 




        // //==============================================================
        // //Code for Checking Intersects -> Don't need it since 
        // //I defined the function.
        // //===============================================================

        // std::cout << "CHECKING FUNCTION AGAINST CODE" << std::endl;

        // double    uxz[2];
        // uxz[0] = dirvec1.X();
        // uxz[1] = dirvec1.Z();
        
        // double    vxz[2];
        // vxz[0] = dirvec2.X();
        // vxz[1] = dirvec2.Z();
        
        // double    wxz[2];
        // wxz[0] = highest1.X() - highest2.X();
        // wxz[1] = highest1.Z() - highest2.Z();
        
        // double I0xz[2], I1xz[2]; // Define array for two intersect points

        // double     Dxz = uxz[0]*vxz[1] - uxz[1]*vxz[0];

        // double highxz1[2], highxz2[2];
        // highxz1[0] = highest1.X();
        // highxz1[1] = highest1.Z();
        // highxz2[0] = highest2.X();
        // highxz2[1] = highest2.Z();
        
        // // the segments are skew and may intersect in a point
        // // get the intersect parameter for S1
        // double     sIxz = (vxz[0]*wxz[1] - vxz[1]*wxz[0]) / Dxz;
        // if (sIxz < 0 || sIxz > 1)                // no intersect with S1
        // {
        //   //std::cout << "Track " << i << " does not intersect at sIxz < 0 || sIxz > 1" << std::endl;
        // }

        // // get the intersect parameter for S2
        // double     tIxz = (uxz[0]*wxz[1] - uxz[1]*wxz[0]) / Dxz;
        // if (tIxz < 0 || tIxz > 1)  // no intersect with S2
        // {
        //   //std::cout << "Track " << j << " does not intersect at tIxz < 0 || tIxz > 1" << std::endl;
        // }
        // else {
     
        //   I0xz[0] = highxz1[0] + sIxz * uxz[0];
        //   I0xz[1] = highxz1[1] + sIxz * uxz[1];// compute S1 intersect point

        //   I1xz[0] = highxz2[0] + tIxz * vxz[0];
        //   I1xz[1] = highxz2[1] + tIxz * vxz[1];// compute S2 intersect point
          
        //   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;


        //   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;

        //   ++ntrackscrossXZ; // Count number of XZ intersecting tracks
        //   ++ntrackscrossXZslice;

        //   continue;

        // } // End of ELSE loop computing the Intersect point on ith and jth track
        



            
        // //===================================================================
        
        // // COMPUTE INTERSECT FOR YZ track segments
        
        // //===================================================================
        // double    uyz[2];
        // uyz[0] = dirvec1.Y();
        // uyz[1] = dirvec1.Z();
        
        // double    vyz[2];
        // vyz[0] = dirvec2.Y();
        // vyz[1] = dirvec2.Z();
        
        // double    wyz[2];
        // wyz[0] = highest1.Y() - highest2.Y();
        // wyz[1] = highest1.Z() - highest2.Z();

        // double highyz1[2], highyz2[2];
        // highyz1[0] = highest1.Y();
        // highyz1[1] = highest1.Z();
        // highyz2[0] = highest2.Y();
        // highyz2[1] = highest2.Z();
        
        // double I0yz[2], I1yz[2];

        // double     Dyz = uyz[0]*vyz[1] - uyz[1]*vyz[0];
        
        // // the segments are skew and may intersect in a point
        // // get the intersect parameter for S1
        // double     sIyz = (vyz[0]*wyz[1] - vyz[1]*wyz[0]) / Dyz;
        // if (sIyz < 0 || sIyz > 1)                // no intersect with S1
        // {
        //   //std::cout << "Track " << i << " does not intersect at tIxz < 0 || tIxz > 1" << std::endl;
        // }

        // // get the intersect parameter for S2
        // double     tIyz = (uyz[0]*wyz[1] - uyz[1]*wyz[0]) / Dyz;
        // if (tIyz < 0 || tIyz > 1)  // no intersect with S2
        // {
        //   //std::cout << "Track " << j << " does not intersect at tIxz < 0 || tIxz > 1" << std::endl;
        // }
        // else {
     
        //   I0yz[0] = highyz1[0] + sIyz * uyz[0];
        //   I0yz[1] = highyz1[1] + sIyz * uyz[1];// compute S1 intersect point
          
        //   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;

        //   I1yz[0] = highyz2[0] + tIyz * vyz[0];
        //   I1yz[1] = highyz2[1] + tIyz * vyz[1];// compute S2 intersect point

        //   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;
                  
        //   ++ntrackscrossYZ; // count number of YZ intersecting tracks
        //   ++ntrackscrossYZslice;

        //   // USE Track.ID() -> Make a boolean function 
        //   // Fine me Track ID -> 
        // } // End of Else statement for checking intersects, computing i'th and j'th track intersect
        


  


DEFINE_ART_MODULE(air::AirKalmanAna)