TBRecoAna_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////
// \file    TBRecoAna_module.cc
// \brief   Analyzer module for basic reco analysis of initial Beamline data files.
// \author  lackey32@fnal.gov
////////////////////////////////////////////////////////////////////////

// C/C++ includes
#include <string>

// ROOT includes
#include "TH1F.h"
#include "TH2F.h"
#include "TTimeStamp.h"
#include "TTree.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/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// NOvASoft includes
#include "BeamlineRecoBase/BeamlineDigit.h"
#include "BeamlineRecoBase/ToF.h"
#include "BeamlineRecoBase/WCTrack.h"
#include "Geometry/LiveGeometry.h"
#include "MCCheater/BackTracker.h"
#include "NovaTimingUtilities/TimingUtilities.h"
#include "RawData/RawTrigger.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/FitSum.h"
#include "RecoBase/Prong.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Shower.h"
#include "RecoBase/Vertex.h"


namespace tbem {


  class TBRecoAna : public art::EDAnalyzer {
  public:
    explicit TBRecoAna(fhicl::ParameterSet const& pset);
    ~TBRecoAna();
    void analyze(const art::Event& evt);
    void reconfigure(const fhicl::ParameterSet& pset);
    void beginJob();
    void endJob();

  private:

    std::string fRawDataLabel;
    std::string fClusterLabel;
    std::string fVertexLabel;
    std::string fProngModLabel;
    std::string fProng2DLabel;
    std::string fProng3DLabel;
    std::string fRecoToFLabel;
    std::string fRecoWCLabel;
    std::string fRecoCherenkovLabel;

    Int_t nhit, run, subrun, event;
    Double_t evtTime;
    Double_t toftime, wcmom, ckovadc;

    std::vector<int> hitPlane, hitCell;
    std::vector<double> hitTime;

    std::vector<int> slice;
    std::vector<int> sliceNhit, prongNhit;
    std::vector<double> cosVtxX, cosVtxY, cosVtxZ;
    std::vector<double> sliceCalE, sliceADCTot, sliceMaxPlane;
    std::vector<int> cosTrkNhit, winTrkNhit;
    std::vector<double> cosTrkStrtX, cosTrkStrtY, cosTrkStrtZ;
    std::vector<double> cosTrkStopX, cosTrkStopY, cosTrkStopZ;
    std::vector<double> cosTrkDirX, cosTrkDirY, cosTrkDirZ;
    std::vector<double> winTrkCalE, winTrkADCTot, winTrkMaxPlane;
    std::vector<double> winTrkStrtX, winTrkStrtY, winTrkStrtZ;
    std::vector<double> winTrkStopX, winTrkStopY, winTrkStopZ;
    std::vector<double> winTrkDirX, winTrkDirY, winTrkDirZ;
    std::vector<double> prongE, prongCalE, prongADCTot;
    std::vector<double> prongStrtX, prongStrtY, prongStrtZ;
    std::vector<double> prongDirX, prongDirY, prongDirZ;
    std::vector<double> prongHitTime;
    std::vector<double> chitHitTime;
    std::vector<double> cosTrkMinTime, winTrkMinTime, bpfTrkMinTime, prongMinTime;
    std::vector<double> cosTrkMeanTime, winTrkMeanTime, bpfTrkMeanTime, prongMeanTime;

    std::vector<int> bpfNhit;
    std::vector<double> bpfStrtX, bpfStrtY, bpfStrtZ;
    std::vector<double> bpfDirX, bpfDirY, bpfDirZ;
    std::vector<double> bpfTrkHitTime;

    TTree *emshower;

  };

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

  TBRecoAna::TBRecoAna(fhicl::ParameterSet const& pset)
    : EDAnalyzer(pset)
  {
    reconfigure(pset);
  }

  //.......................................................................
  
  TBRecoAna::~TBRecoAna()
  {
  }

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

  void TBRecoAna::reconfigure(const fhicl::ParameterSet& pset)
  {
    fRawDataLabel = pset.get<std::string>("RawDataLabel");
    fClusterLabel = pset.get<std::string>("ClusterLabel");
    fVertexLabel  = pset.get<std::string>("VertexLabel");
    fProngModLabel= pset.get<std::string>("ProngModLabel");
    fProng2DLabel = pset.get<std::string>("Prong2DLabel");
    fProng3DLabel = pset.get<std::string>("Prong3DLabel");
    fRecoToFLabel = pset.get<std::string>("RecoToFLabel");
    fRecoWCLabel  = pset.get<std::string>("RecoWCLabel");    
    fRecoCherenkovLabel = pset.get<std::string>("RecoCherenkovLabel");
  }

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

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

    emshower = tfs->make<TTree>("evts","");
    emshower->Branch("run",&run,"run/I");
    emshower->Branch("subrun",&subrun,"subrun/I");
    emshower->Branch("event",&event,"event/I");
    emshower->Branch("slice",&slice);
    emshower->Branch("evtTime",&evtTime,"evtTime/D");
    emshower->Branch("nhit",&nhit,"nhit/I");
    emshower->Branch("hitPlane",&hitPlane);
    emshower->Branch("hitCell",&hitCell);
    emshower->Branch("hitTime",&hitTime);
    emshower->Branch("sliceCalE",&sliceCalE);
    emshower->Branch("sliceADCTot",&sliceADCTot);
    emshower->Branch("sliceMaxPlane",&sliceMaxPlane);
    emshower->Branch("cosVtxX",&cosVtxX);
    emshower->Branch("cosVtxY",&cosVtxY);
    emshower->Branch("cosVtxZ",&cosVtxZ);
    emshower->Branch("prongE",&prongE);
    emshower->Branch("prongCalE",&prongCalE);
    emshower->Branch("prongADCTot",&prongADCTot);
    emshower->Branch("sliceNhit",&sliceNhit);

    emshower->Branch("cosTrkNhit",&cosTrkNhit);
    emshower->Branch("cosTrkMinTime",&cosTrkMinTime);
    emshower->Branch("cosTrkMeanTime",&cosTrkMeanTime);
    emshower->Branch("cosTrkStrtX",&cosTrkStrtX);
    emshower->Branch("cosTrkStrtY",&cosTrkStrtY);
    emshower->Branch("cosTrkStrtZ",&cosTrkStrtZ);
    emshower->Branch("cosTrkStopX",&cosTrkStopX);
    emshower->Branch("cosTrkStopY",&cosTrkStopY);
    emshower->Branch("cosTrkStopZ",&cosTrkStopZ);
    emshower->Branch("cosTrkDirX",&cosTrkDirX);
    emshower->Branch("cosTrkDirY",&cosTrkDirY);
    emshower->Branch("cosTrkDirZ",&cosTrkDirZ);
    
    emshower->Branch("winTrkNhit",&winTrkNhit);
    emshower->Branch("winTrkMinTime",&winTrkMinTime);
    emshower->Branch("winTrkMeanTime",&winTrkMeanTime);
    emshower->Branch("winTrkCalE",&winTrkCalE);
    emshower->Branch("winTrkADCTot",&winTrkADCTot);
    emshower->Branch("winTrkMaxPlane",&winTrkMaxPlane);
    emshower->Branch("winTrkStrtX",&winTrkStrtX);
    emshower->Branch("winTrkStrtY",&winTrkStrtY);
    emshower->Branch("winTrkStrtZ",&winTrkStrtZ);
    emshower->Branch("winTrkStopX",&winTrkStopX);
    emshower->Branch("winTrkStopY",&winTrkStopY);
    emshower->Branch("winTrkStopZ",&winTrkStopZ);
    emshower->Branch("winTrkDirX",&winTrkDirX);
    emshower->Branch("winTrkDirY",&winTrkDirY);
    emshower->Branch("winTrkDirZ",&winTrkDirZ);

    emshower->Branch("prongNhit",&prongNhit);
    emshower->Branch("prongMinTime",&prongMinTime);
    emshower->Branch("prongMeanTime",&prongMeanTime);
    emshower->Branch("prongStrtX",&prongStrtX);
    emshower->Branch("prongStrtY",&prongStrtY);
    emshower->Branch("prongStrtZ",&prongStrtZ);
    emshower->Branch("prongDirX",&prongDirX);
    emshower->Branch("prongDirY",&prongDirY);
    emshower->Branch("prongDirZ",&prongDirZ);
    emshower->Branch("prongHitTime",&prongHitTime);
    emshower->Branch("chitHitTime",&chitHitTime);

    emshower->Branch("bpfNhit",&bpfNhit);
    emshower->Branch("bpfTrkMinTime",&bpfTrkMinTime);
    emshower->Branch("bpfTrkMeanTime",&bpfTrkMeanTime);
    emshower->Branch("bpfStrtX",&bpfStrtX);
    emshower->Branch("bpfStrtY",&bpfStrtY);
    emshower->Branch("bpfStrtZ",&bpfStrtZ);
    emshower->Branch("bpfDirX",&bpfDirX);
    emshower->Branch("bpfDirY",&bpfDirY);
    emshower->Branch("bpfDirZ",&bpfDirZ);
    emshower->Branch("bpfTrkHitTime",&bpfTrkHitTime);

    emshower->Branch("toftime",&toftime,"toftime/D");
    emshower->Branch("wcmom",&wcmom,"wcmom/D");
    emshower->Branch("ckovadc",&ckovadc,"ckovadc/D");

  }

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

  void TBRecoAna::endJob()
  {

  }

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

  void TBRecoAna::analyze(const art::Event& evt)
  {
    nhit = -5;
    toftime = -5;
    wcmom   = -5;
    ckovadc = -5;

    // Get "header" info (not all in header, but general info about event)
    run = evt.run();
    subrun = evt.subRun();
    event = evt.event();

    // Get trigger time
    art::Handle<std::vector<rawdata::RawTrigger>> trigs;
    evt.getByLabel(fRawDataLabel, trigs);
    const rawdata::RawTrigger trig = trigs->at(0);

    struct timespec unixTime;
    novadaq::timeutils::convertNovaTimeToUnixTime(trig.fTriggerTimingMarker_TimeStart, unixTime);

    unsigned long int unixTimeSec = unixTime.tv_sec;
    unsigned long int unixTimeNanoSec = unixTime.tv_nsec;
    evtTime = unixTimeSec + (1e-9)*unixTimeNanoSec;



    art::ServiceHandle<geo::LiveGeometry> fLiveGeom;

    // Get Beamline info --------------------------------------------------
    // ToFs
    art::Handle<std::vector<brb::ToF>> tofHandle;
    std::vector<art::Ptr<brb::ToF>>    tofs;
    if (evt.getByLabel(fRecoToFLabel, tofHandle))
      art::fill_ptr_vector(tofs, tofHandle);

    for (auto tof : tofs) {
      double t = tof->Time();
      //t -= 21.0;
      toftime = t;
    }

    //WCs
    art::Handle<std::vector<brb::WCTrack>> wcHandle;
    std::vector<art::Ptr<brb::WCTrack>>    wcs;
    if (evt.getByLabel(fRecoWCLabel, wcHandle))
      art::fill_ptr_vector(wcs, wcHandle);

    for (auto wc : wcs) {
      wcmom = wc->Momentum();
    }


    //Ckov
    art::Handle<std::vector<brb::BeamlineDigit>> cherenkovHandle;
    std::vector<art::Ptr<brb::BeamlineDigit>>    cherenkovs;
    if (evt.getByLabel(fRecoCherenkovLabel, cherenkovHandle))
      art::fill_ptr_vector(cherenkovs, cherenkovHandle);
    
    for (auto ckov : cherenkovs) {
      ckovadc = ckov->PeakADC();
    }


    // Get Detector info --------------------------------------------------

    // Get list of CellHits for event
    art::Handle< std::vector<rb::CellHit> > cellhits;
    evt.getByLabel("calhit", cellhits);
    nhit = cellhits->size();
    // Loop over all CellHits to get location and timing information
    for (unsigned int hitIdx=0; hitIdx<cellhits->size(); ++hitIdx){
      art::Ptr<rb::CellHit> ihit(cellhits,hitIdx);
      hitPlane.push_back(ihit->Plane());
      hitCell.push_back(ihit->Cell());
      hitTime.push_back(ihit->TNS());
    }


    // Get Slices for event
    art::Handle<std::vector<rb::Cluster>> sliceHandle;
    evt.getByLabel(fClusterLabel, sliceHandle);

    // get associations between slices and other objects
    art::FindManyP<rb::Vertex> cosVtxAssn(sliceHandle, evt, "cosmicrayvertex");
    art::FindManyP<rb::Track> cosTrkAssn(sliceHandle, evt, "cosmictrack");
    art::FindManyP<rb::Track> winTrkAssn(sliceHandle, evt, "windowtrack");
    art::FindManyP<rb::Prong> fuzzykAssn(sliceHandle, evt, "fuzzykvertex");
    art::FindManyP<rb::Track> bpfTrkAssn(sliceHandle, evt, "breakpoint");

    art::PtrVector<rb::Cluster> slices;
    for ( unsigned int i=0; i<sliceHandle->size(); ++i )
      slices.push_back(art::Ptr<rb::Cluster>(sliceHandle, i));


    // Loop over slices
    for ( unsigned int sliceIdx = 0; sliceIdx<slices.size(); ++sliceIdx ){
      slice.push_back(sliceIdx);
      art::Ptr<rb::Cluster> islice = slices.at(sliceIdx);
      
      //sliceNhit.push_back(islice->NCell());
      // slice E
      if (islice->NCell() > 0) {
        //sliceCalE.push_back(islice->CalorimetricEnergy());
          
          // tracks
          TVector3 cosTrkStrt;
          TVector3 cosTrkStop;
          TVector3 winTrkStrt;
          TVector3 winTrkStop;

          if ( cosTrkAssn.isValid() ) {
            std::vector< art::Ptr< rb::Track > > tracks = cosTrkAssn.at(sliceIdx);
            for (unsigned int trackIdx = 0; trackIdx<tracks.size(); ++trackIdx ){
              art::Ptr<rb::Track> cosTrk = tracks.at(trackIdx);
              cosTrkStrt = cosTrk->Start();
              cosTrkStop = cosTrk->Stop();
              // force track to start downstream
              if ( cosTrkStop.Z() < cosTrkStrt.Z() ) {
                cosTrkStrt = cosTrk->Stop();
                cosTrkStop = cosTrk->Start();
              }
              cosTrkStrtX.push_back(cosTrkStrt.X());
              cosTrkStrtY.push_back(cosTrkStrt.Y());
              cosTrkStrtZ.push_back(cosTrkStrt.Z());
              cosTrkStopX.push_back(cosTrkStop.X());
              cosTrkStopY.push_back(cosTrkStop.Y());
              cosTrkStopZ.push_back(cosTrkStop.Z());
              cosTrkDirX.push_back(cosTrk->Dir().X());
              cosTrkDirY.push_back(cosTrk->Dir().Y());
              cosTrkDirZ.push_back(cosTrk->Dir().Z());
              cosTrkNhit.push_back(cosTrk->NCell());
              cosTrkMinTime.push_back(cosTrk->MinTNS());
              cosTrkMeanTime.push_back(cosTrk->MeanTNS());
            } // loop over cosmic tracks
          } // valid cosmic track


          if ( winTrkAssn.isValid() ) {
            std::vector< art::Ptr< rb::Track > > tracks = winTrkAssn.at(sliceIdx);
            for (unsigned int trackIdx = 0; trackIdx<tracks.size(); ++trackIdx ){
              art::Ptr<rb::Track> winTrk = tracks.at(trackIdx);
              winTrkStrt = winTrk->Start();
              winTrkStop = winTrk->Stop();
              // force track to start downstream
              if ( winTrkStop.Z() < winTrkStrt.Z() ) {
                winTrkStrt = winTrk->Stop();
                winTrkStop = winTrk->Start();
              }
              winTrkStrtX.push_back(winTrkStrt.X());
              winTrkStrtY.push_back(winTrkStrt.Y());
              winTrkStrtZ.push_back(winTrkStrt.Z());
              winTrkStopX.push_back(winTrkStop.X());
              winTrkStopY.push_back(winTrkStop.Y());
              winTrkStopZ.push_back(winTrkStop.Z());
              winTrkDirX.push_back(winTrk->Dir().X());
              winTrkDirY.push_back(winTrk->Dir().Y());
              winTrkDirZ.push_back(winTrk->Dir().Z());
              winTrkNhit.push_back(winTrk->NCell());
              winTrkMinTime.push_back(winTrk->MinTNS());
              winTrkMeanTime.push_back(winTrk->MeanTNS());

              winTrkCalE.push_back(winTrk->CalorimetricEnergy());
              winTrkADCTot.push_back(winTrk->TotalADC());
              winTrkMaxPlane.push_back(winTrk->MaxPlane());

              sliceNhit.push_back(islice->NCell());
              sliceCalE.push_back(islice->CalorimetricEnergy());
              sliceADCTot.push_back(islice->TotalADC());
              sliceMaxPlane.push_back(islice->MaxPlane());
            } // loop over window tracks
          } // valid window track


          
          // get prongs
          std::vector<art::Ptr<rb::Prong>> prongs3D;
          art::FindManyP<rb::Prong> fmProng3D(sliceHandle, evt,
                                              art::InputTag(fProngModLabel, fProng3DLabel));
          if ( !fmProng3D.isValid() ) continue;
          else prongs3D = fmProng3D.at(sliceIdx);
          if (prongs3D.size() == 0 ) continue;

          art::FindManyP<rb::Track> fmTrackBPF(prongs3D, evt, "breakpoint");

          // vertices
          std::vector<art::Ptr<rb::Vertex>> cosVtx;
          if (cosVtxAssn.isValid())
            cosVtx = cosVtxAssn.at(sliceIdx);
          if (cosVtx.size() == 1){
            cosVtxX.push_back(cosVtx.at(0)->GetX());
            cosVtxY.push_back(cosVtx.at(0)->GetY());
            cosVtxZ.push_back(cosVtx.at(0)->GetZ());
          }
          

          // prongs
          for (unsigned int prongIdx=0; prongIdx<prongs3D.size(); ++prongIdx){
            art::Ptr<rb::Prong> iprong = prongs3D.at(prongIdx);
            int nhitprong = iprong->NCell();
            prongNhit.push_back(nhitprong);
          
            prongStrtX.push_back(iprong->Start().X());
            prongStrtY.push_back(iprong->Start().Y());
            prongStrtZ.push_back(iprong->Start().Z());
            prongDirX.push_back(iprong->Dir().X());
            prongDirY.push_back(iprong->Dir().Y());
            prongDirZ.push_back(iprong->Dir().Z());     
            prongMinTime.push_back(iprong->MinTNS());
            prongMeanTime.push_back(iprong->MeanTNS());
            
            for (unsigned int hitIdx = 0; hitIdx<iprong->NCell(); ++hitIdx){
              const art::Ptr<rb::CellHit> chit = iprong->Cell(hitIdx);
              const rb::RecoHit rhit = iprong->RecoHit(chit);
              double hitTime = chit->TNS();
              // Events pre-July 3 11:42 CDT (pre run 100192) were 1ms, not 100μs.
              // Need to adjust time for timing cuts to work.
              if (run<100192)
                hitTime -= 450e3;

              chitHitTime.push_back(hitTime);
              if (hitIdx==0)
                prongHitTime.push_back(hitTime);

            } // end loop over hits
            
            prongE.push_back(iprong->TotalGeV());
            prongCalE.push_back(iprong->CalorimetricEnergy());
            prongADCTot.push_back(iprong->TotalADC());


            // Get bpf tracks
            std::vector<art::Ptr<rb::Track>> bpftracks;
            if(fmTrackBPF.isValid()) bpftracks = fmTrackBPF.at(prongIdx);

            art::FindManyP<rb::FitSum> fmFitsum(bpftracks, evt, "breakpoint");

            for (unsigned int trackIdx = 0; trackIdx < bpftracks.size(); ++trackIdx){
              
              std::vector<art::Ptr<rb::FitSum>> fitsums;
              if(fmFitsum.isValid()) fitsums = fmFitsum.at(trackIdx);
              if(fitsums.size()==0) continue;

              // Only look at muon assumption - this should be expanded
              if(fitsums[0]->PDG() == 13) {
                art::Ptr<rb::Track> bpfTrk = bpftracks[trackIdx];
                int nhitbpf = bpfTrk->NCell();
                bpfNhit.push_back(nhitbpf);
                bpfTrkMinTime.push_back(bpfTrk->MinTNS());
                bpfTrkMeanTime.push_back(bpfTrk->MeanTNS());
                
                bpfStrtX.push_back(bpfTrk->Start().X());
                bpfStrtY.push_back(bpfTrk->Start().Y());
                bpfStrtZ.push_back(bpfTrk->Start().Z());
                bpfDirX.push_back(bpfTrk->Dir().X());
                bpfDirY.push_back(bpfTrk->Dir().Y());
                bpfDirZ.push_back(bpfTrk->Dir().Z());       

            
                const art::Ptr<rb::CellHit> bpfhit = bpfTrk->Cell(0);
                double trkTime = bpfhit->TNS();
                // Events pre-July 3 11:42 CDT (pre run 100192) were 1ms, not 100μs.
                // Need to adjust time for timing cuts to work.
                if (run<100192)
                  trkTime -= 450e3;
                bpfTrkHitTime.push_back(trkTime);
              }
            } // loop over bpf tracks

          }//loop over prongs

      } // end non-noise slice
      
    } // end loop over slices


    emshower->Fill();

    // Cleanup vectors
    hitPlane.clear();
    hitCell.clear();
    hitTime.clear();
    
    slice.clear();
    cosVtxX.clear();
    cosVtxY.clear();
    cosVtxZ.clear();

    sliceNhit.clear();
    sliceCalE.clear();
    sliceADCTot.clear();
    sliceMaxPlane.clear();
    cosTrkNhit.clear();
    cosTrkMinTime.clear();
    cosTrkMeanTime.clear();
    cosTrkStrtX.clear();
    cosTrkStrtY.clear();
    cosTrkStrtZ.clear();
    cosTrkStopX.clear();
    cosTrkStopY.clear();
    cosTrkStopZ.clear();
    cosTrkDirX.clear();
    cosTrkDirY.clear();
    cosTrkDirZ.clear();
    winTrkNhit.clear();
    winTrkMinTime.clear();
    winTrkMeanTime.clear();
    winTrkCalE.clear();
    winTrkADCTot.clear();
    winTrkMaxPlane.clear();
    winTrkStrtX.clear();
    winTrkStrtY.clear();
    winTrkStrtZ.clear();
    winTrkStopX.clear();
    winTrkStopY.clear();
    winTrkStopZ.clear();
    winTrkDirX.clear();
    winTrkDirY.clear();
    winTrkDirZ.clear();

    prongNhit.clear();
    prongMinTime.clear();
    prongMeanTime.clear();
    prongE.clear();
    prongCalE.clear();
    prongADCTot.clear();
    prongStrtX.clear();
    prongStrtY.clear();
    prongStrtZ.clear();
    prongDirX.clear();
    prongDirY.clear();
    prongDirZ.clear();
    prongHitTime.clear();
    chitHitTime.clear();

    bpfNhit.clear();
    bpfTrkMinTime.clear();
    bpfTrkMeanTime.clear();
    bpfStrtX.clear();
    bpfStrtY.clear();
    bpfStrtZ.clear();
    bpfDirX.clear();
    bpfDirY.clear();
    bpfDirZ.clear();
    bpfTrkHitTime.clear();

  } // end analyze

  //.......................................................................
  DEFINE_ART_MODULE(tbem::TBRecoAna)
} // end namespace tbem
///////////////////////////////////////////////////////////////////