CosmicRayVertex_module.cc

Go to the documentation of this file.

///
/// \file    CosmicRayVertex_module.cc
/// \brief   Place a vertex on a cosmic ray track.
/// \author  messier@indiana.edu
/// \version $Id:$
///
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/FindMany.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"

#include "RecoBase/FilterList.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Vertex.h"
#include "RecoBase/Track.h"
#include "Utilities/AssociationUtil.h"
#include "GeometryObjects/GeometryBase.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/Geo.h"

/// Place a vertex on the detector wall for a cosmic-ray track
namespace crvtx {
  ///
  /// Module for placing a vertex on a cosmic ray track
  ///
  class CosmicRayVertex : public art::EDProducer {
  private:
    // folder labels
    std::string fSliceLabel;       ///< Where to get slices?
    std::string fTrackLabel;       ///< Where to look for tracks?
    bool        fObeyPreselection; ///< Check rb::IsFiltered?

  public:
    //..................................................................
    explicit CosmicRayVertex(fhicl::ParameterSet const &p) {
      this->reconfigure(p);
      this->produces< std::vector<rb::Vertex> > ();
      this->produces< art::Assns<rb::Vertex, rb::Cluster> >();
    }

    //..................................................................
    void reconfigure(const fhicl::ParameterSet& p) {
#     define SET(N) { f##N = p.get<__typeof__(f##N)>(#N); }
      SET(ObeyPreselection);
      SET(SliceLabel);
      SET(TrackLabel);
#     undef SET
    }
    
    //..................................................................
    virtual ~CosmicRayVertex() {  }
   
    //..................................................................
    void beginJob() { }
 
    //..................................................................
    void produce(art::Event& evt) {
      unsigned int i,j;
      //
      // We produce a vertex with an association to the slice clutser
      // from which it came
      //
      std::unique_ptr< std::vector<rb::Vertex> >
        vtx(new std::vector<rb::Vertex>);
      std::unique_ptr< art::Assns<rb::Vertex, rb::Cluster> > 
        vtx_to_slice(new art::Assns<rb::Vertex, rb::Cluster>);
      
      // 
      // Find the slice we want to analyze and any associated tracks
      //
      art::Handle<std::vector<rb::Cluster> > slice;
      evt.getByLabel(fSliceLabel,slice);
      
      art::FindMany<rb::Track> trackh(slice, evt, fTrackLabel);
      
      art::ServiceHandle<geo::Geometry> geom;
      
      //Vertex location
      double enterp[3];

      for(i = 0; i < slice->size(); ++i) {
        if(fObeyPreselection && rb::IsFiltered(evt, slice, i)) continue;
        //
        // Get a pointer to "this" slice. Skip it if its the noise slice
        //
        const rb::Cluster& s = (*slice)[i];
        if (s.IsNoise()) continue;
        //
        // Find the tracks associated to this slice
        //
        std::vector<const rb::Track*> tracks = trackh.at(i);
        //
        // Determine the locations of the vertices for each track
        // 
        for (j=0; j<tracks.size(); ++j) {
          //(x,y,z) track start point
          TVector3 st = tracks[j]->Start();
          //track direction
          TVector3 direct = tracks[j]->Dir();
          //(x,y,z) track start point (double)
          double start[3] = {st.X(),st.Y(),st.Z()};
          //Backwards along track direction (towards edge of detector)
          double dir[3] = {-direct.X(),-direct.Y(),-direct.Z()};
          //
          //Check if point in detector
          //
          bool is_inside = true;
          is_inside &= (start[0]>-geom->DetHalfWidth())&&(start[0]<geom->DetHalfWidth());
          is_inside &= (start[1]>-geom->DetHalfHeight())&&(start[1]<geom->DetHalfHeight());
          is_inside &= (start[2]>0)&&(start[2]<geom->DetLength());
          //If yes, vertex at extrapolated entrance point of cosmic ray
          if (is_inside) 
            geom->DEdge(start,dir,enterp);
          //If no, check to see if the track intercepts the detector in the opposite direction
          else {
      double oppDir[3] = {direct.X(), direct.Y(), direct.Z()};
      bool intersects = geom->IntersectsDetector(start, oppDir);
      //If it doesn't intersect, vertex at track start                  
      if (!intersects) {
        std::copy(std::begin(start), std::end(start), std::begin(enterp));
      }
      //If it does intersect, vertex at extrapolated entrance point of cosmic ray
      //searched for from the outside
      else {
        TVector3 p1;
        p1.SetX(st.X() + (1000*direct.X()));
        p1.SetY(st.Y() + (1000*direct.Y()));
        p1.SetZ(st.Z() + (1000*direct.Z()));
        geo::ClampRayToDetectorVolume(&st, &p1, &(*geom));
        double newStart[3] = {st.X(), st.Y(), st.Z()};
        std::copy(std::begin(newStart), std::end(newStart), std::begin(enterp));
      }
    }
          //Initial time of track
          double t0 = tracks[j]->MinTNS();
          
          //
          // Add to collection of vertices to store
          //
          rb::Vertex v(enterp[0],enterp[1],enterp[2],t0);
          vtx->push_back(v);
          //
          // Create associations between slice and vertex
          //
          art::Ptr<rb::Cluster> sptr(slice, i);
          util::CreateAssn(*this, evt, *vtx, sptr, *vtx_to_slice);
        }
      }
     
      evt.put(std::move(vtx));
      evt.put(std::move(vtx_to_slice));
    }
    
    //......................................................................
    void endJob() { }
  };
}
////////////////////////////////////////////////////////////////////////
namespace crvtx
{
  DEFINE_ART_MODULE(CosmicRayVertex)
}
////////////////////////////////////////////////////////////////////////