RecoAnalysis_module.cc

Go to the documentation of this file.

#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/Assns.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "HoughTransform.h"
#include "NNbarUtilities.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Prong.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Vertex.h"
#include "Utilities/AssociationUtil.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"

#include <cmath>
#include <iterator>
#include <math.h>
#include <string>
#include <vector>

#include <TCanvas.h>
#include <TColor.h>
#include <TGraph.h>
#include <TH1.h>
#include <TH2.h>
#include <TLatex.h>
#include <TLegend.h>
#include <TMath.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TTree.h>

// struct Point2D {
//   point2d_t point;
//
//   Point2D(double xa, double ya) {
//     point.first = xa;
//     point.second = ya;
//   }
//
//   struct Point2D &operator+=(const Point2D &rhs) {
//     point.first += rhs.point.first;
//     point.second += rhs.point.second;
//     return *this;
//   }
//
//   struct Point2D &operator-=(const Point2D &rhs) {
//     point.first -= rhs.point.first;
//     point.second -= rhs.point.second;
//     return *this;
//   }
//
//   struct Point2D &operator*=(const double &k) {
//     point.first *= k;
//     point.second *= k;
//     return *this;
//   }
//
//   double operator*=(const Point2D &rhs) {
//     return point.first * rhs.point.first + point.second * rhs.point.second;
//   }
//
//   struct Point2D &operator=(const Point2D &rhs) {
//     point.first = rhs.point.first;
//     point.second = rhs.point.second;
//     return *this;
//   }
// };
//
// Point2D operator+(Point2D lhs, const Point2D &rhs) { return lhs += rhs; }
// Point2D operator-(Point2D lhs, const Point2D &rhs) { return lhs -= rhs; }
// double operator*(Point2D lhs, const Point2D &rhs) { return lhs *= rhs; }
// Point2D operator*(const double k, Point2D rhs) { return rhs *= k; }
// Point2D operator*(Point2D lhs, const double k) { return lhs *= k; }
//
// double length(Point2D A) {
//   return TMath::Sqrt(TMath::Power(A.point.first, 2) + TMath::Power(A.point.second, 2));
// }
//
// unsigned int check_isolation(double* array, unsigned int maxIndex, unsigned int hitIndex) {
//   if (TMath::Abs(*(array + hitIndex + 1) - *(array + hitIndex)) > 5.) {
//     return hitIndex;
//   } else {
//     if (hitIndex + 1 == maxIndex) {
//       return hitIndex;
//     } else {
//       return check_isolation(array, maxIndex, hitIndex + 1);
//     }
//   }
// }

namespace nnbar {
class RecoAnalysis;
}

class nnbar::RecoAnalysis : public art::EDAnalyzer {
public:
  explicit RecoAnalysis(fhicl::ParameterSet const &p);
  void beginJob() override;
  void endJob() override;
  void analyze(art::Event const &e) override;
  TGraph *PlotVectorPoint(std::vector<point2d_t> pts);
  double HozVertex(
      std::vector<std::tuple<double, double, double>> HitCollectionSingleView,
      double vertexZ);
  double slopeAngleOfLineConnect2Points(point2d_t pA, point2d_t pB);
  void evd(art::Event const &e, unsigned int clusterIdx);

private:
  double _ExposureTimeInUSec;
  art::InputTag _InputShwrLabel;
  unsigned int _SampleType;

  TTree *_VariableTree;
  unsigned int _tRunNumber;
  unsigned int _tSubrunNumber;
  unsigned int _tEventNumber;
  unsigned int _tSliceNumber;
  unsigned int _tSampleType;
  double _tTotalVisibleEnergy;
  double _tTemporalClusterExpand;
  double _tTotalHitCountInXView;
  double _tTotalHitCountInYView;
  double _tHitCountXYDifference;
  double _tAverageEnergyPerHitXView;
  double _tAverageEnergyPerHitYView;
  double _tHitVsVertexTimingCorrelation;
  double _tTrueVertexX;
  double _tTrueVertexY;
  double _tTrueVertexZ;
  double _tTrueVertexT;
  double _tReconVertexX;
  double _tReconVertexY;
  double _tReconVertexZ;
  double _tReconVertexT;
  double _tEArmVertexX;
  double _tEArmVertexY;
  double _tEArmVertexZ;
  double _tEArmVertexT;

  TH1D *EnergyRatioFromVertexX;
  TH1D *EnergyRatioFromVertexY;
  TH1D *EnergyRatioFromVertexZ;

  TH1D *TimingCorrelation;

  TH1D *DistanceFromReconVertexToTrueVertex;
  TH1D *ZErrorFromReconVertexToTrueVertex;
  TH2D *OrientationFromReconVertexToTrueVertex;
  TH1D *PolarAngleFromReconVertexToTrueVertex;
  TH1D *AzimuthalAngleFromReconVertexToTrueVertex;

  TH1D *DistanceFromElasticArmVertexToTrueVertex;
  TH1D *ZErrorFromElasticArmVertexToTrueVertex;
  TH2D *OrientationFromElasticArmVertexToTrueVertex;
  TH1D *PolarAngleFromElasticArmVertexToTrueVertex;
  TH1D *AzimuthalAngleFromElasticArmVertexToTrueVertex;

  std::vector<std::tuple<double, double, double>> HitsXView, HitsYView;
};

nnbar::RecoAnalysis::RecoAnalysis(fhicl::ParameterSet const &p)
    : EDAnalyzer(p) {
  _InputShwrLabel = "shower:showers";
  _SampleType = p.get<unsigned int>("SampleType");
}

// double nnbar::RecoAnalysis::HozVertex(
//     std::vector<std::tuple<double, double, double>> HitCollectionSingleView,
//     double vertexZ) {
//   std::vector<Point2D> Hits;
//   for (auto &hit : HitCollectionSingleView) {
//     double x = std::get<0>(hit);
//     double y = std::get<1>(hit);
//     Point2D aHit(x, y);
//     Hits.push_back(aHit);
//   }
//
//   return 0;
// }

void nnbar::RecoAnalysis::beginJob() {
  _ExposureTimeInUSec = 0.;

  art::ServiceHandle<art::TFileService> tfs;

  _VariableTree =
      tfs->make<TTree>("VariableTree", "NNbar Analysis Variable Tree");
  _VariableTree->Branch("RunNumber", &_tRunNumber);
  _VariableTree->Branch("SubrunNumber", &_tSubrunNumber);
  _VariableTree->Branch("EventNumber", &_tEventNumber);
  _VariableTree->Branch("SliceNumber", &_tSliceNumber);
  _VariableTree->Branch("SampleType", &_tSampleType);
  _VariableTree->Branch("TotalVisibleEnergy", &_tTotalVisibleEnergy);
  _VariableTree->Branch("TotalHitCountInXView", &_tTotalHitCountInXView);
  _VariableTree->Branch("TotalHitCountInYView", &_tTotalHitCountInYView);
  _VariableTree->Branch("HitCountXYDifference", &_tHitCountXYDifference);
  _VariableTree->Branch("TemporalClusterExpand", &_tTemporalClusterExpand);
  _VariableTree->Branch("AverageEnergyPerHitXView",
                        &_tAverageEnergyPerHitXView);
  _VariableTree->Branch("AverageEnergyPerHitYView",
                        &_tAverageEnergyPerHitYView);
  _VariableTree->Branch("HitVsVertexTimingCorrelation",
                        &_tHitVsVertexTimingCorrelation);
  _VariableTree->Branch("TrueVertexX", &_tTrueVertexX);
  _VariableTree->Branch("TrueVertexY", &_tTrueVertexY);
  _VariableTree->Branch("TrueVertexZ", &_tTrueVertexZ);
  _VariableTree->Branch("TrueVertexT", &_tTrueVertexT);
  _VariableTree->Branch("ReconVertexX", &_tReconVertexX);
  _VariableTree->Branch("ReconVertexY", &_tReconVertexY);
  _VariableTree->Branch("ReconVertexZ", &_tReconVertexZ);
  _VariableTree->Branch("ReconVertexT", &_tReconVertexT);
  _VariableTree->Branch("EArmVertexX", &_tEArmVertexX);
  _VariableTree->Branch("EArmVertexY", &_tEArmVertexY);
  _VariableTree->Branch("EArmVertexZ", &_tEArmVertexZ);
  _VariableTree->Branch("EArmVertexT", &_tEArmVertexT);

  EnergyRatioFromVertexX =
      tfs->make<TH1D>("EnergyRatioFromVertexX", ";Energy Ratio;", 100, 0, 1);
  EnergyRatioFromVertexY =
      tfs->make<TH1D>("EnergyRatioFromVertexY", ";Energy Ratio;", 100, 0, 1);
  EnergyRatioFromVertexZ =
      tfs->make<TH1D>("EnergyRatioFromVertexZ", ";Energy Ratio;", 100, 0, 1);

  TimingCorrelation = tfs->make<TH1D>(
      "TimingCorrelation", ";Timing Correlation Factor;", 100, -1.5, 1.5);

  DistanceFromReconVertexToTrueVertex =
      tfs->make<TH1D>("DistanceFromReconVertexToTrueVertex",
                      ";Vertex Error (cm);", 100, 0, 200);
  ZErrorFromReconVertexToTrueVertex =
      tfs->make<TH1D>("ZErrorFromReconVertexToTrueVertex",
                      ";Vertex Z Error (cm);", 100, 0, 200);
  OrientationFromReconVertexToTrueVertex =
      tfs->make<TH2D>("OrientationFromReconVertexToTrueVertex", ";#phi;#theta",
                      50, -TMath::Pi(), TMath::Pi(), 50, 0, TMath::Pi());
  PolarAngleFromReconVertexToTrueVertex = tfs->make<TH1D>(
      "PolarAngleFromReconVertexToTrueVertex", ";#phi;", 50, -1, 7);
  AzimuthalAngleFromReconVertexToTrueVertex = tfs->make<TH1D>(
      "AzimuthalAngleFromReconVertexToTrueVertex", ";#theta;", 50, -2, 2);

  DistanceFromElasticArmVertexToTrueVertex =
      tfs->make<TH1D>("DistanceFromElasticArmVertexToTrueVertex",
                      ";EArm Vertex Error (cm);", 100, 0, 200);
  ZErrorFromElasticArmVertexToTrueVertex =
      tfs->make<TH1D>("ZErrorFromElasticArmVertexToTrueVertex",
                      ";Vertex Z Error (cm);", 100, 0, 200);
  OrientationFromElasticArmVertexToTrueVertex = tfs->make<TH2D>(
      "OrientationFromElasticArmVertexToTrueVertex", ";#phi;#theta", 50,
      -TMath::Pi(), TMath::Pi(), 50, 0, TMath::Pi());
  PolarAngleFromElasticArmVertexToTrueVertex = tfs->make<TH1D>(
      "PolarAngleFromElasticArmVertexToTrueVertex", ";#phi;", 50, -1, 7);
  AzimuthalAngleFromElasticArmVertexToTrueVertex = tfs->make<TH1D>(
      "AzimuthalAngleFromElasticArmVertexToTrueVertex", ";#theta;", 50, -2, 2);
}

TGraph *nnbar::RecoAnalysis::PlotVectorPoint(std::vector<point2d_t> pts) {
  double *x = new double[pts.size()];
  double *y = new double[pts.size()];

  size_t i = 0;
  for (auto &pt : pts) {
    *(x + i) = pt.first;
    *(y + i) = pt.second;
    i++;
  }

  TGraph *graph = new TGraph(pts.size(), x, y);
  return graph;
}

double nnbar::RecoAnalysis::slopeAngleOfLineConnect2Points(point2d_t pA,
                                                           point2d_t pB) {
  double xa = pA.first;
  double ya = pA.second;
  double xb = pB.first;
  double yb = pB.second;

  return TMath::ATan((ya - yb) / (xa - xb));
}

void nnbar::RecoAnalysis::analyze(art::Event const &e) {

  _ExposureTimeInUSec += 50;

  art::Handle<std::vector<simb::MCTruth>> MCTruthHandle;
  std::vector<art::Ptr<simb::MCTruth>> mctruths;
  if (e.getByLabel("generator", MCTruthHandle))
    art::fill_ptr_vector(mctruths, MCTruthHandle);

  double trueVertexX = 0, trueVertexY = 0, trueVertexZ = 0, trueVertexT = 0;
  for (std::vector<art::Ptr<simb::MCTruth>>::iterator mctruthItr =
           mctruths.begin();
       mctruthItr != mctruths.end(); ++mctruthItr) {
    auto mctruth = *mctruthItr;
    auto particle = (*mctruth).GetParticle(0);
    trueVertexX = particle.Vx();
    trueVertexY = particle.Vy();
    trueVertexZ = particle.Vz();
    trueVertexT = particle.T();
  }

  art::Handle<std::vector<rb::Cluster>> clusterHandle;
  std::vector<art::Ptr<rb::Cluster>> clusters;
  if (e.getByLabel(_InputShwrLabel, clusterHandle))
    art::fill_ptr_vector(clusters, clusterHandle);

  art::FindManyP<rb::Vertex> assnClusterVert(clusterHandle, e, "elasticarmshs");

  //  double elasticArmVertT;
  double elasticArmVertX;
  double elasticArmVertY;
  double elasticArmVertZ;
  unsigned int clusterIdx = 0;
  for (std::vector<art::Ptr<rb::Cluster>>::iterator clusterItr =
           clusters.begin();
       clusterItr != clusters.end(); ++clusterItr) {
    auto cluster = *clusterItr;

    art::PtrVector<rb::CellHit> allCells = (*cluster).AllCells();

    // Conventional Vertex Reconstruction
    //    elasticArmVertT = -1E9;
    elasticArmVertX = -1E9;
    elasticArmVertY = -1E9;
    elasticArmVertZ = -1E9;
    if (assnClusterVert.isValid()) {
      std::vector<art::Ptr<rb::Vertex>> vert = assnClusterVert.at(clusterIdx);
      if (vert.size() != 1)
        continue;
      //      elasticArmVertT = vert.at(0)->GetT();
      elasticArmVertX = vert.at(0)->GetX();
      elasticArmVertY = vert.at(0)->GetY();
      elasticArmVertZ = vert.at(0)->GetZ();
    }

    //--------------------------------------------------------------------------

    double diff_true_elasticArm_vertex_x = elasticArmVertX - trueVertexX;
    double diff_true_elasticArm_vertex_y = elasticArmVertY - trueVertexY;
    double diff_true_elasticArm_vertex_z = elasticArmVertZ - trueVertexZ;
    double distFromElasticArmVertToTrueVert =
        TMath::Sqrt(TMath::Power(diff_true_elasticArm_vertex_x, 2) +
                    TMath::Power(diff_true_elasticArm_vertex_y, 2) +
                    TMath::Power(diff_true_elasticArm_vertex_z, 2));
    double thetaOrientationFromElasticArmVertToTrueVert = TMath::ATan(
        diff_true_elasticArm_vertex_y / diff_true_elasticArm_vertex_x);
    double phiOrientationFromElasticArmVertToTrueVert = TMath::ACos(
        diff_true_elasticArm_vertex_z / distFromElasticArmVertToTrueVert);

    DistanceFromElasticArmVertexToTrueVertex->Fill(
        distFromElasticArmVertToTrueVert);
    OrientationFromElasticArmVertexToTrueVertex->Fill(
        phiOrientationFromElasticArmVertToTrueVert,
        thetaOrientationFromElasticArmVertToTrueVert);
    PolarAngleFromElasticArmVertexToTrueVertex->Fill(
        phiOrientationFromElasticArmVertToTrueVert);
    AzimuthalAngleFromElasticArmVertexToTrueVertex->Fill(
        thetaOrientationFromElasticArmVertToTrueVert);

    // CONSIDERING THE ENERGY DISTRIBUTION TOPOLOGY
    double totalE = (*cluster).TotalGeV();
    double totalEX = 0.;
    double totalEY = 0.;
    double totalEnergyOfHitsAboveVertexX = 0.;
    double totalEnergyOfHitsAboveVertexY = 0.;
    double totalEnergyOfHitsAboveVertexZ = 0.;

    // Consider timing correlation
    std::multimap<double, double> distanceToVertex_relTimingToVertex_Map;
    unsigned int hitIdx = 0;

    // Compare distance between true vertex and reconstructed vertex
    std::multimap<double, double> HitsZ, HitsX, HitsY;

    for (art::PtrVector<rb::CellHit>::iterator itrCell = allCells.begin();
         itrCell != allCells.end(); itrCell++) {
      if ((*itrCell)->View() == geo::View_t::kXorY)
        continue;

      bool cIsYView = (*itrCell)->View() == geo::View_t::kY ? true : false;
      rb::RecoHit recoHit = (*cluster).RecoHit(*itrCell);
      if (!recoHit.IsCalibrated())
        continue;

      double cX = recoHit.X();
      double cY = recoHit.Y();
      double cZ = recoHit.Z();
      //      double cT = recoHit.T();
      double cE = recoHit.GeV();

      HitsZ.insert(std::pair<double, double>(cZ, cE));
      if (!cIsYView) {
        HitsX.insert(std::pair<double, double>(cX, cE));
        totalEX += cE;
      } else {
        HitsY.insert(std::pair<double, double>(cY, cE));
        totalEY += cE;
      }

      // // Consider timing correlation
      // // Timing correlation with true vertex
      // double distanceToVertex = TMath::Sqrt(TMath::Power(cX - trueVertexX, 2)
      // +
      //                                       TMath::Power(cY - trueVertexY, 2)
      //                                       + TMath::Power(cZ - trueVertexY,
      //                                       2));
      // double relTimingToVertex = cT - trueVertexT;
      // distanceToVertex_relTimingToVertex_Map.insert(
      //     std::pair<double, double>(distanceToVertex, relTimingToVertex));

      if (!cIsYView) {
        if (cX > trueVertexX)
          totalEnergyOfHitsAboveVertexX += cE;
      } else {
        if (cY > trueVertexY)
          totalEnergyOfHitsAboveVertexY += cE;
      }
      if (cZ > trueVertexZ)
        totalEnergyOfHitsAboveVertexZ += cE;

      // COUNTING THE NUMBER OF PRONGS
      if (!cIsYView) {
        HitsXView.push_back(std::make_tuple(cZ, cX, cE));
      } else {
        HitsYView.push_back(std::make_tuple(cZ, cY, cE));
      }

      hitIdx++;
    }

    // Compare distance between true vertex and reconstructed vertex
    double lowerBoundOfMedianZVertex = 0.;
    double upperBoundOfMedianZVertex = 0.;
    double halfOfTotalE = 0.;
    for (auto hit = HitsZ.begin(); hit != HitsZ.end(); hit++) {
      halfOfTotalE += (*hit).second;
      if (halfOfTotalE < totalE / 2) {
        lowerBoundOfMedianZVertex = (*hit).first;
        upperBoundOfMedianZVertex = (*(hit++)).first;
        hit--;
      }
    }
    double reconVertexZ =
        (lowerBoundOfMedianZVertex + upperBoundOfMedianZVertex) / 2.;

    double lowerBoundOfMedianXVertex = 0.;
    double upperBoundOfMedianXVertex = 0.;
    double halfOfTotalEX = 0.;
    for (auto hit = HitsX.begin(); hit != HitsX.end(); hit++) {
      halfOfTotalEX += (*hit).second;
      if (halfOfTotalEX < totalEX / 2) {
        lowerBoundOfMedianXVertex = (*hit).first;
        upperBoundOfMedianXVertex = (*(hit++)).first;
        hit--;
      }
    }
    double reconVertexX =
        (lowerBoundOfMedianXVertex + upperBoundOfMedianXVertex) / 2.;

    double lowerBoundOfMedianYVertex = 0.;
    double upperBoundOfMedianYVertex = 0.;
    double halfOfTotalEY = 0.;
    for (auto hit = HitsY.begin(); hit != HitsY.end(); hit++) {
      halfOfTotalEY += (*hit).second;
      if (halfOfTotalEY < totalEY / 2) {
        lowerBoundOfMedianYVertex = (*hit).first;
        upperBoundOfMedianYVertex = (*(hit++)).first;
        hit--;
      }
    }
    double reconVertexY =
        (lowerBoundOfMedianYVertex + upperBoundOfMedianYVertex) / 2.;

    for (art::PtrVector<rb::CellHit>::iterator itrCell = allCells.begin();
         itrCell != allCells.end(); itrCell++) {
      if ((*itrCell)->View() == geo::View_t::kXorY)
        continue;

      rb::RecoHit recoHit = (*cluster).RecoHit(*itrCell);
      if (!recoHit.IsCalibrated())
        continue;

      double cX = recoHit.X();
      double cY = recoHit.Y();
      double cZ = recoHit.Z();
      double cT = recoHit.T();

      // Consider timing correlation
      // Timing correlation with recon vertex
      double distanceToVertex = TMath::Sqrt(TMath::Power(cX - reconVertexX, 2) +
                                            TMath::Power(cY - reconVertexY, 2) +
                                            TMath::Power(cZ - reconVertexZ, 2));
      double relTimingToVertex = cT;
      distanceToVertex_relTimingToVertex_Map.insert(
          std::pair<double, double>(distanceToVertex, relTimingToVertex));
    }

    //--------------------------------------------------------------------------
    // Consider timing correlation
    // Recon Vertex
    const unsigned int mapSize = distanceToVertex_relTimingToVertex_Map.size();
    unsigned int countElem = 0;
    double *distanceToVertex_array = new double[10];
    double *relTimingToVertex_array = new double[10];
    for (auto &elem : distanceToVertex_relTimingToVertex_Map) {
      if (4 < countElem && countElem < mapSize - 5)
        continue;
      *(distanceToVertex_array + countElem) = elem.first;
      *(relTimingToVertex_array + countElem) = elem.second;
      countElem++;
    }
    // The timing of reconstructed vertex is calculated as the mean timing of
    // the 5 hits nearest to the vertex.
    // So this loop find that average timing.
    double reconVertexT = 0.;
    for (unsigned int j = 0; j < 5; j++) {
      reconVertexT += *(relTimingToVertex_array + j);
    }
    reconVertexT = reconVertexT / 5;
    for (unsigned int j = 0; j < 5; j++) {
      *(relTimingToVertex_array + countElem) =
          *(relTimingToVertex_array + countElem) - reconVertexT;
    }
    TGraph *timingCorrelation =
        new TGraph(10, distanceToVertex_array, relTimingToVertex_array);
    double correlationFactor = timingCorrelation->GetCorrelationFactor();
    TimingCorrelation->Fill(correlationFactor);

    //--------------------------------------------------------------------------

    double diff_true_recon_vertex_x = reconVertexX - trueVertexX;
    double diff_true_recon_vertex_y = reconVertexY - trueVertexY;
    double diff_true_recon_vertex_z = reconVertexZ - trueVertexZ;
    double distFromReconVertToTrueVert =
        TMath::Sqrt(TMath::Power(diff_true_recon_vertex_x, 2) +
                    TMath::Power(diff_true_recon_vertex_y, 2) +
                    TMath::Power(diff_true_recon_vertex_z, 2));
    double thetaOrientationFromReconVertToTrueVert =
        TMath::ATan(diff_true_recon_vertex_y / diff_true_recon_vertex_x);
    double phiOrientationFromReconVertToTrueVert =
        TMath::ACos(diff_true_recon_vertex_z / distFromReconVertToTrueVert);
    DistanceFromReconVertexToTrueVertex->Fill(distFromReconVertToTrueVert);
    OrientationFromReconVertexToTrueVertex->Fill(
        phiOrientationFromReconVertToTrueVert,
        thetaOrientationFromReconVertToTrueVert);
    PolarAngleFromReconVertexToTrueVertex->Fill(
        phiOrientationFromReconVertToTrueVert);
    AzimuthalAngleFromReconVertexToTrueVertex->Fill(
        thetaOrientationFromReconVertToTrueVert);

    // // Consider timing correlation
    // // True Vertex
    // const unsigned int mapSize =
    // distanceToVertex_relTimingToVertex_Map.size(); unsigned int countElem =
    // 0; double *distanceToVertex_array = new double[10]; double
    // *relTimingToVertex_array = new double[10]; for (auto &elem :
    // distanceToVertex_relTimingToVertex_Map) {
    //   if (4 < countElem && countElem < mapSize - 5)
    //     continue;
    //   *(distanceToVertex_array + countElem) = elem.first;
    //   *(relTimingToVertex_array + countElem) = elem.second;
    //   countElem++;
    // }
    // TGraph *timingCorrelation =
    //     new TGraph(10, distanceToVertex_array, relTimingToVertex_array);
    // double correlationFactor = timingCorrelation->GetCorrelationFactor();
    // TimingCorrelation->Fill(correlationFactor);

    // CONSIDERING THE ENERGY DISTRIBUTION TOPOLOGY
    EnergyRatioFromVertexX->Fill(totalEnergyOfHitsAboveVertexX / totalE);
    EnergyRatioFromVertexY->Fill(totalEnergyOfHitsAboveVertexY / totalE);
    EnergyRatioFromVertexZ->Fill(totalEnergyOfHitsAboveVertexZ / totalE);

    // Filling analysis tree
    _tRunNumber = e.run();
    _tSubrunNumber = e.subRun();
    _tEventNumber = e.id().event();
    _tSliceNumber = std::distance(clusters.begin(), clusterItr);
    _tSampleType = _SampleType;
    double totVisEX =
        nnbar::NNbarUtilities::totalVisibleEnergyInView(*cluster, geo::kX);
    double totVisEY =
        nnbar::NNbarUtilities::totalVisibleEnergyInView(*cluster, geo::kY);
    _tTotalVisibleEnergy = totVisEX + totVisEY;
    if (_tTotalVisibleEnergy > 0.8) {
      ZErrorFromReconVertexToTrueVertex->Fill(
          TMath::Abs(diff_true_recon_vertex_z));
      ZErrorFromElasticArmVertexToTrueVertex->Fill(
          TMath::Abs(diff_true_elasticArm_vertex_z));
    }

    _tTemporalClusterExpand = nnbar::NNbarUtilities::temporalExpand(*cluster);
    _tTotalHitCountInXView =
        nnbar::NNbarUtilities::hitCountInView(*cluster, geo::kX);
    _tTotalHitCountInYView =
        nnbar::NNbarUtilities::hitCountInView(*cluster, geo::kY);
    _tHitCountXYDifference = _tTotalHitCountInYView - _tTotalHitCountInXView;
    _tAverageEnergyPerHitXView = totVisEX / _tTotalHitCountInXView;
    _tAverageEnergyPerHitYView = totVisEY / _tTotalHitCountInYView;
    _tHitVsVertexTimingCorrelation = correlationFactor;
    _tReconVertexX = reconVertexX;
    _tReconVertexY = reconVertexY;
    _tReconVertexZ = reconVertexZ;
    _tReconVertexT = reconVertexT;
    if (_tSampleType != 0) {
      _tTrueVertexX = trueVertexX;
      _tTrueVertexY = trueVertexY;
      _tTrueVertexZ = trueVertexZ;
      _tTrueVertexT = trueVertexT;
    } else {
      _tTrueVertexX = -1.0E9;
      _tTrueVertexY = -1.0E9;
      _tTrueVertexZ = -1.0E9;
      _tTrueVertexT = -1.0E9;
    }
    _VariableTree->Fill();

    clusterIdx++;
  }

  HitsXView.clear();
  HitsYView.clear();
}

void nnbar::RecoAnalysis::evd(art::Event const &e, unsigned int clusterIdx) {

  // TGRAPH XVIEW
  const unsigned int nhit_xview = HitsXView.size();
  double *x_xv = new double[nhit_xview];
  double *y_xv = new double[nhit_xview];
  unsigned int countHitXView = 0;
  for (auto &hit : HitsXView) {
    *(x_xv + countHitXView) = std::get<0>(hit);
    *(y_xv + countHitXView) = std::get<1>(hit);
    countHitXView++;
  }
  TGraph *graph_xv = new TGraph(nhit_xview, x_xv, y_xv);

  // TGRAPH YVIEW
  const unsigned int nhit_yview = HitsYView.size();
  double *x_yv = new double[nhit_yview];
  double *y_yv = new double[nhit_yview];
  unsigned int countHitYView = 0;
  for (auto &hit : HitsYView) {
    *(x_yv + countHitYView) = std::get<0>(hit);
    *(y_yv + countHitYView) = std::get<1>(hit);
    countHitYView++;
  }
  TGraph *graph_yv = new TGraph(nhit_yview, x_yv, y_yv);

  // DRAW
  TCanvas *c =
      new TCanvas(Form("c_%i_%i", e.event(), clusterIdx), "", 1200, 1200);
  c->Divide(1, 2);
  c->cd(1);
  graph_xv->SetTitle("");
  graph_xv->GetXaxis()->SetLimits(0, 6400);
  graph_xv->GetYaxis()->SetRangeUser(-800, 800);
  graph_xv->GetXaxis()->SetTitle("z");
  graph_xv->GetYaxis()->SetTitle("x");
  graph_xv->GetXaxis()->CenterTitle();
  graph_xv->GetYaxis()->CenterTitle();
  graph_xv->SetMarkerColor(kRed);
  graph_xv->SetMarkerStyle(20);
  graph_xv->SetMarkerSize(1);
  graph_xv->Draw("AP");
  c->cd(2);
  graph_yv->SetTitle("");
  graph_yv->GetXaxis()->SetLimits(0, 6400);
  graph_yv->GetYaxis()->SetRangeUser(-800, 800);
  graph_yv->GetXaxis()->SetTitle("z");
  graph_yv->GetYaxis()->SetTitle("y");
  graph_yv->GetXaxis()->CenterTitle();
  graph_yv->GetYaxis()->CenterTitle();
  graph_yv->SetMarkerColor(kRed);
  graph_yv->SetMarkerStyle(20);
  graph_yv->SetMarkerSize(1);
  graph_yv->Draw("AP");
  c->SaveAs(Form("Plots/EVD_%i_%i.pdf", e.event(), clusterIdx));
}

void nnbar::RecoAnalysis::endJob() {
  std::cout << "Exposure Time In USec: " << _ExposureTimeInUSec << std::endl;
  // gStyle->SetOptStat(0);

  //----------------------------------------------------------------------------

  TCanvas *canvas = new TCanvas("canvas", "Canvas", 2400, 800);
  canvas->Divide(3, 1);
  canvas->cd(1);
  EnergyRatioFromVertexX->SetTitle("x");
  EnergyRatioFromVertexX->Draw();
  canvas->cd(2);
  EnergyRatioFromVertexY->SetTitle("y");
  EnergyRatioFromVertexY->Draw();
  canvas->cd(3);
  EnergyRatioFromVertexZ->SetTitle("z");
  EnergyRatioFromVertexZ->Draw();
  canvas->SaveAs(Form("EnergyVsTopology.pdf"));

  TCanvas *canvas_timing =
      new TCanvas("canvas_timing", "Canvas Timing", 1200, 1200);
  TimingCorrelation->GetXaxis()->SetTitle("Timing Correlation Factor");
  TimingCorrelation->GetXaxis()->CenterTitle();
  TimingCorrelation->GetYaxis()->SetTitle("");
  TimingCorrelation->GetYaxis()->CenterTitle();
  TimingCorrelation->Draw("HIST");
  canvas_timing->SaveAs(Form("TimingCorrelation.pdf"));

  TCanvas *canvas_vertex_error =
      new TCanvas("canvas_vertex_error", "Canvas Vertex Error", 2400, 800);
  canvas_vertex_error->Divide(3, 1);
  canvas_vertex_error->cd(1);
  DistanceFromReconVertexToTrueVertex->GetXaxis()->SetTitle(
      "#deltar_{equale-(v)} (cm)");
  DistanceFromReconVertexToTrueVertex->GetXaxis()->CenterTitle();
  DistanceFromReconVertexToTrueVertex->GetYaxis()->SetTitle("");
  DistanceFromReconVertexToTrueVertex->GetYaxis()->CenterTitle();
  DistanceFromReconVertexToTrueVertex->Draw("HIST");
  canvas_vertex_error->cd(2);
  PolarAngleFromReconVertexToTrueVertex->GetXaxis()->SetTitle(
      "#phi_{equale-(v)} (rad)");
  PolarAngleFromReconVertexToTrueVertex->GetXaxis()->CenterTitle();
  PolarAngleFromReconVertexToTrueVertex->GetYaxis()->SetTitle("");
  PolarAngleFromReconVertexToTrueVertex->GetYaxis()->CenterTitle();
  PolarAngleFromReconVertexToTrueVertex->Draw("HIST");
  canvas_vertex_error->cd(3);
  AzimuthalAngleFromReconVertexToTrueVertex->GetXaxis()->SetTitle(
      "#theta_{equale-(v)} (rad)");
  AzimuthalAngleFromReconVertexToTrueVertex->GetXaxis()->CenterTitle();
  AzimuthalAngleFromReconVertexToTrueVertex->GetYaxis()->SetTitle("");
  AzimuthalAngleFromReconVertexToTrueVertex->GetYaxis()->CenterTitle();
  AzimuthalAngleFromReconVertexToTrueVertex->Draw("HIST");
  canvas_vertex_error->SaveAs(Form("VertexError.pdf"));

  TCanvas *canvas_earm_vertex_error = new TCanvas(
      "canvas_earm_vertex_error", "Canvas EArm Vertex Error", 2400, 800);
  canvas_earm_vertex_error->Divide(3, 1);
  canvas_earm_vertex_error->cd(1);
  DistanceFromElasticArmVertexToTrueVertex->GetXaxis()->SetTitle(
      "#deltar_{earm-(v)} (cm)");
  DistanceFromElasticArmVertexToTrueVertex->GetXaxis()->CenterTitle();
  DistanceFromElasticArmVertexToTrueVertex->GetYaxis()->SetTitle("");
  DistanceFromElasticArmVertexToTrueVertex->GetYaxis()->CenterTitle();
  DistanceFromElasticArmVertexToTrueVertex->Draw("HIST");
  canvas_earm_vertex_error->cd(2);
  PolarAngleFromElasticArmVertexToTrueVertex->GetXaxis()->SetTitle(
      "#phi_{earm-(v)} (rad)");
  PolarAngleFromElasticArmVertexToTrueVertex->GetXaxis()->CenterTitle();
  PolarAngleFromElasticArmVertexToTrueVertex->GetYaxis()->SetTitle("");
  PolarAngleFromElasticArmVertexToTrueVertex->GetYaxis()->CenterTitle();
  PolarAngleFromElasticArmVertexToTrueVertex->Draw("HIST");
  canvas_earm_vertex_error->cd(3);
  AzimuthalAngleFromElasticArmVertexToTrueVertex->GetXaxis()->SetTitle(
      "#theta_{earm-(v)} (rad)");
  AzimuthalAngleFromElasticArmVertexToTrueVertex->GetXaxis()->CenterTitle();
  AzimuthalAngleFromElasticArmVertexToTrueVertex->GetYaxis()->SetTitle("");
  AzimuthalAngleFromElasticArmVertexToTrueVertex->GetYaxis()->CenterTitle();
  AzimuthalAngleFromElasticArmVertexToTrueVertex->Draw("HIST");
  canvas_earm_vertex_error->SaveAs(Form("EArmVertexError.pdf"));

  TCanvas *canvas_zvertex_error =
      new TCanvas("canvas_zvertex_error", "Canvas Z Vertex Error", 2400, 1200);
  canvas_zvertex_error->Divide(2, 1);
  canvas_zvertex_error->cd(1);
  ZErrorFromElasticArmVertexToTrueVertex->GetXaxis()->SetTitle(
      "#deltaz_{earm-(v)} (cm)");
  ZErrorFromElasticArmVertexToTrueVertex->GetXaxis()->CenterTitle();
  ZErrorFromElasticArmVertexToTrueVertex->GetYaxis()->SetTitle("");
  ZErrorFromElasticArmVertexToTrueVertex->GetYaxis()->CenterTitle();
  ZErrorFromElasticArmVertexToTrueVertex->Draw("HIST");
  canvas_zvertex_error->cd(2);
  ZErrorFromReconVertexToTrueVertex->GetXaxis()->SetTitle(
      "#deltaz_{equale-(v)} (cm)");
  ZErrorFromReconVertexToTrueVertex->GetXaxis()->CenterTitle();
  ZErrorFromReconVertexToTrueVertex->GetYaxis()->SetTitle("");
  ZErrorFromReconVertexToTrueVertex->GetYaxis()->CenterTitle();
  ZErrorFromReconVertexToTrueVertex->Draw("HIST");
  canvas_zvertex_error->SaveAs(Form("ZVertexError.pdf"));
}

DEFINE_ART_MODULE(nnbar::RecoAnalysis)