ShowerMuonCoincidence_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       ShowerMuonCoincidence
// Plugin Type: analyzer (art v2_12_01)
// File:        ShowerMuonCoincidence_module.cc
//
// Generated at Mon Jul 15 11:08:30 2019 by Dung Phan using cetskelgen
// from cetlib version v3_06_00.
////////////////////////////////////////////////////////////////////////

#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>

namespace nnbar {
class ShowerMuonCoincidence;
}

class nnbar::ShowerMuonCoincidence : public art::EDAnalyzer {
 public:
  explicit ShowerMuonCoincidence(fhicl::ParameterSet const &p);
  ShowerMuonCoincidence(ShowerMuonCoincidence const &) = delete;
  ShowerMuonCoincidence(ShowerMuonCoincidence &&) = delete;
  ShowerMuonCoincidence &operator=(ShowerMuonCoincidence const &) = delete;
  ShowerMuonCoincidence &operator=(ShowerMuonCoincidence &&) = delete;

  // Required functions.
  void analyze(art::Event const &e) override;
  double distancePointToLine(TVector3 point, TVector3 pointOnLineA, TVector3 pointOnLineB);

  // Selected optional functions.
  void beginJob() override;
  void beginRun(art::Run const &r) override;
  void beginSubRun(art::SubRun const &sr) override;
  void endJob() override;
  void endRun(art::Run const &r) override;
  void endSubRun(art::SubRun const &sr) override;

 private:
  art::InputTag _ShowerClusterLabel;
  art::InputTag _MuonClusterLabel;
  art::InputTag _CosmicTrackLabel;

  TH1D* distFromShowerToCoincidenceMuon;
};

nnbar::ShowerMuonCoincidence::ShowerMuonCoincidence(
    fhicl::ParameterSet const &p)
    : EDAnalyzer(p) {
  _ShowerClusterLabel = "shower:showers";
  _MuonClusterLabel = "shower:cosmics";
  _CosmicTrackLabel = "windowtrack";
}

void nnbar::ShowerMuonCoincidence::beginJob() {
  art::ServiceHandle<art::TFileService> tfs;
  distFromShowerToCoincidenceMuon = tfs->make<TH1D>("distFromShowerToCoincidenceMuon", ";Distance (cm);", 300, 0, 6000);
}

void nnbar::ShowerMuonCoincidence::beginRun(art::Run const &r) {}

void nnbar::ShowerMuonCoincidence::beginSubRun(art::SubRun const &sr) {}

void nnbar::ShowerMuonCoincidence::analyze(art::Event const &e) {
  art::Handle<std::vector<rb::Cluster> > muonClusterHandle;
  e.getByLabel(_MuonClusterLabel, muonClusterHandle);

  art::Handle<std::vector<rb::Cluster> > showerClusterHandle;
  e.getByLabel(_ShowerClusterLabel, showerClusterHandle);

  art::FindManyP<rb::Track> fmtrack(muonClusterHandle, e, "windowtrack");

  for (unsigned int iShower = 0; iShower < showerClusterHandle->size(); iShower++) {
    rb::Cluster showerSlice = (*showerClusterHandle)[iShower];
    auto hitTimeMin = showerSlice.MinTNS();
    auto hitTimeMax = showerSlice.MaxTNS();

    auto coincidenceWindowMin = hitTimeMin - 500;
    auto coincidenceWindowMax = hitTimeMax + 500;
    unsigned int countSlicers = 0;
    std::vector<art::Ptr<rb::Track>> tracks;
    for (unsigned int iMuon = 0; iMuon < muonClusterHandle->size(); iMuon++){
      rb::Cluster muonSlice = (*muonClusterHandle)[iMuon];
      auto muonHitTimeMin = muonSlice.MinTNS();
      auto muonHitTimeMax = muonSlice.MaxTNS();
      bool muonBeforeShower = coincidenceWindowMin <= muonHitTimeMax && muonHitTimeMax <= coincidenceWindowMax;
      bool muonAfterShower = coincidenceWindowMin <= muonHitTimeMin && muonHitTimeMin <= coincidenceWindowMax;
      if (muonBeforeShower && muonAfterShower) {
        countSlicers++;
        if (fmtrack.isValid()){
          tracks = fmtrack.at(iMuon);
        }
      }
    }

    if (countSlicers == 1) {
      // find vertex
      double vertexXinCM;
      double vertexYinCM;
      double vertexZinCM;
      nnbar::NNbarUtilities::energyBalancingVertex(showerSlice, vertexXinCM, vertexYinCM, vertexZinCM);
      TVector3 point(vertexXinCM, vertexYinCM, vertexZinCM);

      // check associated track
      if (tracks.size() == 1) { 
        if (tracks.at(0)->Is3D()) {
          auto startPoint = tracks.at(0)->Start();
          auto endPoint = tracks.at(0)->Stop();
          double dist = distancePointToLine(point, startPoint, endPoint);
          std::cout << "dist: " << dist << std::endl;
          distFromShowerToCoincidenceMuon->Fill(dist);
        }
      }
    }
  }
}

double nnbar::ShowerMuonCoincidence::distancePointToLine(TVector3 point0, TVector3 point1, TVector3 point2) {
  // See: http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html

  double distSq21 = (point2 - point1).Mag2();
  double distSq10 = (point1 - point0).Mag2();
  double scalarProd = (point1 - point0).Dot(point2 - point1);
  double numerator = distSq10 * distSq21 - TMath::Power(scalarProd, 2);
  return TMath::Sqrt(numerator / distSq21);
}

void nnbar::ShowerMuonCoincidence::endJob() {
}

void nnbar::ShowerMuonCoincidence::endRun(art::Run const &r) {}

void nnbar::ShowerMuonCoincidence::endSubRun(art::SubRun const &sr) {}

DEFINE_ART_MODULE(nnbar::ShowerMuonCoincidence)