FinalStateParticleFilter_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
//
// FinalStateParticleFilter class:
// Algoritm to produce a filtered event file having
// events with user-defined final state particles
//
// saima@ksu.edu
//
////////////////////////////////////////////////////////////////////////

//Framework Includes
#include "art/Framework/Core/EDFilter.h"
#include "art/Framework/Core/ModuleMacros.h" 
#include "art/Framework/Principal/Event.h" 
#include "fhiclcpp/ParameterSet.h" 
#include "art/Framework/Principal/Handle.h" 
#include "canvas/Persistency/Common/Ptr.h" 
#include "art/Framework/Services/Registry/ServiceHandle.h" 
#include "art_root_io/TFileService.h"

//Larsoft Includes
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"

// ROOT includes
#include "TH1.h"

namespace filt {

  class FinalStateParticleFilter : public art::EDFilter  {

  public:

    explicit FinalStateParticleFilter(fhicl::ParameterSet const& );

    bool filter(art::Event& evt);
    void beginJob();


  private:

    std::string fGenieModuleLabel;
    std::vector<int> fPDG;
    std::vector<int> fStatusCode;
    TH1D* fSelectedEvents;
    TH1D* fTotalEvents;

    bool isSubset(std::vector<int> const& a, std::vector<int> const& b) const;

  }; // class FinalStateParticleFilter

}

namespace filt{

  //-------------------------------------------------
  FinalStateParticleFilter::FinalStateParticleFilter(fhicl::ParameterSet const & pset)
    : EDFilter{pset}
  {
    fGenieModuleLabel = pset.get< std::string      >("GenieModuleLabel");
    fPDG              = pset.get< std::vector<int> >("PDG");
  }

  //-------------------------------------------------
  void FinalStateParticleFilter::beginJob()
  {
    art::ServiceHandle<art::TFileService const> tfs;
    fSelectedEvents = tfs->make<TH1D>("fSelectedEvents", "Number of Selected Events", 3, 0, 3); //counts the number of selected events
    fTotalEvents = tfs->make<TH1D>("fTotalEvents", "Total Events", 3, 0, 3); //counts the initial number of events in the unfiltered root input file
  }

  //-------------------------------------------------
  bool FinalStateParticleFilter::filter(art::Event &evt)
  {

    //const TDatabasePDG* databasePDG = TDatabasePDG::Instance();

    art::Handle< std::vector<simb::MCTruth> > mclist;
    evt.getByLabel(fGenieModuleLabel,mclist);
    art::Ptr<simb::MCTruth> mc(mclist,0);

    fTotalEvents->Fill(1);

    std::vector<int> finalstateparticles;

    //get a vector of final state particles
    for(int i = 0; i < mc->NParticles(); ++i){
      simb::MCParticle part(mc->GetParticle(i));
      if(part.StatusCode()== 1)
        finalstateparticles.push_back(part.PdgCode());
    }

    if(isSubset(fPDG, finalstateparticles)){
      fSelectedEvents->Fill(1);
      std::cout << "this is a selected event" << std::endl;
    }

    return isSubset(fPDG, finalstateparticles); // returns true if the user-defined fPDG exist(s) in the final state particles

  } // bool
  //} // namespace

//------------------------------------------------

bool FinalStateParticleFilter::isSubset(std::vector<int> const& a, std::vector<int> const& b) const
{
  for (auto const a_int : a) {
    bool found = false;
    for (auto const b_int : b) {
      if (a_int == b_int){
        found = true;
        break;
      }
    }

    if (!found){
      return false;
    }
  }
  return true;
}
}

//--------------------------------------------------

namespace filt {

  DEFINE_ART_MODULE(FinalStateParticleFilter)

} //namespace filt