TTHitFinder_module.cc

Go to the documentation of this file.

#ifndef TTHITFINDER_H
#define TTHITFINDER_H

#include <string>
#include <math.h>

// Framework includes
#include "art/Framework/Core/ModuleMacros.h" 
#include "art/Framework/Principal/Event.h" 
#include "art/Framework/Core/EDProducer.h" 
#include "canvas/Persistency/Common/FindOneP.h"

// LArSoft Includes
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "larcore/Geometry/Geometry.h"
#include "lardataobj/RawData/RawDigit.h"
#include "lardataobj/RecoBase/Wire.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardata/ArtDataHelper/HitCreator.h"

namespace hit{

  class TTHitFinder : public art::EDProducer {
    
  public:
    
    explicit TTHitFinder(fhicl::ParameterSet const& pset); 
    virtual ~TTHitFinder();
         
    void produce(art::Event& evt); 
    void beginJob(); 
    void endJob(); 
    void reconfigure(fhicl::ParameterSet const& p);                

  private:
        
    std::string    fCalDataModuleLabel; 
    float          fMinSigPeakInd;      
    float          fMinSigPeakCol;      
    float          fMinSigTailInd;      
    float          fMinSigTailCol;      
    int            fIndWidth;           
    int            fColWidth;           

    float getTotalCharge(const float*,int,float);
   
  protected:     
  
  }; // class TTHitFinder  
  
  //-------------------------------------------------
  TTHitFinder::TTHitFinder(fhicl::ParameterSet const& pset) {
    this->reconfigure(pset);
    
    // let HitCollectionCreator declare that we are going to produce
    // hits and associations with wires and raw digits
    recob::HitCollectionCreator::declare_products(*this, "uhits");
    recob::HitCollectionCreator::declare_products(*this, "vhits");
    recob::HitCollectionCreator::declare_products(*this, "yhits");
    
  }

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

  //-------------------------------------------------
  void TTHitFinder::reconfigure(fhicl::ParameterSet const& p) {
    fCalDataModuleLabel = p.get< std::string >("CalDataModuleLabel");
    fMinSigPeakInd      = p.get< float       >("MinSigPeakInd");
    fMinSigPeakCol      = p.get< float       >("MinSigPeakCol"); 
    fMinSigTailInd      = p.get< float       >("MinSigTailInd",-99); //defaulting to well-below zero
    fMinSigTailCol      = p.get< float       >("MinSigTailCol",-99); //defaulting to well-below zero
    fIndWidth           = p.get< int         >("IndWidth", 3); //defaulting to 3  
    fColWidth           = p.get< int         >("ColWidth", 3); //defaulting to 3

    //enforce a minimum width
    if(fIndWidth<1){
      mf::LogWarning("TTHitFinder") << "IndWidth must be 1 at minimum. Resetting width to one time tick";
      fIndWidth = 1;
    }
    if(fColWidth<1){
      mf::LogWarning("TTHitFinder") << "ColWidth must be 1 at minimum. Resetting width to one time tick";
      fColWidth = 1;
    }

  }

  //-------------------------------------------------
  void TTHitFinder::beginJob(){}

  //-------------------------------------------------
  void TTHitFinder::endJob(){}

  //-------------------------------------------------
  void TTHitFinder::produce(art::Event& evt)
  { 
    
    // these objects contain the hit collections
    // and their associations to wires and raw digits:
    recob::HitCollectionCreator hitCollection_U(*this, evt, "uhits");
    recob::HitCollectionCreator hitCollection_V(*this, evt, "vhits");
    recob::HitCollectionCreator hitCollection_Y(*this, evt, "yhits");
  
    // Read in the wire List object(s).
    art::Handle< std::vector<recob::Wire> > wireVecHandle;
    evt.getByLabel(fCalDataModuleLabel,wireVecHandle);
    std::vector<recob::Wire> const& wireVec(*wireVecHandle);

    // also get the raw digits associated with wires
    art::FindOneP<raw::RawDigit> WireToRawDigits
      (wireVecHandle, evt, fCalDataModuleLabel);
    
    art::ServiceHandle<geo::Geometry> geom;
   
    //initialize some variables that will be in the loop.
    float threshold_peak = 0;
    float threshold_tail = -99;
    int   width = 3;

    //Loop over wires
    for(unsigned int wireIter = 0; wireIter < wireVec.size(); wireIter++) {
      
      //get our wire
      art::Ptr<recob::Wire> wire(wireVecHandle, wireIter);
      art::Ptr<raw::RawDigit> const& rawdigits = WireToRawDigits.at(wireIter);
      
      std::vector<float> signal(wire->Signal());
      std::vector<float>::iterator timeIter;   // iterator for time bins
      geo::WireID wire_id = (geom->ChannelToWire(wire->Channel())).at(0); //just grabbing the first one
      
      
      //set the thresholds and widths based on wire type
      geo::SigType_t sigType = geom->SignalType(wire->Channel());
      if(sigType == geo::kInduction){
        threshold_peak = fMinSigPeakInd;
        threshold_tail = fMinSigTailInd;
        width          = fIndWidth;
      }
      else if(sigType == geo::kCollection){
        threshold_peak = fMinSigPeakCol;
        threshold_tail = fMinSigTailCol;
        width          = fColWidth;
      }

      //make a half_width variable to be the search window around each time tick.
      float half_width = ((float)width-1)/2.;

      //now do the loop over the time ticks on the wire
      int time_bin = -1;
      float peak_val = 0;
      for(timeIter = signal.begin(); timeIter < signal.end(); timeIter++){
        time_bin++;

        //set the peak value, taking average between ticks if desired total width is even
        if(width%2==1) peak_val = *timeIter;
        else if(width%2==0) peak_val = 0.5 * (*timeIter + *(timeIter+1));

        //continue immediately if we are not above the threshold
        if(peak_val < threshold_peak) continue;

        //continue if we are too close to the edge
        if( time_bin-half_width < 0 ) continue;
        if( time_bin+half_width > signal.size() ) continue;
        
        //if necessary, do loop over hit width, and check tail thresholds
        int begin_tail_tick = std::floor(time_bin-half_width);
        float totalCharge = getTotalCharge(&signal[begin_tail_tick],width,threshold_tail);
        if(totalCharge==-999) {
          LOG_DEBUG("TTHitFinder") << "Rejecting would be hit at (plane,wire,time_bin,first_bin,last_bin)=(" 
                                   << wire_id.Plane << "," << wire_id.Wire << "," << time_bin << "," << begin_tail_tick << "," << begin_tail_tick+width-1 << "): " 
                                   << signal.at(time_bin-1) << " "
                                   << signal.at(time_bin) << " "
                                   << signal.at(time_bin+1);
          continue;
        }
        
        //OK, if we've passed all tests up to this point, we have a hit!

        float hit_time = time_bin;
        if(width%2==0) hit_time = time_bin+0.5;
        
        // hit time region is 2 widths (4 RMS) wide
        const raw::TDCtick_t start_tick = hit_time - width,
          end_tick = hit_time + width;
        
        // make the hit
        recob::HitCreator hit(
          *wire,                     // wire
          wire_id,                   // wireID
          start_tick,                // start_tick
          end_tick,                  // end_tick
          width / 2.,                // rms
          hit_time,                  // peak_time
          0.,                        // sigma_peak_time
          peak_val,                  // peak_amplitude
          0.,                        // sigma_peak_amplitude
          totalCharge,               // hit_integral
          0.,                        // hit_sigma_integral
          totalCharge,               // summedADC
          1,                         // multiplicity (dummy value)
          0,                         // local_index (dummy value)
          1.,                        // goodness_of_fit (dummy value)
          0                          // dof
          );
        if(wire_id.Plane==0)
          hitCollection_U.emplace_back(hit.move(), wire, rawdigits);
        else if(wire_id.Plane==1)
          hitCollection_V.emplace_back(hit.move(), wire, rawdigits);
        else if(wire_id.Plane==2)
          hitCollection_Y.emplace_back(hit.move(), wire, rawdigits);

      }//End loop over time ticks on wire

      LOG_DEBUG("TTHitFinder") << "Finished wire " << wire_id.Wire << " (plane " << wire_id.Plane << ")"
                               << "\tTotal hits (U,V,Y)= (" 
                               << hitCollection_U.size() << ","
                               << hitCollection_V.size() << ","
                               << hitCollection_Y.size() << ")";

    }//End loop over all wires

    // put the hit collection and associations into the event
    mf::LogInfo("TTHitFinder") << "Total TTHitFinder hits (U,V,Y)=(" 
                               << hitCollection_U.size() << ","
                               << hitCollection_V.size() << ","
                               << hitCollection_Y.size() << ")";
    hitCollection_U.put_into(evt); // reminder: instance name was "uhits"
    hitCollection_V.put_into(evt); // reminder: instance name was "vhits"
    hitCollection_Y.put_into(evt); // reminder: instance name was "yhits"

  } // End of produce()  
  

  //-------------------------------------------------
  float TTHitFinder::getTotalCharge(const float* signal_vector,int width=3,float threshold=-99){

    float totalCharge = 0;
    for(int tick=0; tick<width; tick++){
      if(signal_vector[tick] < threshold){
        totalCharge = -999; //special value for being below threshold
        break;
      }
      totalCharge += signal_vector[tick];
    }
    return totalCharge;
    
  }//end getTotalCharge method
  
  DEFINE_ART_MODULE(TTHitFinder)

} // end of hit namespace


#endif // TTHITFINDER_H