CVNAddTrainingData_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    CVNAddTrainingData_module.cc
// \brief   Producer module for adding CVN TrainingData objects
// \author  Grant Nikseresht - gnikseresht@hawk.iit.edu
////////////////////////////////////////////////////////////////////////

// C/C++ includes
#include <iostream>
#include <sstream>
#include <vector>

// ROOT includes
#include "TTree.h"
#include "TH2F.h"

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/FindManyP.h"

// NOvASoft includesOkay
#include "Geometry/Geometry.h"
#include "Calibrator/Calibrator.h"
#include "MCCheater/BackTracker.h"
#include "ShowerLID/ShowerLID.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/FilterList.h"
#include "RecoBase/Prong.h"
#include "RecoBase/Shower.h"
#include "RecoBase/Vertex.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "Utilities/AssociationUtil.h"
#include "canvas/Persistency/Common/Assns.h"

#include "CVN/func/AssignLabels.h"
#include "CVN/func/TrainingData.h"
#include "CVN/func/ProngType.h"
#include "CVN/func/ProngTrainingData.h"
#include "CVN/func/InteractionType.h"
#include "CVN/func/FinalStateType.h"
#include "CVN/func/FinalStateProngsType.h"
#include "CVN/func/ParentParticleType.h"
#include "CVN/func/EventLabeledPMaps.h"

#include "CVN/func/ParticlesTrainingData.h"
#include "CVN/func/ParticlesType.h"

// Required for CAFCutter
#include "RecoValidation/CAFCutter.h"
#include "StandardRecord/StandardRecord.h"

namespace cvn {

  class CVNAddTrainingData : public art::EDProducer {
  public:
    explicit CVNAddTrainingData(fhicl::ParameterSet const& pset);
    ~CVNAddTrainingData();

    void produce(art::Event& evt);
    void reconfigure(const fhicl::ParameterSet& pset);
    void beginJob();
    void endJob();



  private:
    std::string fClusterLabel;
    std::string fPixelMapInput;
    std::string fChooseLabels;
    std::string fProngModLabel;
    std::string fProng3DLabel;
    std::string fProng2DLabel;
    std::string fClusterTDLabel;
    std::string fProngTDLabel;
    std::string fPMapTDLabel;
    std::string fProngPMapTDLabel;
    std::string fVertexLabel;
    std::string fShowerLabel;
    std::string fShowerLIDLabel;
    bool fTDProngs;



    bool interactionTraining = true;
    bool particlesTraining = false;

    TrainingData* fTrain;
    ParticlesTrainingData* fParticlesTrain;

    art::ServiceHandle<cheat::BackTracker> fBT;

    ///  Check rb::IsFiltered?
    bool        fObeyPreselection;
    std::vector<std::string> fPreselectionLabels;

  };

  CVNAddTrainingData::CVNAddTrainingData(fhicl::ParameterSet const& pset)
  {
    this->reconfigure(pset);
    produces< std::vector<cvn::TrainingData> >(fClusterTDLabel);
    produces< std::vector<cvn::ProngTrainingData> >(fProngTDLabel);
    produces< std::vector<cvn::TrainingData> >(fPMapTDLabel);
    produces< std::vector<cvn::ProngTrainingData> >(fProngPMapTDLabel);
    produces< art::Assns<cvn::TrainingData, rb::Cluster> >();
    produces< art::Assns<cvn::ProngTrainingData, rb::Prong> >();
    produces< art::Assns<cvn::TrainingData, cvn::PixelMap> >();
    produces< art::Assns<cvn::ProngTrainingData, cvn::PixelMap> >();
  }
  //.......................................................................
  void CVNAddTrainingData::reconfigure(const fhicl::ParameterSet& pset)
  {
    fClusterLabel     = pset.get<std::string>  ("ClusterLabel");
    fProngModLabel    = pset.get<std::string>  ("ProngModLabel");
    fProng3DLabel     = pset.get<std::string>  ("Prong3DLabel");
    fProng2DLabel     = pset.get<std::string>  ("Prong2DLabel");
    fClusterTDLabel   = pset.get<std::string>  ("ClusterTDLabel");
    fProngTDLabel     = pset.get<std::string>  ("ProngTDLabel");
    fProngPMapTDLabel = pset.get<std::string>  ("ProngPMapTDLabel");
    fPMapTDLabel      = pset.get<std::string>  ("PMapTDLabel");
    fPixelMapInput    = pset.get<std::string>  ("PixelMapInput");
    fShowerLabel      = pset.get<std::string>  ("ShowerLabel");
    fShowerLIDLabel   = pset.get<std::string>  ("ShowerLIDLabel");
    fVertexLabel      = pset.get<std::string>  ("VertexLabel");
    fChooseLabels     = pset.get<std::string>  ("ChooseLabels");
    fObeyPreselection = pset.get<bool>         ("ObeyPreselection");
    fTDProngs         = pset.get<bool>         ("UseTDProngs");
    fPreselectionLabels = pset.get< std::vector<std::string> >("PreselectionLabels");
  }

  //......................................................................
  CVNAddTrainingData::~CVNAddTrainingData()
  {
    //======================================================================
    // Clean up any memory allocated by your module
    //======================================================================
  }

  //......................................................................
  void CVNAddTrainingData::beginJob()
  {
    std::cout << "Beginning job" << '\n';
  }

  //......................................................................
  void CVNAddTrainingData::endJob()
  {
  }

  //......................................................................
  void CVNAddTrainingData::produce(art::Event& evt)
  {
    //std::cout << "Producing an art event in CVNAddTrainingData" << std::endl;

    // This section fills the labels needed for the TrainingData object
    //----- Derived from CVNEventDump_module -------------//
    // get the slices

    art::ServiceHandle<geo::Geometry> geom;

    art::Handle< std::vector< rb::Cluster > > clusterHandle;
    evt.getByLabel(fClusterLabel, clusterHandle);
    const std::vector<rb::Cluster> & clusters = *clusterHandle;
    art::PtrVector<rb::Cluster> slices;
    for(unsigned int i = 0; i < clusterHandle->size(); ++i){
      slices.push_back(art::Ptr<rb::Cluster>(clusterHandle, i));
    }

    std::unique_ptr<std::vector<TrainingData>>
                                  trainData(new std::vector<TrainingData>);
    std::unique_ptr<std::vector<ProngTrainingData>>
                                  prongtrainData(new std::vector<ProngTrainingData>);
    std::unique_ptr<std::vector<TrainingData>>
                                  pmaptrainData(new std::vector<TrainingData>);
    std::unique_ptr<std::vector<ProngTrainingData>>
                                  pmapprongtrainData(new std::vector<ProngTrainingData>);
    std::unique_ptr< art::Assns<TrainingData, rb::Cluster> >
                                  slice_assoc(new art::Assns<TrainingData, rb::Cluster>);
    std::unique_ptr< art::Assns<ProngTrainingData, rb::Prong> >
                                  prong_assoc(new art::Assns<ProngTrainingData, rb::Prong>);
    std::unique_ptr< art::Assns<TrainingData, PixelMap> >
                                  pmap_assoc(new art::Assns<TrainingData, PixelMap>);
    std::unique_ptr< art::Assns<ProngTrainingData, PixelMap> >
                                  prongpmap_assoc(new art::Assns<ProngTrainingData, PixelMap>);
    art::FindManyP<PixelMap> fmPixelMap(clusterHandle, evt, fPixelMapInput);

    std::vector<int> bestNuId;
    std::vector<cheat::NeutrinoWithIndex> nus;
    std::vector< std::vector<cheat::NeutrinoEffPur> > sEffPur;
    //std::cout << "Has truth info: " << fBT->HaveTruthInfo() << std::endl;
    if(fBT->HaveTruthInfo()){
      nus      = fBT->allMCTruth();
      sEffPur  = fBT->SlicesToMCTruthsTable(slices);
      bestNuId = fBT->SliceToOrderedNuIdsByEnergy(nus, sEffPur);
    }

    // for loop over slices
    for(size_t iClust = 0; iClust < clusters.size(); ++iClust){

      if(fObeyPreselection && rb::IsFiltered(evt, clusterHandle, iClust,
                                             fPreselectionLabels)) continue;

      std::vector< art::Ptr<caf::StandardRecord> > records;
      bool pass = true;

      if(!pass) continue;

      const rb::Cluster& cluster = clusters[iClust];
      if(!fmPixelMap.isValid()) continue;

      const std::vector< art::Ptr<PixelMap> > pixelMaps = fmPixelMap.at(iClust);

      if(pixelMaps.empty()) continue;

      InteractionType      interaction      = kOther;
      FinalStateType       finalstate       = kOtherFS;
      FinalStateProngsType finalstateprongs = kOtherFSp;
      ParentParticleType   parentparticle   = kUnknownPP;

      ParticlesType particles_type = kOther_PT;

      float nuEnergy  = 0;
      float lepEnergy = 0;
      double vtxx = std::numeric_limits<double>::lowest();
      double vtxy = std::numeric_limits<double>::lowest();
      double vtxz = std::numeric_limits<double>::lowest();

      art::Ptr<simb::MCTruth> truth;
      if(fBT->HaveTruthInfo()){
        if(bestNuId[iClust] != -1) {
          truth =  sEffPur[iClust][bestNuId[iClust]].neutrinoInt;

          interaction = SliceClassify(truth, cluster, nuEnergy, lepEnergy,
                                  vtxx, vtxy, vtxz);
          particles_type = ParticlesSliceClassify(truth, cluster, nuEnergy,
                                            lepEnergy, vtxx, vtxy, vtxz);
          finalstate       = GetFinalStateType(truth);
          finalstateprongs = GetFinalStateProngsType(truth);
          parentparticle   = GetParentParticleType(truth, evt);
        }
      }
      else particles_type = kCosmic_PT;

      //----- End CVNEventDump_module duplication -------//
      // Now the TrainingData object can be filled with volues
      // set in the previous section

      // std::cout << "Interaction: " << interaction << '\n';
      // std::cout << "ParentParticle: " << parentparticle << '\n';
      // std::cout << "FinalState: " << finalstate << '\n';
      // std::cout << "FinalStateProngs: " << finalstateprongs << '\n';
      // std::cout << "ParticlesType: " << particles_type << '\n';
      // std::cout << "nuEnergy: " << nuEnergy << '\n';
      // std::cout << "lepEnergy: " << lepEnergy << '\n';
      // std::cout << "vtxx: " << vtxx << '\n';
      // std::cout << "vtxy: " << vtxy << '\n';
      // std::cout << "vtxz: " << vtxz << '\n';

      TrainingData train(interaction, parentparticle, finalstate,
                         finalstateprongs, particles_type, nuEnergy,
                         lepEnergy, vtxx, vtxy, vtxz, *pixelMaps[0]);
      trainData->push_back(train);
      pmaptrainData->push_back(train);
      util::CreateAssn(*this, evt, *(trainData.get()), slices[iClust], *(slice_assoc.get()), UINT_MAX, fClusterTDLabel);
      util::CreateAssn(*this, evt, *(pmaptrainData.get()), pixelMaps[0], *(pmap_assoc.get()), UINT_MAX, fPMapTDLabel);
      if(fTDProngs){

        art::FindManyP<rb::Vertex> fmv(clusterHandle, evt, fVertexLabel);
        art::FindManyP<rb::Prong> fmProng3D(clusterHandle, evt,
                              art::InputTag(fProngModLabel, fProng3DLabel));
        art::FindManyP<rb::Prong> fmProng2D(clusterHandle, evt,
                              art::InputTag(fProngModLabel, fProng2DLabel));

        std::vector<art::Ptr<rb::Prong>> prongs3D;
        std::vector<art::Ptr<rb::Prong>> prongs2D;
        std::vector< art::Ptr<rb::Vertex> > vert;

        if( fmProng2D.isValid() ) prongs2D = fmProng2D.at(iClust);
        if( fmProng3D.isValid() ) prongs3D = fmProng3D.at(iClust);
        if(fmv.isValid()       ) vert     = fmv.at(iClust);
        else continue;

        art::FindManyP<PixelMap> fmPixelMap3D(prongs3D, evt, fPixelMapInput);
        art::FindManyP<PixelMap> fmPixelMap2D(prongs2D, evt, fPixelMapInput);
        art::FindOneP<rb::Shower>        foSh(prongs3D, evt, fShowerLabel);

        TVector3 v;
        if(!vert.empty() ) v = vert[0]->GetXYZ();
        else v.SetXYZ(-999.,-999.,-999.);

        double   length;
        TVector3 prongdir;
        TVector3 beamdir  = geom->NuMIBeamDirection();
        float    angle;


        for( unsigned int iProng = 0; iProng < prongs3D.size(); ++iProng ){

          length = -5.0;
          angle  = -5.0;

          if(!fmPixelMap3D.isValid()) continue;

          const std::vector< art::Ptr<PixelMap> > pixelMaps3D = fmPixelMap3D.at(iProng);
          if(pixelMaps3D.empty()) continue;

          double gap = -5.0;
          if(foSh.isValid()){
            cet::maybe_ref< art::Ptr<rb::Shower> const > roSh(foSh.at(iProng));
            art::Ptr<rb::Shower> shower = roSh.ref();

            std::vector< art::Ptr<rb::Shower> > showers;
            showers.push_back(shower);

            art::FindOneP<slid::ShowerLID> foShLID(showers, evt,fShowerLIDLabel);
            if (foShLID.isValid()){
              cet::maybe_ref< art::Ptr<slid::ShowerLID> const > roShLID(foShLID.at(0));
              art::Ptr<slid::ShowerLID> showerLID = roShLID.ref();
              gap = showerLID->Gap();
            }
          }

          ProngTrainingData p_train;
          ProngType ptype3D, ptypeX, ptypeY;
          double purity3D, purityX, purityY, recE;
          unsigned int ncellX, ncellY;
          bool isprimary; 

          ptype3D = ProngClassify(*prongs3D[iProng], &ptype3D, &ptypeX, &ptypeY,
                                  &isprimary, &purity3D, &purityX, &purityY,
                                  &recE, &ncellX, &ncellY);

          p_train = ProngTrainingData(ptype3D, ptypeX, ptypeY, isprimary,
                                      purity3D, purityX, purityY, recE,
                                      ncellX, ncellY, v.X(), v.Y(), v.Z(),
                                      length, (double)angle, gap, *pixelMaps3D[0]);

          prongtrainData->push_back(p_train);
          pmapprongtrainData->push_back(p_train);
          util::CreateAssn(*this, evt, *(prongtrainData.get()),
                           prongs3D[iProng], *(prong_assoc.get()), UINT_MAX, fProngTDLabel); 
          util::CreateAssn(*this, evt, *(pmapprongtrainData.get()),
                           pixelMaps3D[0], *(prongpmap_assoc.get()), UINT_MAX, fProngPMapTDLabel); 

        }
        for( unsigned int iProng2 = 0; iProng2 < prongs2D.size(); ++iProng2 ){

          if(!fmPixelMap2D.isValid()) continue;

          const std::vector< art::Ptr<PixelMap> > pixelMaps2D = fmPixelMap2D.at(iProng2);
          if(pixelMaps2D.empty()) continue;

          ProngTrainingData p2_train;
          ProngType ptype2D, ptypeX2D, ptypeY2D;
          double purity2D, purityX2D, purityY2D, recE;
          unsigned int ncellX, ncellY;
          bool isprimary;

          ptype2D = ProngClassify(*prongs2D[iProng2], &ptype2D, &ptypeX2D, &ptypeY2D,
                                  &isprimary, &purity2D, &purityX2D, &purityY2D,
                                  &recE, &ncellX, &ncellY);

          p2_train = ProngTrainingData(ptype2D, ptypeX2D, ptypeY2D, isprimary,
                                      purity2D, purityX2D, purityY2D, recE,
                                      ncellX, ncellY, v.X(), v.Y(), v.Z(),
                                      0.0, 0.0, 0.0, *pixelMaps2D[0]);

          prongtrainData->push_back(p2_train);
          pmapprongtrainData->push_back(p2_train);
          util::CreateAssn(*this, evt, *(prongtrainData.get()),
                           prongs2D[iProng2], *(prong_assoc.get()), UINT_MAX, fProngTDLabel);
          util::CreateAssn(*this, evt, *(pmapprongtrainData.get()),
                           pixelMaps2D[0], *(prongpmap_assoc.get()), UINT_MAX, fProngPMapTDLabel);

        }
      } //fTDProngs
    } // end for loop over slices
    evt.put(std::move(trainData), fClusterTDLabel);
    evt.put(std::move(prongtrainData), fProngTDLabel);
    evt.put(std::move(pmaptrainData), fPMapTDLabel);
    evt.put(std::move(pmapprongtrainData), fProngPMapTDLabel);
    evt.put(std::move(slice_assoc));
    evt.put(std::move(prong_assoc));
    evt.put(std::move(pmap_assoc));
    evt.put(std::move(prongpmap_assoc));
    //std::cout << "Training data object created." << std::endl;
  //----------------------------------------------------------------------
  }

  DEFINE_ART_MODULE(CVNAddTrainingData)
} // end namespace cvn
////////////////////////////////////////////////////////////////////////