WCTrackReco_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////
/// \file    WCTrackReco_module.cc
/// \module  producer
/// \brief   Reconstruction for the wire chambers used in the beamline.
///          Part of the beamline reconstruction for NOvA test beam.
/// \author  Mike Wallbank (University of Cincinnati) <wallbank@fnal.gov>
///          Fan Gao (University of Pittsburgh) <fag16@pitt.edu>
/// \date    December 2018
////////////////////////////////////////////////////////////////////////////

// framework
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "canvas/Persistency/Common/FindManyP.h"

// nova
#include "BeamlineReco/WCTrackAlg.h"
#include "BeamlineReco/WCHitFinderAlg.h"
#include "BeamlineUtils/BeamlineGeometry.h"
#include "BeamlineUtils/BeamlineMagneticField.h"
#include "BeamlineRecoBase/WCTrack.h"
#include "RawData/RawBeamline.h"
#include "RecoBase/Vertex.h"
#include "SummaryData/MCenterData.h"
#include "Utilities/AssociationUtil.h"

// stl
#include <iostream>

// root

// -----------------------------------------------------------------------
namespace beamlinereco {

  class WCTrackReco : public art::EDProducer {

  public:

    explicit WCTrackReco(const fhicl::ParameterSet& pset);

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

  private:

    // Data labels
    art::InputTag fRawWCDataLabel;
    art::InputTag fRawBeamlineConfigLabel;
    art::InputTag fMCenterDataLabel;

    // Configuration
    bool fVerbose;
    bool fPickyTracks;
    bool fUseIFDB;
    double fVertexTime;

    // Histograms
    TH1F* fReco_P;
    TH1F* fY_Kink;
    TH1F* fX_Dist;
    TH1F* fY_Dist;
    TH1F* fZ_Dist;
    TH1F* fX_Mag_Int;
    TH1F* fY_Mag_Int;
    TH1F* fZ_Mag_Int;
    TH1F* fDist_to_Mag_Axis;
    TH1F* fX_Face_Dist;
    TH1F* fY_Face_Dist;
    TH1F* fTheta_Dist;
    TH1F* fPhi_Dist;
    TH1F* fTrack_Type;
    TH1F* fWCDist;

    // Algorithms and services
    WCHitFinderAlg fWCHitFinderAlg;
    WCTrackAlg fWCTrackAlg;
    art::ServiceHandle<beamlineutil::BeamlineMagneticField> fMagneticField;

  };

  DEFINE_ART_MODULE(WCTrackReco)

}

// -----------------------------------------------------------------------
beamlinereco::WCTrackReco::WCTrackReco(const fhicl::ParameterSet& pset) :
  fWCHitFinderAlg(pset.get<fhicl::ParameterSet>("WCHitFinderAlg")),
  fWCTrackAlg(pset.get<fhicl::ParameterSet>("WCTrackAlg")) {

  this->reconfigure(pset);
  produces<std::vector<brb::WCTrack> >();
  produces<art::Assns<brb::WCTrack, rawdata::RawBeamlineWC> >();
  produces<std::vector<rb::Vertex> > ();
  produces<art::Assns<rb::Vertex, brb::WCTrack> > ();
}

// -----------------------------------------------------------------------
void beamlinereco::WCTrackReco::reconfigure(const fhicl::ParameterSet& pset) {
  fRawWCDataLabel         = pset.get<art::InputTag>("RawWCDataLabel");
  fRawBeamlineConfigLabel = pset.get<std::string>  ("RawBeamlineConfigLabel");
  fMCenterDataLabel       = pset.get<std::string>  ("MCenterDataLabel");
  fVerbose                = pset.get<bool>         ("Verbose", false);
  fUseIFDB                = pset.get<bool>         ("UseIFDB");
  fVertexTime             = pset.get<double>       ("VertexTime");
}

// -----------------------------------------------------------------------
void beamlinereco::WCTrackReco::beginJob() {
  art::ServiceHandle<art::TFileService> tfs;

  fWCDist = tfs->make<TH1F>("WCCond","WC Conditions",7,0,7);
  fReco_P = tfs->make<TH1F>("Reco_P","Reconstructed momentum", 180, 0, 1800);
  fY_Kink = tfs->make<TH1F>("Y_Kink","Angle between US/DS tracks in Y direction (degrees)",200,-10*3.1415926/180,10*3.141592654/180);
  fX_Dist = tfs->make<TH1F>("X_Dist","X distance between US/DS tracks at midplane (mm)",1200,-60,1260); //Do not trust these values
  fY_Dist = tfs->make<TH1F>("Y_Dist","Y distance between US/DS tracks at midplane (mm)",1200,-600,600); //Do not trust these values
  fZ_Dist = tfs->make<TH1F>("Z_Dist","Z distance between US/DS tracks at midplane (mm)",1200,-60,1260); //Do not trust these values
  fX_Mag_Int = tfs->make<TH1F>("X_Mag_Int","X intercept of WC track with magnet front face (cm)",100,-150,-100);
  fX_Mag_Int->GetXaxis()->SetTitle("X intercept [cm]");
  fX_Mag_Int->GetYaxis()->SetTitle("Tracks per 2 cm");
  fY_Mag_Int = tfs->make<TH1F>("Y_Mag_Int","Y intercept of WC track with magnet front face (cm)",100,-10,10);
  fY_Mag_Int->GetXaxis()->SetTitle("Y intercept [cm]");
  fY_Mag_Int->GetYaxis()->SetTitle("Tracks per 0.22 cm");
  fZ_Mag_Int = tfs->make<TH1F>("Z_Mag_Int","Z intercept of WC track with magnet front face (cm)",100,415,425);
  fZ_Mag_Int->GetXaxis()->SetTitle("Z intercept [cm]");
  fZ_Mag_Int->GetYaxis()->SetTitle("Tracks per 0.1 cm");
  fDist_to_Mag_Axis = tfs->make<TH1F>("Dist_to_Mag_Axis","Transverse distance to central axis of magnet [cm]",80,-20,20);
  fDist_to_Mag_Axis->GetXaxis()->SetTitle("Transverse distance from WC track to central axis of magnet [cm]");
  fDist_to_Mag_Axis->GetXaxis()->SetTitle("Tracks per 0.5 cm");
  fX_Face_Dist = tfs->make<TH1F>("X_Face","X Location of Track's Detector Entry (mm)",1600,-200,1400);
  fY_Face_Dist = tfs->make<TH1F>("Y_Face","Y Location of Track's Detector Entry (mm)",800,-400,400);
  fTheta_Dist = tfs->make<TH1F>("Theta","Track Theta (w.r.t. Detector Z axis), (radians),",400,-.4,0.4);
  fPhi_Dist = tfs->make<TH1F>("Phi","Track Phi (w.r.t. Detector X axis), (radians)",2000,-6.28318,6.28318);
  fReco_P->GetXaxis()->SetTitle("Reconstructed momentum (MeV/c)");
  fReco_P->GetYaxis()->SetTitle("Tracks per 10 MeV/c");
  fY_Kink->GetXaxis()->SetTitle("Reconstructed y_kink (radians)");
  fY_Kink->GetYaxis()->SetTitle("Tracks per 0.000872 radians");
  fX_Dist->GetXaxis()->SetTitle("X distance between US and DS track ends");
  fX_Dist->GetYaxis()->SetTitle("Tracks per 1 mm - DO NOT TRUST THESE VALUES");
  fY_Dist->GetXaxis()->SetTitle("Y distance between US and DS track ends");
  fY_Dist->GetYaxis()->SetTitle("Tracks per 1 mm - DO NOT TRUST THESE VALUES");
  fZ_Dist->GetXaxis()->SetTitle("Z distance between US and DS track ends");
  fZ_Dist->GetYaxis()->SetTitle("Tracks per 1 mm - DO NOT TRUST THESE VALUES");
  fX_Face_Dist->GetXaxis()->SetTitle("X (mm)");
  fX_Face_Dist->GetYaxis()->SetTitle("Tracks per 1 mm");
  fY_Face_Dist->GetXaxis()->SetTitle("Y (mm)");
  fY_Face_Dist->GetYaxis()->SetTitle("Tracks per 1 mm");
  fTheta_Dist->GetXaxis()->SetTitle("Theta (radians)");
  fTheta_Dist->GetYaxis()->SetTitle("Tracks per .002 radians");
  fPhi_Dist->GetXaxis()->SetTitle("Phi (radians)");
  fPhi_Dist->GetYaxis()->SetTitle("Tracks per 0.00628 radians");
  fTrack_Type = tfs->make<TH1F>("TrackType","WCTrack conditions: 1=missHit,2=uniqueHits,3=lonelyHit,4=socialHits",4,0,4);
  fTrack_Type->GetYaxis()->SetTitle("# Events");
  fTrack_Type->GetXaxis()->SetTitle("Track Conditions");

}

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

  // data products to be made
  std::unique_ptr<std::vector<brb::WCTrack> > wc_reco(new std::vector<brb::WCTrack>);
  std::unique_ptr<art::Assns<brb::WCTrack, rawdata::RawBeamlineWC> >
    wc_assn(new art::Assns<brb::WCTrack, rawdata::RawBeamlineWC>);
  std::unique_ptr<std::vector<rb::Vertex> > wc_vtx(new std::vector<rb::Vertex>);
  std::unique_ptr<art::Assns<rb::Vertex, brb::WCTrack> >
    vtx_to_wc_trk(new art::Assns<rb::Vertex, brb::WCTrack>);

  // get the raw WC digits
  art::Handle<std::vector<rawdata::RawBeamlineWC> > wcHandle;
  std::vector<art::Ptr<rawdata::RawBeamlineWC> > wcs;
  if (evt.getByLabel(fRawWCDataLabel, wcHandle))
    art::fill_ptr_vector(wcs, wcHandle);

  art::Handle<sumdata::MCenterData> ifdb;
  evt.getByLabel(fMCenterDataLabel, ifdb);
  // get the magnetic field
  float b_field = 0.;
  // If spill is not found, magnet current will be set to -1e10.
  // In this case, fallback on final else to use beamline run to find magnet current in fcl.
  if ( fUseIFDB && ifdb.isValid() && ifdb->mc7magnet!=-1e10 )
    b_field = fMagneticField->FieldFromCurrent(ifdb->mc7magnet);
  else if ( fUseIFDB && !ifdb.isValid() ) {
      mf::LogInfo("TBSpillInfo") << "MW: Change in reconstruction, March 27.\n"
                                 << "The magnet information is now obtained from IFDB via the TBSpillInfo module, which need to be run first.\n"
                                 << "Be sure to have run this module before WCTrackReco.  (Alternatively, turn off the fUseIFDB flag (temporary))."
                                 << std::endl;
      abort();
  }
  else {
    art::Handle<rawdata::RawBeamlineConfig> blConfigHandle;
    unsigned int beamline_run = 0;
    //float b_field = -1.;
    if (evt.getByLabel(fRawBeamlineConfigLabel, blConfigHandle)) {
      //b_field = blConfigHandle->BField();
      beamline_run = blConfigHandle->BeamlineRun();
    }
    try {
      b_field = fMagneticField->MagneticField(beamline_run);
    }
    catch (cet::exception& e) {
      e.what();
    }
  }
  fWCTrackAlg.setMagneticField(b_field);

  // Find hits
  std::vector<std::vector<WCHitList> > good_hits; // Two vectors: WC#, axis
  fWCHitFinderAlg.createHits(wcs,
                             good_hits,
                             fVerbose);

  // Find tracks
  std::vector<double> reco_p_list;
  std::vector<TVector3> mag_int_list;
  std::vector<double> dist_to_mag_axis_list;
  std::vector<TVector2> face_list;
  std::vector<std::vector<TVector3> > wc_hit_pos_list;
  std::vector<TVector3> dir_list;
  std::vector<double> y_kink_list;
  std::vector<TVector3> dist_list;
  std::vector<WCHitList> final_tracks;
  float residual;
  int WCMissed;
  fWCTrackAlg.reconstructTracks(reco_p_list,
                                mag_int_list,
                                dist_to_mag_axis_list,
                                face_list,
                                wc_hit_pos_list,
                                dir_list,
                                final_tracks,
                                good_hits,
                                y_kink_list,
                                dist_list,
                                WCMissed,
                                fWCDist,
                                residual);

  for (unsigned int iNewTrack = 0; iNewTrack < final_tracks.size(); ++iNewTrack) {

    WCHitList final_track = final_tracks[iNewTrack];

    // std::cout << "Making track with p " << reco_p_list[iNewTrack] << std::endl
    //        << "y kink " << y_kink_list[iNewTrack] << std::endl
    //        << "dist (" << dist_list[iNewTrack].X() << ", " << dist_list[iNewTrack].Y() << ", " << dist_list[iNewTrack].Z() << ")" << std::endl
    //        << "face (" << face_list[iNewTrack].X() << ", " << face_list[iNewTrack].Y() << ")" << std::endl
    //        << "wc 0 hit " << wc_hit_pos_list[iNewTrack][0].X() << ", " << wc_hit_pos_list[iNewTrack][0].Y() << ", " << wc_hit_pos_list[iNewTrack][0].Z() << ")" << std::endl
    //        << "wc 1 hit " << wc_hit_pos_list[iNewTrack][1].X() << ", " << wc_hit_pos_list[iNewTrack][1].Y() << ", " << wc_hit_pos_list[iNewTrack][1].Z() << ")" << std::endl
    //        << "wc 2 hit " << wc_hit_pos_list[iNewTrack][2].X() << ", " << wc_hit_pos_list[iNewTrack][2].Y() << ", " << wc_hit_pos_list[iNewTrack][2].Z() << ")" << std::endl
    //        << "wc 3 hit " << wc_hit_pos_list[iNewTrack][3].X() << ", " << wc_hit_pos_list[iNewTrack][3].Y() << ", " << wc_hit_pos_list[iNewTrack][3].Z() << ")" << std::endl
    //        << "dir " << dir_list[iNewTrack].X() << ", " << dir_list[iNewTrack].Y() << ", " << dir_list[iNewTrack].Z() << ")" << std::endl
    //        << "theta " << dir_list[iNewTrack].Theta() << std::endl
    //        << "phi " << dir_list[iNewTrack].Phi() << std::endl
    //        << "residual " << residual << std::endl;

    brb::WCTrack the_track(reco_p_list[iNewTrack],
                           y_kink_list[iNewTrack],
                           dist_list[iNewTrack],
                           mag_int_list[iNewTrack],
                           dist_to_mag_axis_list[iNewTrack],
                           face_list[iNewTrack],
                           wc_hit_pos_list[iNewTrack],
                           dir_list[iNewTrack],
                           dir_list[iNewTrack].Theta(),
                           dir_list[iNewTrack].Phi(),
                           residual);

    fReco_P->Fill(reco_p_list[iNewTrack]);
    fY_Kink->Fill(y_kink_list[iNewTrack]);
    fX_Dist->Fill(dist_list[iNewTrack].X());
    fY_Dist->Fill(dist_list[iNewTrack].Y());
    fZ_Dist->Fill(dist_list[iNewTrack].Z());
    fX_Mag_Int->Fill(mag_int_list[iNewTrack].X());
    fY_Mag_Int->Fill(mag_int_list[iNewTrack].Y());
    fZ_Mag_Int->Fill(mag_int_list[iNewTrack].Z());
    fDist_to_Mag_Axis->Fill(dist_to_mag_axis_list[iNewTrack]);
    fX_Face_Dist->Fill(face_list[iNewTrack].X());
    fY_Face_Dist->Fill(face_list[iNewTrack].Y());
    fTheta_Dist->Fill(dir_list[iNewTrack].Theta());
    fPhi_Dist->Fill(dir_list[iNewTrack].Phi());
    
    // create vertex for detector
    rb::Vertex vtx(face_list[iNewTrack].X(),face_list[iNewTrack].Y(),0,fVertexTime);

    // add to output data product and make associations
    (*wc_reco).push_back(the_track);
    util::CreateAssn(*this, evt, *(wc_reco.get()), wcs, *(wc_assn.get()));
    (*wc_vtx).push_back(vtx);
    util::CreateAssn(*this, evt, *(vtx_to_wc_trk.get()), *(wc_vtx.get()), *(wc_reco.get()));

  }

  evt.put(std::move(wc_reco));
  evt.put(std::move(wc_assn));
  evt.put(std::move(wc_vtx));
  evt.put(std::move(vtx_to_wc_trk));

  return;

}