FillNusSandbox_module.cc

Go to the documentation of this file.

/////////////////////////////////////////////////////////
/// \brief  Sandbox to make objects for nu_s analysis
///         in CAF
/// \author gsdavies@iu.edu - Gavin Davies
////////////////////////////////////////////////////////
#include <cmath>
#include <string>
#include <vector>

// ROOT includes
#include "TMath.h"
#include "TVector3.h"

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "NovaDAQConventions/DAQConventions.h"
#include "NusSandbox/NusSandObj.h"
#include "Geometry/Geometry.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Track.h"
#include "RecoBase/Vertex.h"
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/MathUtil.h"


namespace nussand {

  class FillNusSandbox : public art::EDProducer {
  public:
    explicit FillNusSandbox(fhicl::ParameterSet const & pset);
    virtual ~FillNusSandbox();

    void produce(art::Event & evt);

  private:

    std::string fSlicerLabel;
    std::string fVertexLabel;
    std::string fTrackLabel;

  };
}


namespace nussand
{

  //----------------------------------------------------------------------
  FillNusSandbox::FillNusSandbox(fhicl::ParameterSet const & pset):
    fSlicerLabel (pset.get< std::string >("SlicerLabel")),
    fVertexLabel (pset.get< std::string >("VertexLabel")),
    fTrackLabel  (pset.get< std::string >("TrackLabel"))
  {
    produces< std::vector<nussand::NusSandObj> >();
    produces< art::Assns<nussand::NusSandObj, rb::Cluster> >();
  }


  //----------------------------------------------------------------------
  FillNusSandbox::~FillNusSandbox()
  {
  }

  //-------------------------------------------------------------------
  void FillNusSandbox::produce(art::Event& evt)
  {

    art::Handle<std::vector<rb::Cluster> > slicevec;
    evt.getByLabel(fSlicerLabel,slicevec);

    if(slicevec->empty()) {
      LOG_WARNING("NusSandbox") << "No Slices in the input file";
      return;
    }

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

    art::FindManyP<rb::Track>  fmpTrack(slicevec, evt, fTrackLabel);
    art::FindManyP<rb::Vertex> fmpVert(slicevec, evt, fVertexLabel);

    std::unique_ptr< std::vector<nussand::NusSandObj> > outputObjects(new std::vector<nussand::NusSandObj>);
    std::unique_ptr< art::Assns<nussand::NusSandObj, rb::Cluster> > assoc(new art::Assns<nussand::NusSandObj, rb::Cluster>);

    const TVector3 beamdir = geom->NuMIBeamDirection();

    const double beamx = beamdir.X();
    const double beamy = beamdir.Y();
    const double beamz = beamdir.Z();

    const double xbeam  = beamx/util::pythag(beamx, beamz);
    const double ybeam  = beamy/util::pythag(beamy, beamz);
    const double zxbeam = beamz/util::pythag(beamx, beamz);
    const double zybeam = beamz/util::pythag(beamy, beamz);

    for(size_t slcIdx = 0; slcIdx < slicevec->size(); ++slcIdx){

      art::Ptr<rb::Cluster> slice(slicevec,slcIdx);

      if(slice->IsNoise()) continue;

      nussand::NusSandObj sandbox;

      //Make sure there is a vertex associated with this slice
      if(!fmpVert.isValid()) { 
        LOG_WARNING("NusSandbox") << "No " << fVertexLabel.c_str() << "Vertices!";
        outputObjects->push_back(sandbox);
        util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
        continue;
      }
      
      std::vector< art::Ptr<rb::Vertex> > v = fmpVert.at(slcIdx);
      if(v.size() != 1) continue;// Continue if there is no vertex or more then 1 vertex

      if(!fmpTrack.isValid()) {
        LOG_WARNING("NusSandbox") << "No " << fTrackLabel.c_str() << " Tracks!";
        outputObjects->push_back(sandbox);
        util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
        continue;
      }

      double xFinal   = 0.;
      double yFinal   = 0.;
      double zxFinal  = 0.;
      double zyFinal  = 0.;
      double exFinal  = 0.;
      double eyFinal  = 0.;
      double ezxFinal = 0.;
      double ezyFinal = 0.;
      
      std::vector< art::Ptr<rb::Track> > tracks = fmpTrack.at(slcIdx);
      for(size_t trkIdx = 0; trkIdx < tracks.size(); ++trkIdx){
        
        const art::Ptr<rb::Track> track = tracks[trkIdx];

        TVector3 tDir = track->Dir().Unit();

        double xtrkE = 0.;
        double ytrkE = 0.;
        // Loop over XCells
        if(track->View() == geo::kX){
          for(size_t hitIdx = 0; hitIdx < track->NXCell(); ++hitIdx){
            const art::Ptr<rb::CellHit> chit = track->XCell(hitIdx);
            const rb::RecoHit rhit = track->RecoHit(chit);
            if(rhit.IsCalibrated()) xtrkE += rhit.GeV()*track->Weight(geo::kX, hitIdx);
          } // loop over xhits
          
          TVector3 tdirx  = track->Dir().Unit();
          xFinal   += fabs(tDir.X());
          zxFinal  += tDir.Z();
          exFinal  += fabs(tDir.X())*xtrkE;
          ezxFinal += tDir.Z()*xtrkE;
        }
        else if(track->View() == geo::kY){
          // Loop over YCells
          for(size_t hitIdx = 0; hitIdx < track->NYCell(); ++hitIdx){
            const art::Ptr<rb::CellHit> chit = track->YCell(hitIdx);
            const rb::RecoHit rhit = track->RecoHit(chit);
            if(rhit.IsCalibrated()) ytrkE += rhit.GeV()*track->Weight(geo::kY, hitIdx);
          } // loop over yhits 
          yFinal   += fabs(tDir.Y());
          zyFinal  += tDir.Z();
          eyFinal  += fabs(tDir.Y())*ytrkE;
          ezyFinal += tDir.Z()*ytrkE;
        }
        else{
          LOG_INFO("NusSandbox") << "Not a geo::kX or geo::kY view. Don't do anything!";
          continue;
        }

      } // loop over tracks

      double fTXZ  = util::pythag(xFinal, zxFinal);
      double fTYZ  = util::pythag(yFinal, zyFinal);
      double fTEXZ = util::pythag(exFinal, ezxFinal);
      double fTEYZ = util::pythag(eyFinal, ezyFinal);

      sandbox.SetSumTX(fTXZ);
      sandbox.SetSumTY(fTYZ);
      sandbox.SetEWSumTX(fTEXZ);
      sandbox.SetEWSumTY(fTEYZ);
      
      double beamTrkXZ  = ((xFinal*xbeam + zxFinal*zxbeam)/fTXZ);
      double beamTrkYZ  = ((yFinal*ybeam + zyFinal*zybeam)/fTYZ);
      double beamTrkEXZ = ((exFinal*xbeam + ezxFinal*zxbeam)/fTEXZ);
      double beamTrkEYZ = ((eyFinal*ybeam + ezyFinal*zybeam)/fTEYZ);

      sandbox.SetCosSumTX(beamTrkXZ);
      sandbox.SetCosEWSumTX(beamTrkEXZ);
      sandbox.SetCosSumTY(beamTrkYZ);
      sandbox.SetCosEWSumTY(beamTrkEYZ);

      sandbox.SetAngSumTX(acos(beamTrkXZ)*180.0/TMath::Pi());
      sandbox.SetAngEWSumTX(acos(beamTrkEXZ)*180.0/TMath::Pi());
      sandbox.SetAngSumTY(acos(beamTrkYZ)*180.0/TMath::Pi());
      sandbox.SetAngEWSumTY(acos(beamTrkEYZ)*180.0/TMath::Pi());

      outputObjects->push_back(sandbox);
      util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
    } //end loop over slices

    evt.put(std::move(outputObjects));
    evt.put(std::move(assoc));

  } //end produce

} //end namespace


DEFINE_ART_MODULE(nussand::FillNusSandbox)