G4Gen_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   A module for propagating particles through Geant4 and
///          generating energy depositions
/// \author  messier@indiana.edu brebel@fnal.gov rhatcher@fnal.gov
////////////////////////////////////////////////////////////////////////


// NOvA includes 
#include "g4nova/G4Alg.h"
#include "g4nova/PrintTables.h"
#include "g4nova/ParticleListAction.h"
#include "Simulation/FLSHitList.h"
#include "Simulation/Particle.h"
#include "Simulation/Simulation.h"
#include "Simulation/TrueEnergy.h"
#include "Simulation/G4Status.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "Utilities/AssociationUtil.h"

#include <cassert>
#include <cstdlib>
#include <string>
#include <sstream>
#include <sys/stat.h>
#include <vector>
#include <map>
#include <unistd.h>

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

// Framework includes
#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 "art/Framework/Core/ModuleMacros.h"

namespace art { 
  class Event; 
  class ParameterSet;
}

namespace g4n {

  /// A module to pass interaction information to Geant4 and 
  /// create hits and particle lists
  class G4Gen : public art::EDProducer {
  public:
    explicit G4Gen(fhicl::ParameterSet const &pset);
    virtual ~G4Gen();                        

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

  private:
    std::string         fGeneratorLabel;

    g4n::G4Alg*         fG4Alg;     ///< G4Helper object
    fhicl::ParameterSet fG4AlgPSet; ///< parameter set to configure the G4Alg object
    TStopwatch          fStopwatch; ///< keep track of how long it takes to run the job

    bool                fPrintTables; ///< Print the dE/dx tables and other info for debugging etc.?

    //bool                fDoFLSHitListDebug; ///< Are we going to debug the FLSHitList?
    
  };
}

namespace g4n {

  //___________________________________________________________________________
  // Constructor
  G4Gen::G4Gen(fhicl::ParameterSet const& pset)
  : fG4Alg(0)
  , fG4AlgPSet(pset.get< fhicl::ParameterSet >("G4AlgPSet") )
  {
    fGeneratorLabel = fG4AlgPSet.get<std::string>("GenModuleLabel");
    fPrintTables    = fG4AlgPSet.get<bool>       ("PrintTables"   );

    // get the random number seed, use a random default if not specified
    // in the configuration file. 
    unsigned int seed = pset.get< unsigned int >("Seed", sim::GetRandomNumberSeed());
    // setup the random number service for Geant4, the "G4Engine" label is a 
    // special tag setting up a global engine for use by Geant4/CLHEP
    createEngine(seed, "G4Engine");

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

    produces< std::vector<sim::G4Status>                 >();
    produces< std::vector<sim::FLSHitList >              >();
    produces< std::vector<sim::Particle>                 >();
    produces< art::Assns<sim::Particle, simb::MCTruth>   >();
    produces< std::vector<sim::TrueEnergy>               >();
    produces< art::Assns<sim::TrueEnergy, sim::Particle> >();

    /// Setup the FLSHitList debugger
    //fDoFLSHitListDebug = doG4Debug<FLSHitListDebug>();
    // if(fDoFLSHitListDebug){
//         mf::LogVerbatim("G4Gen") << "Will be Debugging FLSHitLists";
//     }
//     else{
//         mf::LogVerbatim("G4Gen") << "Will NOT be Debugging FLSHitLists";
//     }


  }// end of constructor

  //___________________________________________________________________________
  // Destructor
  G4Gen::~G4Gen()  
  {  
    fStopwatch.Stop();
    mf::LogVerbatim("G4Gen") << "real time to run through Geant4: " 
                             << fStopwatch.RealTime();
    if(fG4Alg) delete fG4Alg;
  }// end of destructor

  //___________________________________________________________________________
  void G4Gen::beginRun(art::Run& run) 
  {
    fG4Alg = new G4Alg(fG4AlgPSet);
    
    return;
  }// end of begin job

  //___________________________________________________________________________
  void G4Gen::endJob() 
  {
    if(fPrintTables) {
      PrintTables PT;
      PT.Materials();
      PT.dEdx();
    }

    return;
  }// end of endJob

  //___________________________________________________________________________
  void G4Gen::produce(art::Event& evt)  
  {
    // need to instantiate the G4Helper in order to make the 
    // G4RunManager and pass it the physics list - G4 really dislikes
    // it if you try to make a G4UserRunAction (ie UserActionManager)
    // before passing a physics list to the run manager
  
    // Define the FLSHitList and Particle vectors.
    std::unique_ptr<std::vector<sim::G4Status>     >            g4status(new std::vector<sim::G4Status>(1) );
    std::unique_ptr<std::vector<sim::FLSHitList>   >            flshlcol(new std::vector<sim::FLSHitList>  );
    std::unique_ptr<std::vector<sim::Particle>     >            pcol    (new std::vector<sim::Particle>    );
    std::unique_ptr< art::Assns<sim::Particle, simb::MCTruth> > tpassn  (new art::Assns<sim::Particle, simb::MCTruth>);

    std::unique_ptr<std::vector<sim::TrueEnergy> >               trueens (new std::vector<sim::TrueEnergy>);
    std::unique_ptr<art::Assns<sim::TrueEnergy, sim::Particle> > treassn (new art::Assns<sim::TrueEnergy, sim::Particle>);

    //look to see if there is any MCTruth information for this event
    art::Handle<std::vector<simb::MCTruth>> mclist;
    evt.getByLabel(fGeneratorLabel, mclist);

    // the next steps are a little clunky.  Make a vector of art::Ptr<MCTruth> 
    // objects to use when making the associations and another of just
    // const* MCTruth objects to pass to G4Alg
    std::vector< art::Ptr< simb::MCTruth > > mctp;
    std::vector< const simb::MCTruth* > mct;

    for(size_t i = 0; i < mclist->size(); ++i){
      art::Ptr<simb::MCTruth> ptr(mclist, i);
      mctp.push_back(ptr);
      mct .push_back(ptr.get());
    }

    // make a map to keep track of which G4 track ID goes with which MCTruth object
    std::map<int, size_t> trackIDToMCTruthIndex;
    fG4Alg->RunGeant(mct, *flshlcol, *pcol, trackIDToMCTruthIndex);

    // Check to see if Geant4 aborted during processing
    g4status->at(0).SetAborted(fG4Alg->IsAborted());
    if (g4status->at(0).IsAborted()) {
      std::cout << "Geant4 aborted during processing. Flagged in G4Status data product." << std::endl;
    }

    // Karl::Containment -- Let's look at the map which I have just made...
    // --- Make an iterator to loop over fTrueDepMap
    std::map< int, std::vector<float> >::iterator TrueMapIt;
    // --- Iterate through my fTrueDepMap and create assns to the MCParticle list.
    for (TrueMapIt = ParticleListAction::fTrueDepMap.begin(); TrueMapIt != ParticleListAction::fTrueDepMap.end(); ++TrueMapIt) {
      // --- Loop through pcol until I find the right MCParticle.
      for (unsigned int pp=0; pp<pcol->size(); ++pp) {
        // --- Check if I have the correct MCParticle.
        if ( (*pcol)[pp].TrackId() != TrueMapIt->first ) continue;
        // --- If I'm here I should have the right particle so lets get it
        sim::Particle &part = (*pcol)[ pp ];
        // --- Check that I have done this correctly.
        if (part.TrackId() != TrueMapIt->first) {
          mf::LogVerbatim("G4Gen") << "Why doesn't TrueMapIt->first " << TrueMapIt->first << " equal part.TrackId() " << part.TrackId() << std::endl;
          continue;
        }
        // Make a sim::TrueEnergy( TrackID, Entering Energy, Escaping Energy, Total Deposited Energy, Total Escaping energy )
        sim::TrueEnergy ThisEn(TrueMapIt->first, TrueMapIt->second[0], TrueMapIt->second[1], TrueMapIt->second[2], TrueMapIt->second[3]);
        (*trueens).push_back( ThisEn );
        // And make the association!
        util::CreateAssn( *this, evt, *treassn, *trueens, *pcol, pp, (trueens->size()-1) );
        // --- Don't want to carry on looping through all of the particles!
        break;
      } // Loop through particles.
    } // Iterator through my map.
    // Karl::Containment -- Let's look at the map which I have just made...

    // make associations for the particles and MCTruth objects
    int    trackID = INT_MAX;
    size_t mctidx  = 0;

    for(size_t p = 0; p < pcol->size(); ++p){

      sim::Particle &part = (*pcol)[p];
      trackID = part.TrackId();

      if( trackIDToMCTruthIndex.count(trackID) > 0){
        mctidx = trackIDToMCTruthIndex.find(trackID)->second;
        util::CreateAssn(*this, evt, *pcol, mctp[mctidx], *tpassn, p);
      }
      else
      throw cet::exception("G4Gen") << "Cannot find MCTruth for Track Id: "
                                    << trackID
                                    << " to create association between Particle and MCTruth";

    }// end loop over handles

    // Put the data products into the event
    evt.put(std::move(g4status));
    evt.put(std::move(flshlcol));
    evt.put(std::move(pcol));
    evt.put(std::move(tpassn));
    
    evt.put(std::move(trueens));
    evt.put(std::move(treassn));

    return;
  }// end of produce
  
} // end namespace

namespace g4n { DEFINE_ART_MODULE(G4Gen) }