HitHandScan_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       HitHandScan
// Module Type: producer
// File:        HitHandScan_module.cc
//
// Generated at Tue May  3 01:14:23 2016 by Matthew Thiesse using artmod
// from cetpkgsupport v1_10_01.
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "canvas/Persistency/Common/FindOneP.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Wire.h"
#include "lardataobj/RawData/RawDigit.h"
#include "lardata/ArtDataHelper/HitCreator.h"
#include "lardata/Utilities/AssociationUtil.h"

#include <sstream>
#include <string>
#include <fstream>
#include <vector>
#include <memory>
#include <map>

#include "TMath.h"

class HitHandScan;

typedef struct {
  unsigned int run;
  unsigned int event;
  unsigned int tpc;
  unsigned int plane;
  double wire;// wire number in this tpc on this plane
  double tick;// IN RELATION TO EXTERNAL TRIGGER, NOT NECESSARILY THE FIRST TICK IN THE EVENT
  unsigned int numpoly;
} Point;

class HitHandScan : public art::EDProducer {
public:
  explicit HitHandScan(fhicl::ParameterSet const & p);
  // The destructor generated by the compiler is fine for classes
  // without bare pointers or other resource use.

  // Plugins should not be copied or assigned.
  HitHandScan(HitHandScan const &) = delete;
  HitHandScan(HitHandScan &&) = delete;
  HitHandScan & operator = (HitHandScan const &) = delete;
  HitHandScan & operator = (HitHandScan &&) = delete;

  // Required functions.
  void produce(art::Event & e) override;

  // Selected optional functions.
  void reconfigure(fhicl::ParameterSet const & p) ;

private:
  bool insidePolygon(std::vector<Point> polygon, Point ptest);

  std::string fPreviousHitModuleLabel;
  std::string fHandScanFileName;

};


HitHandScan::HitHandScan(fhicl::ParameterSet const & p) : EDProducer{p}
{
  this->reconfigure(p);
  recob::HitCollectionCreator::declare_products(producesCollector());
}

void HitHandScan::produce(art::Event & e)
{
  auto prevHitHandle = e.getHandle< std::vector< recob::Hit> >(fPreviousHitModuleLabel);

  art::FindOneP<recob::Wire> wires(prevHitHandle,e,fPreviousHitModuleLabel);
  art::FindOneP<raw::RawDigit> rawdigits(prevHitHandle,e,fPreviousHitModuleLabel);

  recob::HitCollectionCreator hcol(e, wires.isValid(), rawdigits.isValid());

  std::vector<Point> trackBox;
  
  std::string readrun,readevent,readtpc,readplane,readwire,readtick,readnumpoly;
  unsigned int run,event,tpc,plane,numpoly;
  double tick,wire;
  int linenumber = 0;
  std::string line;
  std::ifstream handscanfile(fHandScanFileName.c_str(),std::ifstream::in);
  if (handscanfile.is_open())
    {
      while (std::getline(handscanfile,line))
        {
          linenumber++;
          std::stringstream ss(line);
          ss >> std::ws >> readrun >> readevent >> readtpc >> readplane >> readwire >> readtick >> readnumpoly;
          try 
            {
              run = stoul(readrun);
              event = stoul(readevent);
              tick = stod(readtick);
              tpc = stoul(readtpc);
              plane = stoul(readplane);
              wire = stod(readwire);
              numpoly = stoul(readnumpoly);
              
              Point p;
              p.run = run;
              p.event = event;
              p.tpc = tpc;
              p.plane = plane;
              p.wire = wire;
              p.tick = tick;
              p.numpoly = numpoly;
              trackBox.push_back(p);
            } 
          catch (const std::exception& ex)
            {
              std::cout << "Line # " << linenumber << " is invalid, skipping..." << std::endl;
            }
        }
      handscanfile.close();
    }
  else
    {
      std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
                << "NO FILE FOUND WITH HAND SCANNED HIT DATA\n"
                << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
      return;
    }

  for (size_t i_hit = 0; i_hit < prevHitHandle->size(); i_hit++)
    {
      art::Ptr<recob::Hit> phit(prevHitHandle,i_hit);

      art::Ptr<recob::Wire> wire(wires.at(i_hit));
      art::Ptr<raw::RawDigit> rawdigit(rawdigits.at(i_hit));

      Point ptest;
      ptest.run = e.run();
      ptest.event = e.event();
      ptest.tpc = phit->WireID().TPC;
      ptest.plane = phit->WireID().Plane;
      ptest.wire = phit->WireID().Wire;
      ptest.tick = phit->PeakTime();
            
      std::map< unsigned int, std::vector<Point> > polygons;
      for (auto const &poly : trackBox)
        {
          if (poly.run == ptest.run 
              && poly.event == ptest.event
              && poly.tpc == ptest.tpc
              && poly.plane == ptest.plane)
            {
              polygons[poly.numpoly].push_back(poly);
            }
        }
      bool acceptPoint = false;
      for (auto const &p : polygons)
        {
          if (insidePolygon(p.second,ptest)) 
            {
              acceptPoint = true;
            }
        }
      if (acceptPoint) 
        {
          hcol.emplace_back(*phit,wire,rawdigit);
        }
    }
  
  hcol.put_into(e);

}

void HitHandScan::reconfigure(fhicl::ParameterSet const & p)
{
  fHandScanFileName = p.get<std::string>("HandScanFileName");
  fPreviousHitModuleLabel = p.get<std::string>("PreviousHitModuleLabel");
}

bool HitHandScan::insidePolygon(std::vector<Point> polygon, Point ptest)
{
  size_t n = polygon.size();
  size_t i, j;
  bool c = false;
  for (i = 0, j = n - 1; i < n; j = i++)
    {
      if ( ((polygon[i].tick > ptest.tick) != (polygon[j].tick > ptest.tick)) &&
           (ptest.wire < (polygon[j].wire - polygon[i].wire) * 
            (ptest.tick - polygon[i].tick) / (polygon[j].tick - polygon[i].tick) + polygon[i].wire) )
        c = !c;
    }
  return c;
}

DEFINE_ART_MODULE(HitHandScan)