GenieOutput_module.cc

Go to the documentation of this file.

#include "art/Framework/Core/EDAnalyzer.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 "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "nusimdata/SimulationBase/MCTruth.h"
#include "nusimdata/SimulationBase/MCFlux.h"
#include "nusimdata/SimulationBase/GTruth.h"

#include "nugen/EventGeneratorBase/GENIE/GENIE2ART.h"
#include "nugen/EventGeneratorBase/GENIE/MCTruthAndFriendsItr.h"

// GENIE includes
#ifdef GENIE_PRE_R3
  #include "GENIE/Messenger/Messenger.h"
  #include "Ntuple/NtpMCFormat.h"
  #include "Ntuple/NtpWriter.h"
  #include "Ntuple/NtpMCEventRecord.h"
  //#include "Ntuple/NtpMCTreeHeader.h"
  #include "PDG/PDGLibrary.h"
  #include "GHEP/GHepRecord.h"
#else
  #include "GENIE/Framework/Messenger/Messenger.h"
  // careful: potential conflict LOG_INFO w/ messagefacility
  #include "GENIE/Framework/GHEP/GHepRecord.h"
  #include "GENIE/Framework/Ntuple/NtpMCFormat.h"
  #include "GENIE/Framework/Ntuple/NtpWriter.h"
  #include "GENIE/Framework/Ntuple/NtpMCEventRecord.h"
  // #include "GENIE/Framework/Ntuple/NtpMCTreeHeader.h"
#endif

#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>

// Necessary because the GENIE LOG_* macros don't fully qualify Messenger
using genie::Messenger;

namespace evg {
  class GenieOutput;

  struct GenieOutputParams {
    // hold/document fcl parameters
    template<class T> using Atom     = fhicl::Atom<T>;
    template<class T> using Sequence = fhicl::Sequence<T>;
    template<class T> using Table    = fhicl::Table<T>;
    using Comment = fhicl::Comment;
    using Name    = fhicl::Name;

    Sequence<std::string> inputModuleLabels {
      Name("inputModuleLabels"),
      Comment("list of input module labels to use, if empty use everything")
        // { }  // how to default this to empty ??
    };
    Atom<std::string> outputGHEPFilePattern {
      Name("outputGHEPFile"),
      Comment("name of file to write in std GENIE gntp.*.ghep.root format\n"
              "if blank don't write file, if name contains '%l' then\n"
              "events from different input labels will go to separate outputs\n"),
      ""
    };
    Atom<std::string> dumpFilePattern {
      Name("dumpFilePattern"),
      Comment("name of file for formatted dumps; if name contains '%l' then\n"
              "events from different input labels will go to separate streams\n"
              "if blank use std::cout"),
      ""
    };
    Atom<int>  dumpGeniePrintLevel {
      Name("dumpGeniePrintLevel"),
      Comment("print fetched genie::EventRecord -1=no, 13=max info\n"
              "see GENIE manual for legal values"),
      -1
    };
    Atom<bool> dumpMCTruth {
      Name("dumpMCTruth"),
      Comment("dump the MCTruth objects (std:cout) as they're retrieved"),
      false
    };
    Atom<bool> dumpGTruth {
      Name("dumpGTruth"),
      Comment("dump the GTruth objects (std:cout) as they're retrieved"),
      false
    };
    Atom<bool> dumpMCFlux {
      Name("dumpMCFlux"),
      Comment("dump the MCFlux objects (std:cout) as they're retrieved"),
      false
    };
  };  // end-of GenieOutputParams
}

class evg::GenieOutput : public art::EDAnalyzer {

public:

  // Allow 'art --print-description' to work
  using Parameters = art::EDAnalyzer::Table<evg::GenieOutputParams>;

  //explicit GenieOutput(fhicl::ParameterSet const & p);
  explicit GenieOutput(const Parameters& params);

  // The destructor generated by the compiler is fine for classes
  // without bare pointers or other resource use.
  ~GenieOutput();

  // Plugins should not be copied or assigned.
  GenieOutput(GenieOutput const &) = delete;
  GenieOutput(GenieOutput &&) = delete;
  GenieOutput & operator = (GenieOutput const &) = delete;
  GenieOutput & operator = (GenieOutput &&) = delete;

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

  // Selected optional functions.
  /*
  void beginJob() override;
  void beginRun(art::Run const & r) override;
  void beginSubRun(art::SubRun const & sr) override;
  void endJob() override;
  void endRun(art::Run const & r) override;
  void endSubRun(art::SubRun const & sr) override;
  void reconfigure(fhicl::ParameterSet const & p) override;
  void respondToCloseInputFile(art::FileBlock const & fb) override;
  void respondToCloseOutputFiles(art::FileBlock const & fb) override;
  void respondToOpenInputFile(art::FileBlock const & fb) override;
  void respondToOpenOutputFiles(art::FileBlock const & fb) override;
  */

private:

  // private methods
  genie::NtpWriter*          FetchNtpWriter(const std::string& label);
  std::ostream*              FetchDumpStream(const std::string& label);


  // member data here.

  Parameters                 fParams;

  std::vector<std::string>   fInputModuleLabels;  
  std::string                fOutputGHEPFilePattern;
  std::string                fDumpFilePattern;
  int                        fDumpGeniePrintLevel;

  bool                       fDumpMCTruth;
  bool                       fDumpGTruth;
  bool                       fDumpMCFlux;

  bool                       fSeparateOutputNtpWriters;
  bool                       fSeparateDumpStreams;

  std::map<std::string,genie::NtpWriter*>  fOutputNtpWriters;
  std::map<std::string,std::ostream*>      fDumpStreams;

};


//evg::GenieOutput::GenieOutput(fhicl::ParameterSet const & pset)
// : EDAnalyzer(pset)
evg::GenieOutput::GenieOutput(const Parameters& params)
  : EDAnalyzer(params)
  , fParams(params)
  , fSeparateOutputNtpWriters(false)
  , fSeparateDumpStreams(false)
{
#ifdef GENIE_PRE_R3
  // trigger early initialization of PDG database & GENIE message service
  // just to get it out of the way and not intermixed with other output
  genie::PDGLibrary::Instance();
#else
  // get the GENIE banner out of the way
  // no longer can use genie::PDGLibrary::Instance() to do this
  // because that must happen, in some cases in v3_02_xx, after the tune
  // is determined
  // banner is triggered by first use of GENIE Messenger
  // avoid using GENIE macros (possible conflict with mf macros)
  // LOG("GENIE",pInfo) << "Trigger GENIE banner";
  (*genie::Messenger::Instance())("GENIE") << log4cpp::Priority::INFO
                                           << "Trigger GENIE banner";
#endif


  fInputModuleLabels     = fParams().inputModuleLabels();
  fOutputGHEPFilePattern = fParams().outputGHEPFilePattern();
  fDumpFilePattern       = fParams().dumpFilePattern();
  fDumpGeniePrintLevel   = fParams().dumpGeniePrintLevel();

  fSeparateOutputNtpWriters =
    ( fOutputGHEPFilePattern.find("%l") != std::string::npos );
  fSeparateDumpStreams =
    ( fDumpFilePattern.find("%l") != std::string::npos );

  fDumpMCTruth           = fParams().dumpMCTruth();
  fDumpGTruth            = fParams().dumpGTruth();
  fDumpMCFlux            = fParams().dumpMCFlux();

  /*
  mf::LogInfo("GenieOutput") << "##### Dump options "
                             << fDumpMCTruth << " "
                             << fDumpGTruth << " "
                             << fDumpMCFlux << " " ;
  */

}

evg::GenieOutput::~GenieOutput()
{

  // release resources
  std::map<std::string,genie::NtpWriter*>::iterator mitro =
    fOutputNtpWriters.begin();
  for ( ; mitro != fOutputNtpWriters.end(); ++mitro ) {
    std::string       label = mitro->first;
    genie::NtpWriter*  ntpw = mitro->second;
    if ( ntpw ) {
      // close out the file
      ntpw->Save();
      delete ntpw;
    }
    mitro->second = 0;
  }

  std::map<std::string,std::ostream*>::iterator mitrd = fDumpStreams.begin();
  for ( ; mitrd != fDumpStreams.end(); ++mitrd ) {
    std::string    label = mitrd->first;
    std::ofstream* fout = dynamic_cast<std::ofstream*>(mitrd->second);
    if ( fout ) {
      fout->flush();
      fout->close();
      delete fout;
    }
    mitrd->second = 0;
  }
}

void evg::GenieOutput::analyze(art::Event const & evt)
{
  //std::cerr << "evg::GenieOutput::analyze " << std::endl;

  evgb::MCTruthAndFriendsItr mcitr(evt,fInputModuleLabels);

  bool flag = true;
  int  indx = -1;
  while ( ( flag = mcitr.Next() ) ) {
    ++indx;
    std::string label             = mcitr.GetLabel();
    const simb::MCTruth* pmctruth = mcitr.GetMCTruth();
    const simb::GTruth*  pgtruth  = mcitr.GetGTruth();
    const simb::MCFlux*  pmcflux  = mcitr.GetMCFlux();

    static int ievt = -1; // ! bah!
    ++ievt;

    static const simb::GTruth nullGTruth;
    if ( ! pgtruth ) {
      mf::LogInfo("GenieOutput") << "##### no GTruth ";
      pgtruth = &nullGTruth;
    }

    genie::EventRecord* grec = evgb::RetrieveGHEP(*pmctruth,*pgtruth);

    genie::NtpWriter* ntpWriter = FetchNtpWriter(label);
    if ( ntpWriter ) {
      // check ownership! (might need copy)
      ntpWriter->AddEventRecord(ievt,grec);
    }

    std::ostream* osdump = FetchDumpStream(label);
    if ( osdump ) {
      *osdump
        << " ** Event: GenieOutput_module " << ievt
        << *grec;
      osdump->flush();
    }

    if ( fDumpMCTruth || fDumpGTruth || fDumpMCFlux ) {
      std::ostringstream dumpSimBaseObj;
      dumpSimBaseObj << " after Next() " << indx << " " << flag
                     << std::endl;;
      if ( fDumpMCTruth ) {
        if ( pmctruth ) dumpSimBaseObj << *pmctruth << std::endl;
        else            dumpSimBaseObj << "no simb::MCTruth available" << std::endl;
      }
      if ( fDumpGTruth ) {
        if ( pgtruth ) {
          dumpSimBaseObj << *pgtruth << std::endl;
          //dumpSimBaseObj << "sorry no operator<< exists for simb::GTruth"
          //               << " - someone should write one" << std::endl;
        } else dumpSimBaseObj << "no simb::GTruth available" << std::endl;
      }
      if ( fDumpMCFlux ) {
        if ( pmcflux ) dumpSimBaseObj << *pmcflux << std::endl;
        else           dumpSimBaseObj << "no simb::MCFlux available" << std::endl;
      }
      mf::LogInfo("GenieOutput") << dumpSimBaseObj.str();
    }

    delete grec;  // don't leak stuff
  } // loop over MCTruthAndFriends

}

genie::NtpWriter* evg::GenieOutput::FetchNtpWriter(const std::string& label) {

  if ( fOutputGHEPFilePattern == "" ) return 0;

  genie::NtpWriter* ntpwret = 0;
  if ( fSeparateOutputNtpWriters ) ntpwret = fOutputNtpWriters[label];
  else                             ntpwret = fOutputNtpWriters["*"];

  if ( ntpwret ) return ntpwret; // already openned

  // nope?? okay

  std::string finalFileName = fOutputGHEPFilePattern;
  if ( fSeparateOutputNtpWriters ) {
    size_t posl = finalFileName.find("%l");
    if ( posl != std::string::npos ) {
      finalFileName.replace(posl,2,label);
    }
  }
  ntpwret = new genie::NtpWriter(genie::kNFGHEP,0);
  ntpwret->CustomizeFilename(finalFileName);
  ntpwret->Initialize();

  if ( fSeparateOutputNtpWriters ) fOutputNtpWriters[label] = ntpwret;
  else                             fOutputNtpWriters["*"]   = ntpwret;

  return ntpwret;


}


std::ostream* evg::GenieOutput::FetchDumpStream(const std::string& label) {
  if (fDumpGeniePrintLevel < 0 ) return 0;
  genie::GHepRecord::SetPrintLevel(fDumpGeniePrintLevel);

  std::ostream* osret = 0;
  if ( fSeparateDumpStreams ) osret = fDumpStreams[label];
  else                        osret = fDumpStreams["*"];

  if ( osret ) return osret; // already openned

  // nope?? okay

  if ( fDumpFilePattern == ""          ||
       fDumpFilePattern == "--"        ||
       fDumpFilePattern == "cout"      ||
       fDumpFilePattern == "std::cout"     ) {
    // standardize so we don't check all these again
    fDumpFilePattern = "std::cout";
    osret = &(std::cout);
  } else {
    std::string finalFileName = fDumpFilePattern;
    if ( fSeparateDumpStreams ) {
      size_t posl = finalFileName.find("%l");
      if ( posl != std::string::npos ) {
        finalFileName.replace(posl,2,label);
      }
    }
    osret = new std::ofstream(finalFileName.c_str(),
                              std::ios_base::trunc|std::ios_base::out);

  }
  if ( fSeparateDumpStreams ) fDumpStreams[label] = osret;
  else                        fDumpStreams["*"]   = osret;

  return osret;

}

/*
void evg::GenieOutput::beginJob()
{
  // Implementation of optional member function here.
}

void evg::GenieOutput::beginRun(art::Run const & r)
{
  // Implementation of optional member function here.
}

void evg::GenieOutput::beginSubRun(art::SubRun const & sr)
{
  // Implementation of optional member function here.
}

void evg::GenieOutput::endJob()
{
  // Implementation of optional member function here.
}

void evg::GenieOutput::endRun(art::Run const & r)
{
  // Implementation of optional member function here.
}

void evg::GenieOutput::endSubRun(art::SubRun const & sr)
{
  // Implementation of optional member function here.
}

void evg::GenieOutput::reconfigure(fhicl::ParameterSet const & p)
{
  // Implementation of optional member function here.
}

void evg::GenieOutput::respondToCloseInputFile(art::FileBlock const & fb)
{
  // Implementation of optional member function here.
}

void evg::GenieOutput::respondToCloseOutputFiles(art::FileBlock const & fb)
{
  // Implementation of optional member function here.
}

void evg::GenieOutput::respondToOpenInputFile(art::FileBlock const & fb)
{
  // Implementation of optional member function here.
}

void evg::GenieOutput::respondToOpenOutputFiles(art::FileBlock const & fb)
{
  // Implementation of optional member function here.
}
*/

DEFINE_ART_MODULE(evg::GenieOutput)