SNSlicerAna_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       SNSlicerAna
// Module Type: analyzer
// File:        SNSlicerAna_module.cc
//
// Generated at Tue May 23 12:56:03 2017 by Justin Vasel using artmod
// from cetpkgsupport v1_10_02.
////////////////////////////////////////////////////////////////////////

// ROOT includes
#include <TTree.h>
#include <TVector3.h>

// Framework includes
#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/Services/Optional/TFileService.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// NOvASoft includes
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/PlaneGeo.h"
#include "MCCheater/BackTracker.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "Simulation/Particle.h"


namespace sn {
  class SNSlicerAna;
}

class sn::SNSlicerAna : public art::EDAnalyzer {
public:
  explicit SNSlicerAna(fhicl::ParameterSet const & p);

  SNSlicerAna(SNSlicerAna const &) = delete;
  SNSlicerAna(SNSlicerAna &&) = delete;
  SNSlicerAna & operator = (SNSlicerAna const &) = delete;
  SNSlicerAna & operator = (SNSlicerAna &&) = delete;

  // Required functions.
  void analyze(art::Event const & e) override;

  // Selected optional functions.
  void beginJob() override;

private:
  art::ServiceHandle<art::TFileService> tfs;
  art::ServiceHandle<geo::Geometry>     fGeom;


  std::string fClusterLabel;
  std::string fCellHitLabel;

  TTree* fClusterTree;
  TTree* fEventTree;
  TTree* fHitTree;
  TTree* fHitPairTree;

  // ART stuff
  int fEventNumber;              /// ART event number
  int fClusterId;                /// Unique ID for each cluster in the event
  double fClusterTNS;            /// Min time of the cluster

  // Hit metrics
  int fNHitsMC;                  /// Number of MC hits in cluster
  int fNHitsBG;                  /// Number of non-MC hits in cluster
  int fAdcMc;                    /// ADC for MC hits
  int fAdcBg;                    /// ADC for BG hits
  int fNFEB;                     /// Number of distinct FEBs in the cluster

  // Cluster metrics
  double fSumADC;                /// Total cluster ADC
  int fSpanADC;                  /// Difference between max and min ADCs
  double fMeanADC;               /// Mean cluster ADC
  double fStdADC;                /// Std. Dev. of cluster ADCs
  double fSumMeV;                /// Cluster total energy (MeV)
  double fSpanTNS;               /// Cluster time span (ns)
  int fSpanPlane;                /// Cluster plane span (planes)
  int fSpanCellX;                /// Cluster cell span in X-view (cells)
  int fSpanCellY;                /// Cluster cell span in Y-view (cells)
  int fSpanViews;                /// Cluster view span (1 or 2)
  float fPairwiseADC;            /// Average of hit-pairwise ADC differences in a cluster
  float fPairwiseTNS;            /// Average of hit-pairwise TNS differences in a cluster

  // Means
  double fMeanTNS;               /// Cluster time mean (ns)
  double fMeanX;                 /// Cluster X mean
  double fMeanY;                 /// Cluster Y mean
  double fMeanZ;                 /// Cluster Z mean

  float fClusterEfficiency;      /// Number of MC hits (in all clusters / total in event)

  float fPurity;                 /// Cluster Purity (Nsig / Ntotal)
  int fAdcDiff;                  /// Difference in ADC between pair-wise hits
  int fPlaneDiff;
  int fCellDiff;
  float fXdiff;
};


sn::SNSlicerAna::SNSlicerAna(fhicl::ParameterSet const & p):
  EDAnalyzer(p),
  fClusterLabel(p.get<std::string>("ClusterLabel")),
  fCellHitLabel(p.get<std::string>("CellHitLabel")),
  fEventNumber(0),
  fClusterId(0),
  fClusterTNS(0),
  fNHitsMC(0),
  fNHitsBG(0),
  fAdcMc(0),
  fAdcBg(0),
  fNFEB(0),
  fSumADC(0),
  fSpanADC(0),
  fMeanADC(0),
  fStdADC(0),
  fSumMeV(0),
  fSpanTNS(0),
  fSpanPlane(0),
  fSpanCellX(0),
  fSpanCellY(0),
  fSpanViews(0),
  fPairwiseADC(0),
  fPairwiseTNS(0),
  fMeanTNS(0),
  fMeanX(0),
  fMeanY(0),
  fMeanZ(0),
  fClusterEfficiency(0),
  fPurity(0),
  fAdcDiff(0),
  fPlaneDiff(0),
  fCellDiff(0),
  fXdiff(0)
{
}


void sn::SNSlicerAna::analyze(art::Event const & e)
{
  fEventNumber = e.id().event();

  // Get the clusters from the ART event
  art::Handle<std::vector<rb::Cluster>> clusters;
  e.getByLabel(fClusterLabel, clusters);

  int fClusterHitsMC = 0;

  for (size_t iCluster = 0; iCluster < clusters->size(); ++iCluster) {
    rb::Cluster aCluster = clusters->at(iCluster);
    fClusterId = iCluster;

    if (aCluster.NCell()==0) continue;

    fClusterTNS = aCluster.MinTNS();

    fSumADC    = aCluster.TotalADC();
    fSumMeV    = aCluster.TotalGeV() * 1.0e3;
    fMeanADC   = aCluster.TotalADC() / (float)aCluster.NCell();
    fSpanTNS   = aCluster.ExtentTNS();
    fSpanPlane = aCluster.ExtentPlane();
    fSpanCellX = (aCluster.NXCell()) ? aCluster.ExtentCell(geo::kX) : 0;
    fSpanCellY = (aCluster.NYCell()) ? aCluster.ExtentCell(geo::kY) : 0;

    fSpanViews = (fSpanCellX>0 && fSpanCellY>0) ? 2 : 1;

    fMeanTNS = aCluster.MeanTNS();
    fMeanX   = aCluster.MeanX();
    fMeanY   = aCluster.MeanY();
    fMeanZ   = aCluster.MeanZ();


    // Loop over the contained cell hits
    std::map<int, int> mapInteractionCounts;
    fNHitsMC = 0;
    fNHitsBG = 0;
    std::vector<unsigned int> fFEBs;

    int fMinADC=0;
    int fMaxADC=0;
    float fVarADC=0;
    art::PtrVector<rb::CellHit> allHits = aCluster.AllCells();
    for (art::Ptr<rb::CellHit> aHit : allHits) {
      fAdcMc = 0;
      fAdcBg = 0;
      (aHit->IsMC()) ? fNHitsMC++ : fNHitsBG++;
      (aHit->IsMC()) ? fAdcMc = aHit->ADC() : fAdcBg = aHit->ADC();

      if (fMinADC==0 || aHit->ADC() < fMinADC) fMinADC = aHit->ADC();
      if (fMaxADC==0 || aHit->ADC() > fMaxADC) fMaxADC = aHit->ADC();

      fVarADC += pow(aHit->ADC() - fMeanADC, 2) / (float)aCluster.NCell();

      unsigned int fFEB = aHit->DaqChannel() & 0xFFFFFF00;
      if (std::find(fFEBs.begin(), fFEBs.end(), fFEB) == fFEBs.end()) fFEBs.push_back(fFEB);

      // fHitTree->Fill();
    }

    fStdADC = sqrt(fVarADC);

    int n = 0;
    float totalPairwiseADC = 0;
    float totalPairwiseTNS = 0;
    for (size_t i=0; i<allHits.size()-1; ++i) {
      for (size_t j=i+1; j<allHits.size(); ++j) {
        ++n;

        totalPairwiseADC = std::abs(allHits[i]->ADC() - allHits[j]->ADC());
        totalPairwiseTNS = std::abs(allHits[i]->TNS() - allHits[j]->TNS());
        fAdcDiff = std::abs(allHits[i]->ADC() - allHits[j]->ADC());
        fPlaneDiff = std::abs(allHits[i]->Plane() - allHits[j]->Plane());
        fCellDiff = std::abs(allHits[i]->Cell() - allHits[j]->Cell());
        fPurity = fNHitsMC / (fNHitsMC + fNHitsBG);

        TVector3 x1, x2;
        fGeom->Plane(allHits[i]->Plane())->Cell(allHits[i]->Cell())->GetCenter(x1);
        fGeom->Plane(allHits[j]->Plane())->Cell(allHits[j]->Cell())->GetCenter(x2);

        fXdiff = (x1 - x2).Mag();

        fHitPairTree->Fill();
      }
    }
    fPairwiseADC = totalPairwiseADC / (float)n;
    fPairwiseTNS = totalPairwiseTNS / (float)n;

    fSpanADC = fMaxADC - fMinADC;

    fNFEB = fFEBs.size();
    fClusterHitsMC += fNHitsMC;

    fClusterTree->Fill();
  }

  // Get the hits from the ART event
  art::Handle<std::vector<rb::CellHit>> hits;
  e.getByLabel(fCellHitLabel, hits);

  int fEventHitsMC = 0;
  for (rb::CellHit aHit : *hits) {
    if (aHit.IsMC()) ++fEventHitsMC;
  }

  fClusterEfficiency = (float)fClusterHitsMC / (float)fEventHitsMC;

  fEventTree->Fill();
}


void sn::SNSlicerAna::beginJob()
{
  fClusterTree = tfs->make<TTree>("ClusterTree", "ClusterTree");
  fClusterTree->Branch("EventNumber",            &fEventNumber);
  fClusterTree->Branch("ClusterID",              &fClusterId);
  fClusterTree->Branch("ClusterTNS",             &fClusterTNS);
  fClusterTree->Branch("NHitsMC",                &fNHitsMC);
  fClusterTree->Branch("NHitsBG",                &fNHitsBG);
  fClusterTree->Branch("NFEB",                   &fNFEB);
  fClusterTree->Branch("SumADC",                 &fSumADC);
  fClusterTree->Branch("SpanADC",                &fSpanADC);
  fClusterTree->Branch("MeanADC",                &fMeanADC);
  fClusterTree->Branch("StdADC",                 &fStdADC);
  fClusterTree->Branch("SumMeV",                 &fSumMeV);
  fClusterTree->Branch("PairwiseADC",            &fPairwiseADC);
  fClusterTree->Branch("PairwiseTNS",            &fPairwiseTNS);
  fClusterTree->Branch("SpanTNS",                &fSpanTNS);
  fClusterTree->Branch("SpanPlane",              &fSpanPlane);
  fClusterTree->Branch("SpanCellX",              &fSpanCellX);
  fClusterTree->Branch("SpanCellY",              &fSpanCellY);
  fClusterTree->Branch("SpanViews",              &fSpanViews);
  fClusterTree->Branch("MeanTNS",                &fMeanTNS);
  fClusterTree->Branch("MeanX",                  &fMeanX);
  fClusterTree->Branch("MeanY",                  &fMeanY);
  fClusterTree->Branch("MeanZ",                  &fMeanZ);

  fEventTree = tfs->make<TTree>("EventTree", "EventTree");
  fEventTree->Branch("EventNumber",              &fEventNumber);
  fEventTree->Branch("ClusterEfficiency",        &fClusterEfficiency);

  fHitTree = tfs->make<TTree>("HitTree", "HitTree");
  fHitTree->Branch("EventNumber",                &fEventNumber);
  fHitTree->Branch("AdcMc",                      &fAdcMc);
  fHitTree->Branch("AdcBg"  ,                    &fAdcBg);

  fHitPairTree = tfs->make<TTree>("HitPairTree", "HitPairTree");
  fHitPairTree->Branch("EventNumber",            &fEventNumber);
  fHitPairTree->Branch("ClusterID",              &fClusterId);
  fHitPairTree->Branch("ClusterTNS",             &fClusterTNS);
  fHitPairTree->Branch("Purity",                 &fPurity);
  fHitPairTree->Branch("PlaneDiff",              &fPlaneDiff);
  fHitPairTree->Branch("CellDiff",               &fCellDiff);
  fHitPairTree->Branch("ADCdiff",                &fAdcDiff);
  fHitPairTree->Branch("Xdiff",                  &fXdiff);
}

DEFINE_ART_MODULE(sn::SNSlicerAna)