GENIEFilter_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       GENIEFilter
// Module Type: filter
// File:        GENIEFilter_module.cc
//
// Generated at Tue Nov 17 11:20:20 2015 by Kanika Sachdev using artmod
// from cetpkgsupport v1_08_07.
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDFilter.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 <memory>

#include "GENIE/Framework/GHEP/GHepStatus.h"

#include "nusimdata/SimulationBase/MCTruth.h"

namespace evgen {
  
  /// A filter module to throw away spills where genie output
  /// has undesired characteristics. Currently only throws away
  /// events with undecayed 2-nucleon final states
  class GENIEFilter;
}

class evgen::GENIEFilter : public art::EDFilter {
public:
  explicit GENIEFilter(fhicl::ParameterSet const & p);

  GENIEFilter(GENIEFilter const &) = delete;
  GENIEFilter(GENIEFilter &&) = delete;
  GENIEFilter & operator = (GENIEFilter const &) = delete;
  GENIEFilter & operator = (GENIEFilter &&) = delete;


  bool filter(art::Event & e) override;
  void reconfigure(const fhicl::ParameterSet& pset);

private:
  const art::ProductToken<std::vector<simb::MCTruth>> fGENIEToken;

};

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

evgen::GENIEFilter::GENIEFilter(fhicl::ParameterSet const & p)
  : fGENIEToken(consumes<std::vector<simb::MCTruth>>(p.get<std::string>("GENIELabel")))
{
  reconfigure(p);
}

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

void evgen::GENIEFilter::reconfigure(const fhicl::ParameterSet& pset)
{
}

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

bool evgen::GENIEFilter::filter(art::Event & e)
{
  art::Handle< std::vector<simb::MCTruth> > truthcol;
  e.getByToken(fGENIEToken, truthcol);
  
  unsigned int ntruths = truthcol->size();

  for(unsigned int itruth = 0; itruth < ntruths; ++itruth){
    
    unsigned int nparts = truthcol->at(itruth).NParticles();

    for(unsigned int ipart = 0; ipart < nparts; ++ipart){
      const simb::MCParticle&  particle= truthcol->at(itruth).GetParticle(ipart);

      // Criteria to throw away dinucleon states 
      // that GENIE wasn't able to decay into appropriate
      // final states. From Robert H., doc 14388, dated Nov 16, 2015.
      // Needed for genie 2_8_6. Might not be needed for future
      // revisions.
      if( (particle.PdgCode() == 2000000200  ||  // n+n final state
           particle.PdgCode() == 2000000201  ||  // n+p final state
           particle.PdgCode() == 2000000202) &&  // p+p final state
          (particle.StatusCode() == genie::kIStIntermediateState ||
           particle.StatusCode() == genie::kIStDecayedState )){
        
        // loop over all the other particles and see if anyone 
        // claims to have this dinucleon as its mother
        // We can directly probe the daughter info directly because 
        // GENIE doesn't seem to be filling it in properly
        bool foundDaughter = false;

        for(unsigned int jpart = 0; jpart < nparts; ++jpart){
          const simb::MCParticle&  jparticle= truthcol->at(itruth).GetParticle(jpart);
          if( jparticle.Mother() == particle.TrackId() ){
            foundDaughter = true;
            break;
          }
        }// end loop over j particles

        if(!foundDaughter)
          return Fail;
        
      }// end if found a dinucleon
    }// end loop over MCParticles
  }// end loop over MCTruths
  
  return Pass;
}

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

DEFINE_ART_MODULE(evgen::GENIEFilter)