MergeGenCollections_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////
// \file    MergeGenCollections.cxx
// \brief
// \authors janow@fnal.gov
///////////////////////////////////////////////////////////////////////////  

#include <math.h>
#include <string>

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "dk2nu/tree/dk2nu.h"
#include "dk2nu/tree/NuChoice.h"

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

#include "Simulation/FLSHit.h"
#include "Simulation/FLSHitList.h"
#include "Simulation/Particle.h"
#include "SummaryData/POTSum.h"
#include "SummaryData/SpillData.h"
#include "Utilities/AssociationUtil.h"


/// Merge generation collections MCTruth/MCFlux/GTruth/Dk2Nu/NuChoice
namespace MergeGenCollections
{  
  class MergeGenCollections : public art::EDProducer {
  public:

    explicit MergeGenCollections(fhicl::ParameterSet const &pset);
    virtual ~MergeGenCollections();
    
    void produce(art::Event& evt);
    void endSubRun(art::SubRun& sr);

  protected:

    // configuration settings
    std::string  fFirstCollection;    /// MCTruth label for main collection
    std::string  fSecondCollection;   /// Process label for the second collection
    bool fFirstColIsMCNu;             
    bool fSecondColIsMCNu;            
    bool fMergeDk2Nu;                 /// Merge Dk2Nu/NuChoice objects? 
    int  fSecondaryModeOffset;        /// Flag a MCNeutrino as secondary
  };


  //----------------------------------------------------------------------
  MergeGenCollections::MergeGenCollections(fhicl::ParameterSet const &pset) :
    fFirstCollection         (pset.get< std::string >("FirstGenCollection","")),
    fSecondCollection        (pset.get< std::string >("SecondCollection","")),
    fFirstColIsMCNu          (!pset.get< bool       >("IsFirstCollectionCosmic")),
    fSecondColIsMCNu         (!pset.get< bool       >("IsSecondCollectionCosmic")),
    fMergeDk2Nu              (pset.get< bool        >("MergeDk2Nu")),
    fSecondaryModeOffset     (pset.get< int         >("SecondaryModeOffset",0))
  {
    if(fFirstCollection=="" && fSecondCollection==""){
      std::cerr << "MergeGenCollections: "
                << "Neither collection is specified to merge from. Abort." 
                << std::endl;
      std::abort();
    }
    if(!fFirstColIsMCNu && !fSecondColIsMCNu){
      std::cerr << "MergeGenCollections: "
                << "No MC Collections to merge. Abort." 
                << std::endl;
      std::abort();
    }
    
    produces< sumdata::POTSum, art::InSubRun               >();
    produces< sumdata::SpillData                           >();
    produces< std::vector<simb::MCTruth>                   >();
    produces< std::vector<simb::MCFlux>                    >();
    produces< std::vector<simb::GTruth>                    >();
    produces< std::vector<bsim::Dk2Nu>                     >();
    produces< std::vector<bsim::NuChoice>                  >();
    produces< art::Assns<simb::MCTruth, simb::MCFlux>      >();
    produces< art::Assns<simb::MCTruth, simb::GTruth>      >();  
  }
  
  //----------------------------------------------------------------------
  MergeGenCollections::~MergeGenCollections()
  {
  }
  
  //----------------------------------------------------------------------
  void MergeGenCollections::produce(art::Event& evt)
  {
    std::unique_ptr< sumdata::SpillData          > spillcol (new sumdata::SpillData);
    std::unique_ptr< std::vector<simb::MCTruth>  > truthcol (new std::vector<simb::MCTruth >);
    std::unique_ptr< std::vector<simb::MCFlux>   > fluxcol  (new std::vector<simb::MCFlux  >);
    std::unique_ptr< std::vector<simb::GTruth>   > gtruthcol(new std::vector<simb::GTruth  >);
    std::unique_ptr< std::vector<bsim::Dk2Nu>    > dk2nucol (new std::vector<bsim::Dk2Nu   >);
    std::unique_ptr< std::vector<bsim::NuChoice> > nuccol   (new std::vector<bsim::NuChoice>);
    
    std::unique_ptr< art::Assns<simb::MCTruth, simb::MCFlux> > assns (new art::Assns<simb::MCTruth, simb::MCFlux>);
    std::unique_ptr< art::Assns<simb::MCTruth, simb::GTruth> > gassns(new art::Assns<simb::MCTruth, simb::GTruth>);



    //// First Collection
    //
    if(fFirstCollection!=""){ // Dont look for MC if not there

      // Pull SpillData from the first collection only
      art::Handle<sumdata::SpillData> spill;
      evt.getByLabel(fFirstCollection, spill);
      if(spill.isValid()) *spillcol = *spill;
      
      // MCTruth from generator -- not necessarily neutrino MC
      art::Handle< std::vector<simb::MCTruth> > mclist;
      evt.getByLabel(fFirstCollection, mclist);
      
      art::Handle< std::vector<simb::MCFlux>   > fllist;
      art::Handle< std::vector<simb::GTruth>   > glist;
      art::Handle< std::vector<bsim::Dk2Nu>    > dklist;
      art::Handle< std::vector<bsim::NuChoice> > nulist;    
      if(fFirstColIsMCNu) {
        evt.getByLabel(fFirstCollection, fllist);
        evt.getByLabel(fFirstCollection, glist);
        if(fMergeDk2Nu){
          evt.getByLabel(fFirstCollection, dklist);
          evt.getByLabel(fFirstCollection, nulist);
        }
      }
      
      // Loop over MCTruths and associate to Flux/GTruth
      for(unsigned int i = 0; i < mclist->size(); ++i){ 
        simb::MCTruth truth ((*mclist)[i]);
        truthcol->push_back(truth);
        
        if(fFirstColIsMCNu){
          if(i < fllist->size() && i < glist->size()){
            simb::MCFlux  flux  ((*fllist)[i]);
            simb::GTruth gtruth ((*glist)[i]);
            fluxcol  ->push_back(flux);   
            gtruthcol->push_back(gtruth);
            
            util::CreateAssn(*this, evt, *truthcol, *fluxcol,   *assns,  fluxcol  ->size()-1, fluxcol  ->size());
            util::CreateAssn(*this, evt, *truthcol, *gtruthcol, *gassns, gtruthcol->size()-1, gtruthcol->size());
          }
          if(i < dklist->size() && i < nulist->size()){
            bsim::Dk2Nu dk2nu   ((*dklist)[i]);
            dk2nucol->push_back(dk2nu);
            bsim::NuChoice nuc  ((*nulist)[i]);
            nuccol  ->push_back(nuc);
          }
        } // end fFirstColIsMCNu
      } // end loop over MCTruths
    } // If FirstCollection specified



    //// Second Collection
    //
    if(fSecondCollection!=""){ // dont use label if not there
      art::Handle< std::vector<simb::MCTruth> > mc2list;
      evt.getByLabel(fSecondCollection, mc2list);
      
      art::Handle< std::vector<simb::GTruth>   > g2list;
      art::Handle< std::vector<simb::MCFlux>   > fl2list;
      art::Handle< std::vector<bsim::Dk2Nu>    > dk2list;
      art::Handle< std::vector<bsim::NuChoice> > nu2list;
      if(fSecondColIsMCNu) {
        evt.getByLabel(fSecondCollection, g2list);
        evt.getByLabel(fSecondCollection, fl2list);
        evt.getByLabel(fSecondCollection, dk2list);
        evt.getByLabel(fSecondCollection, nu2list);
      }

      for(unsigned int i = 0; i < mc2list->size(); ++i){
        simb::MCTruth tru2 ((*mc2list)[i]);
        
        if(fSecondaryModeOffset){
          // Make new truth so MCNeutrino mode can be
          // offset and used to flag overlaid MCTruth
          simb::MCNeutrino nu2 = tru2.GetNeutrino();
          simb::MCTruth truth;
          truth.SetOrigin(tru2.Origin());
          for( int p=0; p<tru2.NParticles(); p++ ){
            simb::MCParticle mcp2 = tru2.GetParticle(p);
            truth.Add(mcp2);
          }
          // set nu must be after MCParticle list is made
          truth.SetNeutrino( nu2.CCNC(),
                             nu2.Mode() + fSecondaryModeOffset,
                             nu2.InteractionType(),
                             nu2.Target(),    nu2.HitNuc(),
                             nu2.HitQuark(),  nu2.W(), 
                             nu2.X(), nu2.Y(),nu2.QSqr() );
          truthcol->push_back(truth);
          
        } else {
          // don't bother if no offset
          truthcol->push_back(tru2);
        }
        
                
        if(fSecondColIsMCNu) {
          if(i < fl2list->size() && i < g2list->size()){
            simb::MCFlux  flux  ((*fl2list)[i]);
            simb::GTruth gtruth ((*g2list) [i]);
            gtruthcol->push_back(gtruth);
            fluxcol->push_back(flux);

            util::CreateAssn(*this, evt, *truthcol, *fluxcol,   *assns,  fluxcol  ->size()-1, fluxcol  ->size());
            util::CreateAssn(*this, evt, *truthcol, *gtruthcol, *gassns, gtruthcol->size()-1, gtruthcol->size()); 
          }
          if(i < dk2list->size() && i < nu2list->size()){
            bsim::Dk2Nu dk2nu   ((*dk2list)[i]);
            dk2nucol->push_back(dk2nu);
            bsim::NuChoice nuc  ((*nu2list)[i]);
            nuccol  ->push_back(nuc);
          }
        } // end fSecondColIsMCNu
      } // end loop over MCTruths
    } // If SecondCollection specified

    // put the collections in the event
    evt.put(std::move(spillcol));
    evt.put(std::move(truthcol));
    evt.put(std::move(fluxcol));
    evt.put(std::move(gtruthcol));
    evt.put(std::move(dk2nucol));
    evt.put(std::move(nuccol));
    evt.put(std::move(assns));
    evt.put(std::move(gassns));
  }

  //--------------------------------------------------------------------
  void MergeGenCollections::endSubRun(art::SubRun &sr)
  {
    std::unique_ptr< sumdata::POTSum > potcol(new sumdata::POTSum);

    if(fFirstColIsMCNu) {
      art::Handle< sumdata::POTSum > pots;
      sr.getByLabel(fFirstCollection, pots);

      if(pots.isValid()) *potcol = *pots;      
    }
    else if(fSecondColIsMCNu) {
      mf::LogWarning("MergeGenCollections") 
        << "Ignoring second-collection POTSum, using only first";
    }
    sr.put(std::move(potcol));
  }
  
  DEFINE_ART_MODULE(MergeGenCollections)
} // end namespace merge
/////////////////////////////////////////////////////////////////////////