MichelHitCounter_module.cc

Go to the documentation of this file.

//  Example module showing how to access and count
//  michel elctron hits in reconstruction

#include "lardataobj/Simulation/SimChannel.h"
#include "larsim/Simulation/LArG4Parameters.h"
#include "larsim/Simulation/LArVoxelData.h"
#include "larsim/Simulation/LArVoxelList.h"
#include "larsim/Simulation/SimListUtils.h"
#include "larsim/MCCheater/BackTracker.h"
#include "larcore/Geometry/Geometry.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "larreco/Calorimetry/CalorimetryAlg.h"
#include "larreco/RecoAlg/PMAlg/Utilities.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h"
#include "lardata/Utilities/DatabaseUtil.h"

#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/FindManyP.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 "art/Framework/Services/Registry/ServiceHandle.h"
#include "art_root_io/TFileService.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "lardata/ArtDataHelper/MVAReader.h"

#include "TTree.h"

#include <fstream>

namespace MichelReco {

class MichelReco;

class MichelReco : public art::EDAnalyzer {
public:
  explicit MichelReco(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.
  MichelReco(MichelReco const &) = delete;
  MichelReco(MichelReco &&) = delete;
  MichelReco & operator = (MichelReco const &) = delete;
  MichelReco & operator = (MichelReco &&) = delete;

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

  // Selected optional functions
  void beginJob() override;
  void endJob() override;
  void beginRun(const art::Run& run) override;
  void reconfigure(fhicl::ParameterSet const& p) ;

private:
  // My parameters
  size_t fEvNumber;   // Event number for each event
  int fNMichelHits;   // Total number of Michel hits in events
  
  // Input parameters
  art::InputTag fNNetModuleLabel;   // label of the module used for CNN tagging
};


MichelReco::MichelReco(fhicl::ParameterSet const & p)
  : 
  EDAnalyzer(p)
{
  reconfigure(p);
}

void MichelReco::beginJob()
{
  art::ServiceHandle<art::TFileService> tfs;
}

void MichelReco::endJob()
{
  mf::LogVerbatim("MichelReco") << "MichelHitCounter finished job";
  mf::LogVerbatim("MichelReco") << fNMichelHits << " michel hits found";
}

void MichelReco::beginRun(const art::Run&)
{
  art::ServiceHandle<sim::LArG4Parameters> larParameters;
}

void MichelReco::analyze(art::Event const & evt)
{
  // keep track of progress
  fEvNumber = evt.id().event();
  mf::LogVerbatim("MichelReco") << "MichelReco module on event " << fEvNumber;
  
  // Helper to get hits and the 4 associated CNN outputs
  // CNN Outputs: EM, Track, Michel, Empty
  anab::MVAReader<recob::Hit,4> hitResults(evt, fNNetModuleLabel);
  
  // loop over hits 
  for (size_t h = 0; h < hitResults.size(); ++h) {
    
    // Get cnn output for hit h 
    std::array<float,4> cnn_out = hitResults.getOutput(h);

    // Compare michel cnn output to treshold to decide if Michel hit 
    // Example using 0.5, needs to be optimised with clustering
    // e.g. 0.9 used for event selection
    if (cnn_out[hitResults.getIndex("michel")] > 0.5) {
      fNMichelHits += 1;
    }

  } // End of loop over hits

}

void MichelReco::reconfigure(fhicl::ParameterSet const& p)
{
  
  fNMichelHits = 0; 
  fNNetModuleLabel = p.get<std::string>("NNetModuleLabel");
  return;
}

} // MichelReco namespace

DEFINE_ART_MODULE(MichelReco::MichelReco)