TruthSlim_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief  TruthSlim -- remove generated objects that don't contribute
/// \author dwwhitti@indiana.edu
////////////////////////////////////////////////////////////////////////
#include <cstdlib>
#include <string>
#include <sstream>
#include <sys/stat.h>
#include <vector>
#include <set>

// ROOT include
#include <TDatabasePDG.h>
#include <TSystem.h>
#include "TStopwatch.h"

// Framework includes
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Principal/Event.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "cetlib/search_path.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "Simulation/Particle.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "nusimdata/SimulationBase/MCFlux.h"
#include "nusimdata/SimulationBase/GTruth.h"
#include "Utilities/AssociationUtil.h"
#include "canvas/Persistency/Common/FindManyP.h"

namespace g4n {

  /// A module to slim the generator output
  class TruthSlim : public art::EDProducer {
  public:
    explicit TruthSlim(fhicl::ParameterSet const &pset);
    virtual ~TruthSlim();                        


    void produce (art::Event& evt);
    void beginRun(art::Run& run);

  private:

    TStopwatch          fStopwatch; ///< keep track of how long it takes to run the job
    std::string         fTruthLabel; ///< label used for the Monte Carlo generator module
    std::string         fGeantLabel; ///< label used for the Geant simulation module

  };

  //____________________________________________________________________________
  // Constructor
  TruthSlim::TruthSlim(fhicl::ParameterSet const& pset)
    : fTruthLabel(pset.get< std::string >("TruthLabel"))
    , fGeantLabel(pset.get< std::string >("GeantLabel"))
  {

    // Start counting the time
    fStopwatch.Start();

    // Produce new (slimmed) truth containers
    produces< std::vector<simb::MCTruth> >();
    produces< std::vector<simb::MCFlux>  >();
    produces< std::vector<simb::GTruth>  >();
    produces< art::Assns<simb::MCTruth, simb::MCFlux> >();
    produces< art::Assns<simb::MCTruth, simb::GTruth> >();

  }// end of constructor

  //____________________________________________________________________________
  // Destructor
  TruthSlim::~TruthSlim()  
  {  
    fStopwatch.Stop();
    mf::LogVerbatim("TruthSlim") << "real time to run truth slimming: " << fStopwatch.RealTime();
  }// end of destructor

  //______________________________________________________________________________
  void TruthSlim::beginRun(art::Run& run) 
  {
    return;
  }// end of begin job

  //______________________________________________________________________________
  void TruthSlim::produce(art::Event& evt)  
  {
    // Get the Particle vector.
    art::Handle< std::vector<sim::Particle> > particles;
    evt.getByLabel(fGeantLabel,particles);

    // Compile a list of MCTruth associated with the
    //   sim::Particles in this event
    std::set<art::Ptr<simb::MCTruth> > truthToKeep;
    // Retrieve the associated MCTruth(s)
    /* each element is a *list* of MCTruths associated with each particle */
    art::FindManyP<simb::MCTruth> truthAssn(particles,evt,fGeantLabel);
    // Loop through MCTruth objects to accumulate list of keepers
    for ( size_t p = 0; p < particles->size(); ++p ) {
      for ( size_t t = 0; t < truthAssn.at(p).size(); ++t ) {
        art::Ptr<simb::MCTruth> mcTru = truthAssn.at(p)[t];
        truthToKeep.insert(mcTru);
      }
    }

    // Retrieve the MCTruth information for this event
    std::vector< art::Handle< std::vector<simb::MCTruth> > > mclists;
    evt.getManyByType(mclists);

    // Create the associated (slimmed) truth data containers
    for(size_t mcl = 0; mcl < mclists.size(); ++mcl){
      std::unique_ptr<std::vector<simb::MCTruth>   >             slimMCTruth   (new std::vector<simb::MCTruth> );
      std::unique_ptr<std::vector<simb::MCFlux>    >             slimMCFlux    (new std::vector<simb::MCFlux>  );
      std::unique_ptr<std::vector<simb::GTruth>    >             slimGTruth    (new std::vector<simb::GTruth>  );
      std::unique_ptr<art::Assns<simb::MCTruth, simb::MCFlux > > slimFluxAssns (new art::Assns<simb::MCTruth, simb::MCFlux> );
      std::unique_ptr<art::Assns<simb::MCTruth, simb::GTruth > > slimGTruAssns (new art::Assns<simb::MCTruth, simb::GTruth> );
      // Loop through MCTruth container(s)
      art::Handle< std::vector<simb::MCTruth> > mclistHandle = mclists[mcl];
      // Get MCFlux and GTruth associations
      art::FindManyP<simb::MCFlux> fluxAssn(mclistHandle,evt,fTruthLabel);
      art::FindManyP<simb::GTruth> gtruAssn(mclistHandle,evt,fTruthLabel);

      // Loop through MCTruth objects to veto those with no associated sim::Particles.
      for(size_t i = 0; i < mclistHandle->size(); ++i){
        art::Ptr<simb::MCTruth> mct(mclistHandle, i);
        // Skip this MCTruth if it's not in the list of keepers
        //   (Only works if fSlimSimParticles was true)
        if ( truthToKeep.find(mct) == truthToKeep.end() ) continue;
        // Pack everything into slimmed containers and build new associations
        slimMCTruth->push_back(*mct);
        for ( size_t j = 0; j < fluxAssn.at(i).size(); ++j ) {
          art::Ptr<simb::MCFlux> theFlux = fluxAssn.at(i)[j];
          slimMCFlux->push_back(*theFlux);
          util::CreateAssn(*this, evt, *slimMCTruth, theFlux, *slimFluxAssns, i);
        }
        for ( size_t j = 0; j < gtruAssn.at(i).size(); ++j ) {
          art::Ptr<simb::GTruth> theGTru = gtruAssn.at(i)[j];
          slimGTruth->push_back(*theGTru);
          util::CreateAssn(*this, evt, *slimMCTruth, theGTru, *slimGTruAssns, i);
        }
      }
      evt.put(std::move(slimMCTruth));
      evt.put(std::move(slimMCFlux ));
      evt.put(std::move(slimGTruth ));
      evt.put(std::move(slimFluxAssns));
      evt.put(std::move(slimGTruAssns));
    }// end loop through MCTruth container(s)

    return;
  }// end of produce

} // end namespace

DEFINE_ART_MODULE(g4n::TruthSlim)