MCHitFinder_module.cc

Go to the documentation of this file.

// Class:       MCHitFinder
// Module Type: producer
// File:        MCHitFinder_module.cc
//
// Generated at Tue Jun 24 07:30:08 2014 by Kazuhiro Terao using artmod
// from cetpkgsupport v1_05_04.

#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 <memory>

#include "lardataobj/Simulation/SimChannel.h"

#include "larcore/Geometry/Geometry.h"
#include "lardataobj/MCBase/MCDataHolder.h"
#include "lardataobj/MCBase/MCHitCollection.h"
#include "lardataobj/MCBase/MCWireCollection.h"

namespace hit {

  class MCHitFinder;
  
  class MCHitFinder : public art::EDProducer {
  public:
    explicit MCHitFinder(fhicl::ParameterSet const & p);
    virtual ~MCHitFinder();
    
    void produce(art::Event & e) override;
    
  private:
    
    std::string fLArG4ModuleName;
    bool fVerbose;
    bool fMakeMCWire;
    
  };
  
  
  MCHitFinder::MCHitFinder(fhicl::ParameterSet const & p)
  {
    fLArG4ModuleName = p.get<std::string>("LArG4ModuleName");
    fMakeMCWire      = p.get<bool>("MakeMCWire");
    fVerbose         = p.get<bool>("Verbose");
    produces< std::vector< sim::MCHitCollection> >();
    if(fMakeMCWire)
      produces< std::vector< sim::MCWireCollection> >();
  }
  
  MCHitFinder::~MCHitFinder()
  {
  }
  
  void MCHitFinder::produce(art::Event & e)
  {
    
    art::ServiceHandle<geo::Geometry> geo;
    
    const unsigned int nch = geo->Nchannels();
    
    std::unique_ptr< std::vector<sim::MCHitCollection> > hits_v ( new std::vector<sim::MCHitCollection>() );
    std::unique_ptr< std::vector<sim::MCWireCollection> > wires_v ( new std::vector<sim::MCWireCollection>() );
    
    hits_v->reserve(nch);
    wires_v->reserve(nch);
    for(size_t ch=0; ch<nch; ++ch) {
      
      hits_v->push_back(sim::MCHitCollection(ch));
      wires_v->push_back(sim::MCWireCollection(ch));
      
    }
    
    art::Handle< std::vector<sim::SimChannel> > simchArray;
    e.getByLabel(fLArG4ModuleName,simchArray);
    if(!simchArray.isValid())
      throw cet::exception(__PRETTY_FUNCTION__) 
        << "Did not find sim::SimChannel with a label: " << fLArG4ModuleName.c_str() << std::endl;
    
    // Loop over SimChannel
    for(size_t simch_index=0; simch_index<simchArray->size(); ++simch_index) {
      
      const art::Ptr<sim::SimChannel> simch_ptr(simchArray,simch_index);
      
      size_t ch = simch_ptr->Channel();
      
      if(ch >= hits_v->size())
        
        throw cet::exception(__PRETTY_FUNCTION__)
          << "Channel number " << ch << " exceeds total # of channels: " << nch << std::endl;
      
      auto &mchits  = hits_v->at(ch);
      auto &mcwires = wires_v->at(ch);
      
      std::map<sim::MCEnDep,sim::MCWire> edep_wire_map;
      auto tdc_ide_map = simch_ptr->TDCIDEMap();
      
      // Loop over stored data for a particular sim::SimChannel //
      
      //
      // Read & Convert
      //
      if(fVerbose)
        std::cout<<std::endl<<"Processing Ch: "<<ch<<std::endl;
      
      for(auto const& tdc_ide_pair : tdc_ide_map) {
        
        auto const& tdc   = tdc_ide_pair.first;
        auto const& ide_v = tdc_ide_pair.second;
        
        for(auto const& ide : ide_v) {
          sim::MCEnDep edep;
          edep.SetVertex(ide.x,ide.y,ide.z);
          edep.SetEnergy(ide.energy);
          edep.SetTrackId(ide.trackID);
          
          sim::MCWire wire;
          wire.SetStartTDC(tdc);
          
          auto edep_iter = edep_wire_map.insert(std::make_pair(edep,wire));
          
          //auto edep_iter = edep_wire_map.find(edep);
          
          auto last_tdc = (edep_iter).first->second.StartTDC() + (edep_iter).first->second.size() - 1;
          
          if(fVerbose) {
            
            if( edep_iter.second ) std::cout<<std::endl;
            
            std::cout<<"  Track: "<<ide.trackID
                     <<" Vtx: " <<ide.x << " " << ide.y << " " <<ide.z << " " <<ide.energy
                     << " ...  @ TDC = "<<tdc<<" ... "<<ide.numElectrons << std::endl;
          }
          
          if( !(edep_iter.second) ) {
            
            if( last_tdc+1 != tdc ) {
              /*
              std::cerr 
                << "Found discontinuous TDC! "
                << " Last (ADC @ TDC): " << (*(edep_iter.first->second.rbegin())) << " @ " << last_tdc 
                << " while current (ADC @ TDC): " << (ide.numElectrons * detp->ElectronsToADC()) << " @ " << tdc
                << " ... skipping to store!"
                << std::endl;
              */
              continue;
            }
          }
          (edep_iter).first->second.push_back(ide.numElectrons);
        }// Done looping over IDEs in a particular TDC
      }// Done looping over TDC
      
      // 
      // Store
      //
      for(auto const& edep_wire_pair : edep_wire_map) {
        
        auto const& edep = edep_wire_pair.first;
        auto const& wire = edep_wire_pair.second;
        
        // Create MCHit
        ::sim::MCHit hit;
        float vtx[3] = {float(edep.Vertex()[0]),
                        float(edep.Vertex()[1]),
                        float(edep.Vertex()[2])};
        hit.SetParticleInfo(vtx, edep.Energy(), edep.TrackId());
        
        double max_time = 0;
        double qsum     = 0;
        double qmax     = 0;
        for(size_t wire_index=0; wire_index < wire.size(); ++wire_index) {
          
          auto q = wire.at(wire_index);
          
          qsum += q;
          
          if( q > qmax) { qmax = q; max_time = wire.StartTDC() + wire_index; }
          
        }       
        hit.SetCharge(qsum,qmax);
        hit.SetTime(max_time,0);
        
        mchits.push_back(hit);
        
        if(!fMakeMCWire) continue;
        
      mcwires.push_back(wire);
      } // End looping over all MCEnDep-MCWire pairs
    } // End looping over all SimChannels

    e.removeCachedProduct(simchArray);
    
    std::sort((*hits_v).begin(),(*hits_v).end());
    e.put(std::move(hits_v));
    
    if(fMakeMCWire) {
      std::sort((*wires_v).begin(),(*wires_v).end());
      e.put(std::move(wires_v));
    }
    
  }
}
DEFINE_ART_MODULE(hit::MCHitFinder)