AlignLibToVtx_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    AlignLibToVtx_module.cc
// \brief   Rewrite a library file with events aligned according to vertex
// \author  Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/FileBlock.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 "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "TFile.h"
#include "TTree.h"

#include "Geometry/Geometry.h"
#include "GeometryObjects/OfflineChan.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Vertex.h"

#include "LEM/func/EventSummary.h"

#include "LEM/Util.h"

namespace lem
{
  /// Rewrite a library file with events aligned according to vertex
  class AlignLibToVtx: public art::EDProducer
  {
  public:
    explicit AlignLibToVtx(const fhicl::ParameterSet& pset);
    ~AlignLibToVtx();

    virtual void reconfigure(const fhicl::ParameterSet& pset);
    virtual void beginJob();
    virtual void respondToOpenInputFile(const art::FileBlock& fb);
    virtual void produce(art::Event& evt);

  protected:
    TTree* fTreeIn;
    TTree* fTreeOut;

    std::string fSlicerLabel;
    std::string fVertexLabel;
  };

  //.......................................................................
  AlignLibToVtx::AlignLibToVtx(const fhicl::ParameterSet& pset)
    : fTreeIn(0), fTreeOut(0)
  {
    reconfigure(pset);

    // Just to mark where the vertex is so you can check it's right
    produces<std::vector<geo::OfflineChan> >();
  }

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

  //......................................................................
  void AlignLibToVtx::reconfigure(const fhicl::ParameterSet& pset)
  {
    fVertexLabel = pset.get<std::string>("VertexLabel");
    fSlicerLabel = pset.get<std::string>("SlicerLabel");

    // We want to store exactly the values of the hits, only rescale them later
    // when actually doing matching.
    EventSummary::fgDecayRate = 0;
    EventSummary::fgChargePower = 1;
  }

  //......................................................................
  void AlignLibToVtx::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;
    fTreeOut = tfs->make<TTree>("tree", "tree");
    EventSummary::InitToTree(fTreeOut);
  }

  //......................................................................
  void AlignLibToVtx::respondToOpenInputFile(const art::FileBlock& fb)
  {
    TFile* fin = new TFile(fb.fileName().c_str());
    assert(!fin->IsZombie());
    fTreeIn = (TTree*)fin->Get("makelib/tree");
    assert(fTreeIn);

    EventSummary::InitFromTree(fTreeIn);
  }

  //......................................................................
  void AlignLibToVtx::produce(art::Event& evt)
  {
    std::unique_ptr<std::vector<geo::OfflineChan> > chancol(new std::vector<geo::OfflineChan>);

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

    art::Handle<std::vector<rb::Vertex> > vtxs;
    evt.getByLabel(fVertexLabel, vtxs);

    art::Handle<std::vector<rb::Cluster> > slices;
    evt.getByLabel(fSlicerLabel, slices);
    assert(slices->size() == 2); // The real slice, plus the noise slice
    const rb::Cluster& slice = (*slices)[1];
    assert(!slice.IsNoise());

    int plane, cell, cellOther;
    if(vtxs->size() == 1){
      VertexToPlaneAndCell((*vtxs)[0].GetXYZ(), slice, plane, cell, cellOther);
    }
    else{
      DefaultVertex(slice, plane, cell, cellOther);
    }

    chancol->push_back(geo::OfflineChan(plane, cell));
    // Should really be nextPlane, but that's lost inside of Util.cxx. This is
    // almost always right, and only for display purposes anyway.
    chancol->push_back(geo::OfflineChan(plane+1, cellOther));

    evt.put(std::move(chancol));


    const int trIdx = evt.event()-1;
    // Should be exactly the same file
    assert(trIdx < fTreeIn->GetEntries());

    // This is the offset GenFromLib applied to put the event in the middle of
    // the detector. Take it out again here.
    const int plane0 = geom->NPlanes()/2-128;
    const int cell0 = geom->Plane(0)->Ncells()/2-128;

    // Doesn't matter if we say it's enriched or not
    EventSummary sum = EventSummary::FromTree(fTreeIn, trIdx, false);

    // Occasionally fails, should just be a warning.
    //    assert(sum.hits.size() == slice.NCell());

    const int recoVtxPlane = plane-plane0;
    const int recoVtxCell = cell-cell0;
    const int recoVtxCellOtherView = cellOther-cell0;

    for(int i = 0; i < sum.nhits; ++i){
      lem::LiteHit& h = sum.hits[i];

      const int planeDelta = h.Plane()-recoVtxPlane;
      const int cellDelta = ((abs(planeDelta)%2) ?
                             h.Cell()-recoVtxCellOtherView :
                             h.Cell()-recoVtxCell);

      h.SetPlane(planeDelta+lem::kVertexPlane);
      h.SetCell(cellDelta+lem::kVertexCell);
    }

    // Transform the true vertex the same way
    const int planeDelta = sum.trueVtxPlane-recoVtxPlane;
    const int cellDelta = ((abs(planeDelta)%2) ?
                           sum.trueVtxCell-recoVtxCellOtherView :
                           sum.trueVtxCell-recoVtxCell);
    const int cellDeltaOther = ((abs(planeDelta)%2) ?
                                sum.trueVtxCellOther-recoVtxCell :
                                sum.trueVtxCellOther-recoVtxCellOtherView);

    sum.trueVtxPlane = planeDelta+lem::kVertexPlane;
    sum.trueVtxCell = cellDelta+lem::kVertexCell;
    sum.trueVtxCellOther = cellDeltaOther+lem::kVertexCell;

    sum.ToTree(fTreeOut);
  }

  DEFINE_ART_MODULE(AlignLibToVtx)

} // end namespace lem
////////////////////////////////////////////////////////////////////////