ClusterMergeAlg.h

Go to the documentation of this file.

// \file ClusterMergeAlg.h
//
// \brief ClusterMergeAlg header file
//
// \author david.kaleko@gmail.com, kazuhiro@nevis.columbia.edu
//

#ifndef CLUSTERMERGEALG_H
#define CLUSTERMERGEALG_H

// ART includes
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Registry/ServiceMacros.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"

// LArSoft
#include "nusimdata/SimulationBase/MCTruth.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "larreco/RecoAlg/SpacePointAlg.h"
#include "larcore/Geometry/Geometry.h"
// STL
#include <set>
#include <vector>
#include <sstream>

// ROOT
#include <TString.h>
#include <TTree.h>

namespace cluster
{

  struct cluster_merge_info {
    
    unsigned int cluster_index; 
    geo::View_t  view;          
    geo::PlaneID planeID;       
    
    float  start_wire;          
    float  start_time;          
    float  end_wire;            
    float  end_time;            
    
    double start_wire_err;      
    double start_time_err;      
    double end_wire_err;        
    double end_time_err;        
    
    float  angle;               
    
    cluster_merge_info(): planeID() {
      
      cluster_index = 0xffffffff;
      view = geo::kUnknown;
      start_wire = start_time = end_wire = end_time = -1;
      start_wire_err = start_time_err = end_wire_err =end_time_err = -1;
      
    };
    
    explicit cluster_merge_info(const recob::Cluster& cl)
      :cluster_index(cl.ID())
      ,view(cl.View())
      ,planeID(cl.Plane())
      ,start_wire(cl.StartWire())
      ,start_time(cl.StartTick())
      ,end_wire(cl.EndWire())
      ,end_time(cl.EndTick())
      ,angle(cl.StartAngle())
      {}
  };

  class ClusterMergeAlg {
    
  public:

    ClusterMergeAlg(fhicl::ParameterSet const& pset);
    //ClusterMergeAlg(fhicl::ParameterSet const& pset, art::ActivityRegistry& reg);

    virtual ~ClusterMergeAlg(){};

    void VerboseMode(bool on) { _verbose = on; }

    void ReportConfig() const;

    void SetAngleCut(double angle) { _max_allowed_2D_angle_diff = angle; }

    void SetSquaredDistanceCut(double d) { _max_2D_dist2 = d; }

    void AppendClusterInfo(const recob::Cluster &in_cluster,
                           const std::vector<art::Ptr<recob::Hit> > &in_hit_v);

    void AppendClusterInfo(const art::Ptr<recob::Cluster> in_cluster,
                           const std::vector<art::Ptr<recob::Hit> > &in_hit_v);

    void ClearEventInfo();

    void ProcessMergeAlg();

    const std::vector<std::vector<unsigned int> > GetClusterSets () const {return _cluster_sets_v;};

    bool CompareClusters(const cluster_merge_info &clus_info_A,
                         const cluster_merge_info &clus_info_B);

    bool Angle2DCompatibility(const cluster_merge_info &clus_info_A,
                              const cluster_merge_info &clus_info_B) const;

    bool ShortestDistanceCompatibility(const cluster_merge_info &clus_info_A,
                                       const cluster_merge_info &clus_info_B) const;

    double ShortestDistanceSquared(double point_x, double point_y, 
                                   double start_x, double start_y,
                                   double end_x,   double end_y  ) const;

    void PrintClusterVars(cluster_merge_info clus_info) const;

    int isInClusterSets(unsigned int index) const;

  protected:

    void SetWire2Cm(double f) { _wire_2_cm = f; }

    void SetTime2Cm(double f) { _time_2_cm = f; }

    void ClearOutputInfo();

    void ClearInputInfo();

    void ClearTTreeInfo();

    void PrepareTTree();

    void PrepareDetParams();

    void AppendHitInfo(cluster_merge_info &ci,
                       const std::vector<art::Ptr<recob::Hit> > &in_hit_v);

    void BuildClusterSets(const cluster_merge_info &clus_info_A,
                          const cluster_merge_info &clus_info_B);

    void FinalizeClusterSets();

    int AppendToClusterSets(unsigned int cluster_index, int merged_index=-1);

  protected:

    bool _verbose;             
    bool _det_params_prepared; 
    TTree* _merge_tree;        

    std::vector<int> _cluster_merged_index;

    std::vector<std::vector<unsigned int> > _cluster_sets_v;

    std::vector<cluster::cluster_merge_info> _u_clusters; 
    std::vector<cluster::cluster_merge_info> _v_clusters; 
    std::vector<cluster::cluster_merge_info> _w_clusters; 

    double _wire_2_cm; 
    double _time_2_cm; 
    double _max_allowed_2D_angle_diff; //in degrees
    double _max_2D_dist2;              //in cm^2
    double _min_distance_unit;         //in cm^2
  }; // class ClusterMergeAlg
  
} //namespace cluster
#endif