TrackOverlapECalc_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// brief:       Cacluate the overlapping energy on each track using the 
//              TrackCleanUpAlg class, and put that number in the event.
// Module Type: producer
// File:        TrackOverlapECalc_module.cc
//
// Generated on Thu Aug 03 2017, at high noon by Michael Baird
////////////////////////////////////////////////////////////////////////

#include <string>
#include "TMath.h"

// Framework includes
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Principal/Event.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/Assns.h"

// NOvA includes 
#include "Geometry/Geometry.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "RecoBase/Energy.h"
#include "RecoBase/FitSum.h"
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/MathUtil.h"
#include "MuonRemove/art/TrackCleanUpAlg.h"



namespace numue {

  class TrackOverlapECalc : public art::EDProducer {
  public:
    explicit TrackOverlapECalc(fhicl::ParameterSet const & pset);
    virtual ~TrackOverlapECalc();
    
    void produce(art::Event & evt);
    void beginRun(art::Run& run);
    
  private:
    
    murem::TrackCleanUpAlg*  fTrackCleanUpAlg;     ///< Track Clean Up Algorithm Object
    fhicl::ParameterSet      fTrackCleanUpAlgPSet; ///< Parameter Set to configure the TrackCleanUpAlg object
    std::string              fSlicerModuleLabel;   ///< Label for the slicing module
    std::string              fTrackModuleLabel_1;  ///< Label for the first tracking module (and also the instance label for the energy objects)
    std::string              fTrackModuleLabel_2;  ///< Label for the second tracking module (and also the instance label for the energy objects)

    double fMUmipHigh; ///< high mip value for muon tracks
    double fMUmipLow;  ///< low mip value for muon tracks
    double fMUmip;     ///< mip value for muon tracks
    double fMUvertex;  ///< dE/dx cut off value for determining the vertex region for muon tracks

    double fPImipHigh; ///< high mip value for pion tracks
    double fPImipLow;  ///< low mip value for pion tracks
    double fPImip;     ///< mip value for pion tracks
    double fPIvertex;  ///< dE/dx cut off value for determining the vertex region for pion tracks

    double fPRmipHigh; ///< high mip value for proton tracks
    double fPRmipLow;  ///< low mip value for proton tracks
    double fPRmip;     ///< mip value for proton tracks
    double fPRvertex;  ///< dE/dx cut off value for determining the vertex region for proton tracks

  };
  
  //----------------------------------------------------------------------
  TrackOverlapECalc::TrackOverlapECalc(fhicl::ParameterSet const & pset):
    fTrackCleanUpAlg     ( 0 ),
    fTrackCleanUpAlgPSet ( pset.get< fhicl::ParameterSet > ("TrackCleanUpAlgPSet") ),
    fSlicerModuleLabel   ( pset.get< std::string         > ("SlicerModuleLabel"  ) ),
    fTrackModuleLabel_1  ( pset.get< std::string         > ("TrackModuleLabel_1" ) ),
    fTrackModuleLabel_2  ( pset.get< std::string         > ("TrackModuleLabel_2" ) ),
    fMUmipHigh           ( pset.get< double > ("MUmipHigh"            ) ),
    fMUmipLow            ( pset.get< double > ("MUmipLow"             ) ),
    fMUmip               ( pset.get< double > ("MUmip"                ) ),
    fMUvertex            ( pset.get< double > ("MUvertexRegionCutOff" ) ),
    fPImipHigh           ( pset.get< double > ("PImipHigh"            ) ),
    fPImipLow            ( pset.get< double > ("PImipLow"             ) ),
    fPImip               ( pset.get< double > ("PImip"                ) ),
    fPIvertex            ( pset.get< double > ("PIvertexRegionCutOff" ) ),
    fPRmipHigh           ( pset.get< double > ("PRmipHigh"            ) ),
    fPRmipLow            ( pset.get< double > ("PRmipLow"             ) ),
    fPRmip               ( pset.get< double > ("PRmip"                ) ),
    fPRvertex            ( pset.get< double > ("PRvertexRegionCutOff" ) )
  { 
    this->produces< std::vector<rb::Energy> >           (fTrackModuleLabel_1); // using track label as instance label too...
    this->produces< art::Assns<rb::Energy, rb::Track> > (fTrackModuleLabel_1); // using track label as instance label too...
    this->produces< std::vector<rb::Energy> >           (fTrackModuleLabel_2); // using track label as instance label too...
    this->produces< art::Assns<rb::Energy, rb::Track> > (fTrackModuleLabel_2); // using track label as instance label too...
  }
  
  
  //----------------------------------------------------------------------
  TrackOverlapECalc::~TrackOverlapECalc()
  {
    if (fTrackCleanUpAlg) delete fTrackCleanUpAlg;
  }
  
  //----------------------------------------------------------------------
  void TrackOverlapECalc::beginRun(art::Run& run)
  {
    // Have to recreate these because behavior can depend on run number
    // (apparently...)
    if(fTrackCleanUpAlg) delete fTrackCleanUpAlg;
    fTrackCleanUpAlg = new murem::TrackCleanUpAlg(fTrackCleanUpAlgPSet);
  }

  //----------------------------------------------------------------------
  void TrackOverlapECalc::produce(art::Event& evt)
  {
    
    // Load slices list from the event.
    art::Handle<std::vector<rb::Cluster> > slicecol;
    evt.getByLabel(fSlicerModuleLabel, slicecol);
        
    // Get the first set of tracks associated with the slices.
    art::FindManyP<rb::Track> sliceToTracks_1(slicecol, evt, fTrackModuleLabel_1);

    // Get the second set of tracks associated with the slices.
    art::FindManyP<rb::Track> sliceToTracks_2(slicecol, evt, fTrackModuleLabel_2);

    // Declare containers for things to be stored in event.
    std::unique_ptr< std::vector<rb::Energy> >           trackECol_1 (new std::vector<rb::Energy>);
    std::unique_ptr< art::Assns<rb::Energy, rb::Track> > trackAssoc_1(new art::Assns<rb::Energy, rb::Track>);
    std::unique_ptr< std::vector<rb::Energy> >           trackECol_2 (new std::vector<rb::Energy>);
    std::unique_ptr< art::Assns<rb::Energy, rb::Track> > trackAssoc_2(new art::Assns<rb::Energy, rb::Track>);



    // Loop over the slices
    for(size_t islice = 0; islice < slicecol->size(); ++islice){
      const art::Ptr<rb::Cluster> slice(slicecol, islice);

      // As usual, skip the noise slice.
      if(slice->IsNoise()) continue;

      // Get all of the tracks from the first tracker in the slice.
      const std::vector< art::Ptr<rb::Track> > sliceTracks_1 = sliceToTracks_1.at(islice);

      // Get all of the tracks from the second tracker in the slice.
      const std::vector< art::Ptr<rb::Track> > sliceTracks_2 = sliceToTracks_2.at(islice);

      // Create the FindMany object for getting FitSum objects created by the second tracker.
      art::FindManyP<rb::FitSum> fitsumFMP(sliceTracks_2, evt, fTrackModuleLabel_2);



      // Loop over the first set of tracks:
      //
      // NOTE: These tracks are assumed to be the Kalman tracks, therefore
      //       we will use the muon assumption for all tracks.
      if( !(sliceTracks_1.empty()) ) {
        for(unsigned int itrack = 0; itrack < sliceTracks_1.size(); ++itrack) {
          
          // make the energy object
          rb::Energy trackE;
          
          // call TrackCleanUpAlg and set the energy
          trackE.SetE( (float)(fTrackCleanUpAlg->ExtraEOnTrackInGeV(*sliceTracks_1[itrack], *slice,
                                                                    &fMUmipHigh, &fMUmipLow, &fMUmip, &fMUvertex)));
          trackECol_1->push_back(trackE);

          // make the association
          util::CreateAssn(*this, evt, *(trackECol_1.get()), sliceTracks_1[itrack], *(trackAssoc_1.get()),
                           UINT_MAX, fTrackModuleLabel_1); // using track label as instance label too...
          
        } // End loop over first set of tracks (itrack)
      } // End if( !(sliceTracks_1.empty()) )



      // Loop over the second set of tracks:
      //
      // NOTE: These tracks are assumed to be the BPF tracks, therefore
      //       we will use the appropriate assumption for each track.
      if( !(sliceTracks_2.empty()) ) {
        for(unsigned int itrack = 0; itrack < sliceTracks_2.size(); ++itrack) {

          // Get the FitSum object associated with this track.
          std::vector< art::Ptr< rb::FitSum > > fitsums;
          if(fitsumFMP.isValid()) fitsums = fitsumFMP.at(itrack);
          if(fitsums.size() == 0) continue; // this should never happen, but just in case...
          
          // make the energy object
          rb::Energy trackE;
          
          // Determine the particle assumption and set the mip values:
          double mipHigh = 0.0;
          double mipLow  = 0.0;
          double mip     = 0.0;
          double vertex  = 0.0;
          if(fitsums[0]->PDG() == 13) { // muon...
            mipHigh = fMUmipHigh;
            mipLow  = fMUmipLow;
            mip     = fMUmip;
            vertex  = fMUvertex;
          }
          else if(fitsums[0]->PDG() == 211) { // pion...
            mipHigh = fPImipHigh;
            mipLow  = fPImipLow;
            mip     = fPImip;
            vertex  = fPIvertex;
          }
          else if(fitsums[0]->PDG() == 2212) { // proton...
            mipHigh = fPRmipHigh;
            mipLow  = fPRmipLow;
            mip     = fPRmip;
            vertex  = fPRvertex;
          }
          else {
            std::cout << "\n\n\nWARNING! PDG code: "
                      << fitsums[0]->PDG()
                      << " not currently supported!\n\n\n";
            continue;
          }

          // call TrackCleanUpAlg and set the energy
          trackE.SetE( (float)(fTrackCleanUpAlg->ExtraEOnTrackInGeV(*sliceTracks_2[itrack], *slice,
                                                                    &mipHigh, &mipLow, &mip, &vertex)));
          trackECol_2->push_back(trackE);

          // make the association
          util::CreateAssn(*this, evt, *(trackECol_2.get()), sliceTracks_2[itrack], *(trackAssoc_2.get()),
                           UINT_MAX, fTrackModuleLabel_2); // using track label as instance label too...
          
        } // End loop over second set of tracks (itrack)
      } // End if( !(sliceTracks_2.empty()) )
      
    } // End of loop over slices (islice).



    // Put objects and associations into the file.
    evt.put(std::move(trackECol_1),  fTrackModuleLabel_1); // using track label as instance label too...
    evt.put(std::move(trackAssoc_1), fTrackModuleLabel_1); // using track label as instance label too...
    evt.put(std::move(trackECol_2),  fTrackModuleLabel_2); // using track label as instance label too...
    evt.put(std::move(trackAssoc_2), fTrackModuleLabel_2); // using track label as instance label too...

    return;

  } // End of produce function.

  //----------------------------------------------------------------------
  DEFINE_ART_MODULE(TrackOverlapECalc)

} // End of namespace