RecVarPID_module.cc

Go to the documentation of this file.

// NOvA includes
#include "RecVarPID/RVP.h"
#include "Calibrator/Calibrator.h"
#include "RecoBase/Track.h"
#include "RecoBase/Prong.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/FilterList.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Vertex.h"
#include "Utilities/func/EnvExpand.h"
#include "Utilities/AssociationUtil.h"

// Framework includes
#include "art/Framework/Core/ModuleMacros.h"
#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 "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "canvas/Persistency/Common/Assns.h"
#include "canvas/Persistency/Common/FindManyP.h"

#include <vector>
#include <string>

// ROOT includes
#include "TObject.h"
#include "TNtuple.h"
#include "TMVA/Factory.h"
#include "TMVA/Tools.h"
#include "TMVA/Reader.h"
#include "TString.h"
#include "TMath.h"


#include <boost/algorithm/string.hpp>
#include <boost/foreach.hpp>
#include <boost/tokenizer.hpp>

/// for std::out
#define MAX_BUFFER         100000

/// In case this will be changed to whatever in future
//#define FUZZYKVERTEXOUTPUT rb::Track
//#define FUZZYKVERTEXOUTPUT rb::Shower
#define FUZZYKVERTEXOUTPUT rb::Prong


namespace rvp {
  class RecVarPID : public art::EDProducer {
  public:
    
    explicit RecVarPID(fhicl::ParameterSet const& pset);
    virtual ~RecVarPID();
    void beginJob();
    void produce(art::Event& evt);
    
  private:

    rb::RecoHit MakeRecoHit(art::Ptr<rb::CellHit> const& chit, 
                            double const& w);
    
    /// Book rvp reader
    bool                     bookVariables(const std::string& weight_file, TMVA::Reader& rvp, std::string& method_name);
    
    
    /// Method to get Variable names from XML
    std::vector<std::string> getVarNamesFromXML(const std::string weights_file_name) const;
    
    
    /// Get TMVA method name  from XML
    std::string              getMethodNameFromXML(const std::string weights_file_name, const bool get_full_name=true) const;
    
  protected:
    std::string                fRVPWeightFile;         ///< RVP weight file name
    std::string                fRVP12WeightFile;       ///< RVP12 weight file name
    //    std::vector<std::string>   fExtendedRVPWeightFile; ///< Extended RVP weight file name
    std::string                fSliceLabel;            ///< Label for slices
    std::string                fTrackLabel;            ///< Label for tracks
    std::string                fProngLabel;            ///< Label for prong
    std::string                fInstLabel3D;           ///< Instance label for 3d prongs
    std::string                fInstLabel2D;           ///< Instance label for 2d prongs
    std::string                fVertexLabel;           ///< Label for vertex
    
    std::vector<Float_t>       fInputVarsRVP;          ///< Vector of the variables that are input for TMVA
    
    TMVA::Reader*              fbdtRVP;                ///< Reader for RVP15
    TMVA::Reader*              fbdtRVP12;              ///< Reader for RVP15
    std::string                fbdtRVPMethodName;
    std::string                fbdtRVP12MethodName;

    bool                       fRunTMVA;               ///< Run the actual (memory-hungry) TMVA step?
  };

  //----------------------------------------------------------------------------
  RecVarPID::RecVarPID(fhicl::ParameterSet const& pset):
    fbdtRVP         (new TMVA::Reader()),
    fbdtRVP12       (new TMVA::Reader())
  {
    /// Reserve the number of variables to be something big
    fInputVarsRVP.resize(10000);

    const std::string bdt  = pset.get< std::string >("BDTweight");
    fRVPWeightFile         = pset.get< std::string >("RVPWeightFile");
    fRVP12WeightFile       = pset.get< std::string >("RVP12WeightFile");
    fSliceLabel            = pset.get< std::string >("SliceLabel");
    fTrackLabel            = pset.get< std::string >("TrackLabel");
    fProngLabel            = pset.get< std::string >("ProngLabel");
    fInstLabel3D             = pset.get< std::string >("InstLabel3D");
    fInstLabel2D             = pset.get< std::string >("InstLabel2D");
    fVertexLabel           = pset.get< std::string >("VertexLabel");

    fRunTMVA               = pset.get< bool >("RunTMVA");

    /// Expand fcl input path to BDTs
    std::string bdtpath = util::EnvExpansion(bdt);
    
    /// Add BDT path to the files
    fRVPWeightFile   = bdtpath + fRVPWeightFile;
    fRVP12WeightFile = bdtpath + fRVP12WeightFile;

    
    produces< std::vector<rvp::RVP> >();
    produces< art::Assns<rvp::RVP,rb::Cluster> >();
    
  }
  
  //...........................................................................
  RecVarPID::~RecVarPID()
  {
    if(fbdtRVP)        delete fbdtRVP;
    if(fbdtRVP12)      delete fbdtRVP12;
    
  }
  
  //...........................................................................
  void RecVarPID::beginJob()
  {
    if(fRunTMVA){
      /// Need to set Weight file
      fbdtRVPMethodName   = "RVP ";
      bookVariables(fRVPWeightFile  , *fbdtRVP  , fbdtRVPMethodName);

      /// Need to set Weight file
      fbdtRVP12MethodName = "RVP12 ";
      bookVariables(fRVP12WeightFile, *fbdtRVP12, fbdtRVP12MethodName);
    }
  }

  //......................................................................
  rb::RecoHit RecVarPID::MakeRecoHit(art::Ptr<rb::CellHit> const& chit, 
                                     double const& w)
  {
    art::ServiceHandle<calib::Calibrator> cal;
    rb::RecoHit rhit(cal->MakeRecoHit(*chit, w));
    
    return rhit;
  }
  
  //...........................................................................
  void RecVarPID::produce(art::Event&  evt)
  {
    art::ServiceHandle<calib::Calibrator> cal;
    
    typedef std::unique_ptr<std::vector<rvp::RVP> >            RVPs_t;
    typedef std::unique_ptr<art::Assns<rvp::RVP,rb::Cluster> > RVPAssn_t;
    RVPs_t                 pidcol (new std::vector<rvp::RVP>);
    RVPAssn_t              assncol(new art::Assns<rvp::RVP,rb::Cluster>);
    
    double BDT12 = -1.0;//set the default 12 variable BDT output to -1
    double BDT15 = -1.0;//set the default 15 variable BDT output to -1
    
    art::Handle< std::vector<rb::Cluster> > slices;
    evt.getByLabel(fSliceLabel, slices);

    
    const art::FindManyP<rb::Track>  fmTrack  (slices, evt, fTrackLabel);
    const art::FindManyP<rb::Vertex> fmVertex (slices, evt, fVertexLabel);
    
    //loop over all slices
    for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
      
      if(rb::IsFiltered(evt, slices, sliceIdx)) continue;
      const rb::Cluster& slice = (*slices)[sliceIdx];
      
      if(slice.IsNoise()) continue;
      
      const unsigned int nCells       = slice.NCell();//# of cells
      const int          lengthPlanes = slice.ExtentPlane();//# of extended planes
      
      //select the longest track
      double longestTrack        = 0.0;//length of longest track
      int    longestTrack_ncells = 0;  //# of cells associated to the longest track

      if (fmTrack.isValid()){
        std::vector<art::Ptr<rb::Track>> tracks = fmTrack.at(sliceIdx);      
        
        for(unsigned int n = 0; n < tracks.size(); ++n){
          const art::Ptr<rb::Track> current_track = tracks[n];
          const double              L = current_track->TotalLength();
          if(L > longestTrack){
            longestTrack        = L;
            longestTrack_ncells = current_track->NCell();
          }
        }
      }
      
      //mip cells and reco energy for slicer
      double recoEnergy = 0.0;//reco energy of slicer
      double nmip       = 0.0;//# of mip cells
      double Ntotal     = 0.0;//# of total cells
      
      const unsigned int ncell_hits = slice.NCell();
      double             wx         = (slice.NXCell() > 0) ? slice.W(slice.XCell(0).get()) : 0;
      double             wy         = (slice.NYCell() > 0) ? slice.W(slice.YCell(0).get()) : 0;
      std::vector<rb::RecoHit> all_reco_hits;
      
      for( unsigned int icell=0; icell<ncell_hits; ++icell){
        const art::Ptr<rb::CellHit>& chit = slice.Cell(icell);
        all_reco_hits.push_back(this->MakeRecoHit(chit, chit->View() == geo::kX ? wx : wy));
        const rb::RecoHit& rhit = all_reco_hits[icell];
        
        if( !rhit.IsCalibrated() ) continue;
        recoEnergy += rhit.GeV();
        
        if( rhit.PECorr()>100. && rhit.PECorr()<245. ) ++nmip;
        ++Ntotal;
      }

      //reco vertex
      //only continue if this get by label succeeded
      if (fmVertex.isValid()){
        std::vector<art::Ptr<rb::Vertex>> vert = fmVertex.at(sliceIdx);
        
        //pass over slice if there was no vertex
        //NOTE: right now it is impossible for a slice to have more then one vertex, so the size
        //should be exactly one
        if (vert.size() != 1) continue;
        //this is extra preselection that is being removed so nothing is hidden from caf
        //if( recoEnergy>0.05 && nCells>5 && lengthPlanes<140 && vert.size()>0 ){
        //if( recoEnergy>0.05 && nCells>5 && lengthPlanes<140 ){
        //calculate per plane cell energy
        std::vector<double> PlaneSumEnergy;
        double totalE_road2sigma = 0.0; //total energy in 2 sigma road
        double Esig_3            = 0.0; //energy not in a 3 sigma window of shower core
        
        for( int iplane=0; iplane<lengthPlanes; ++iplane ){
          double cellX  = 0.0;//energy-weighted cell X in each plane
          double cellY  = 0.0;//energy-weighted cell Y in each plane
          double Ecell  = 0.0;//sum cell energy in each plane
          double EXcell = 0.0;
          double EYcell = 0.0;
          
          std::vector<double> Xxview;
          std::vector<double> Exview;
          std::vector<double> Yyview;
          std::vector<double> Eyview;
          Xxview.clear();
          Exview.clear();
          Yyview.clear();
          Eyview.clear();

          for( unsigned int icell=0; icell<ncell_hits; ++icell){
            const art::Ptr<rb::CellHit>& chit = slice.Cell(icell);
            const rb::RecoHit& rhit = all_reco_hits[icell];
            if( !rhit.IsCalibrated() ) continue;
            
            const int Dplane=chit->Plane()-slice.MinPlane();
            if( Dplane<0 ) std::cout<<"plane # of this cell is beyond the minimal plane."<<std::endl;
            
            if( Dplane==iplane ){
              Ecell += rhit.GeV();
              EXcell += rhit.GeV()*rhit.X();
              EYcell += rhit.GeV()*rhit.Y();
              
              if( chit->View() == geo::kX ){
                 Xxview.push_back(rhit.X());
                 Exview.push_back(rhit.GeV());
              }//in X view
              else {
                 Yyview.push_back(rhit.Y());
                 Eyview.push_back(rhit.GeV());
              }//in Y view
              
            }
          }//loop all cells to calculate the centroid cell position in each plane
          
          if( Ecell>0. ){
            cellX = EXcell/Ecell;
            cellY = EYcell/Ecell;
          }
          
          PlaneSumEnergy.push_back(Ecell);
          
          const int Nxviews=Xxview.size();
          const int Nyviews=Yyview.size();
          if( Nxviews>0 ){
            for( int ixv=0; ixv<Nxviews; ++ixv ){
              if( fabs(Xxview[ixv]-cellX) < 2.*2. ) totalE_road2sigma+=Exview[ixv];
              if( fabs(Xxview[ixv]-cellX) < 2.*3. ) Esig_3+=Exview[ixv];
            }
          }//in X view
          if( Nyviews>0 ){
            for( int iyv=0; iyv<Nyviews; ++iyv ){
              if( fabs(Yyview[iyv]-cellY) < 2.*2. ) totalE_road2sigma+=Eyview[iyv];
              if( fabs(Yyview[iyv]-cellY) < 2.*3. ) Esig_3+=Eyview[iyv];
            }
          }//in Y view

        }//loop all planes
        
        double efrac_2slides = -1.0;//energy fraction in 2 sliding window
        double efrac_6slides = -1.0;//energy fraction in 6 sliding window
        double Eplane20      =  0.0;
        int nplanes=PlaneSumEnergy.size();
        if (recoEnergy > 0.0){
          for( int i=0; i<nplanes; ++i ){
            if( i<20 ) Eplane20 += PlaneSumEnergy[i];
            
            //maximal energy fraction in 2 sliding window
            if( i<nplanes-2 ){
              const double e_frac2=(PlaneSumEnergy[i]+PlaneSumEnergy[i+1])/recoEnergy;
              if( efrac_2slides<e_frac2 ) efrac_2slides=e_frac2;        }
            else if( nplanes<=2 ) efrac_2slides=1.;
            //maximal energy fraction in 6 sliding window
            if( i<nplanes-6 ){
              const double e_frac6=(PlaneSumEnergy[i]+PlaneSumEnergy[i+1]+PlaneSumEnergy[i+2]+PlaneSumEnergy[i+3]+PlaneSumEnergy[i+4]+PlaneSumEnergy[i+5])/recoEnergy;
              if( efrac_6slides<e_frac6 ) efrac_6slides=e_frac6; }
            else if( nplanes<=6 ) efrac_6slides=1.;
          }//loop all planes
        }//end if
        //reco prongs
        //Pulling Tracks from event since this is what FuzzyK makes for others
        //Tracks are then converted to prongs for the rest of this code
        //art::Handle< std::vector<rb::Track> > pngcol;
        //evt.getByLabel(fProngLabel, fInstLabel, pngcol);
        //std::vector<const rb::Prong*> SelectedProng;
        art::FindManyP<FUZZYKVERTEXOUTPUT>        fmProng3D(vert, evt, art::InputTag(fProngLabel, fInstLabel3D));
        art::FindManyP<FUZZYKVERTEXOUTPUT>        fmProng2D(vert, evt, art::InputTag(fProngLabel, fInstLabel2D));
        std::vector<art::Ptr<FUZZYKVERTEXOUTPUT>> prongs = fmProng3D.at(0);
        std::vector<art::Ptr<FUZZYKVERTEXOUTPUT>> SelectedProng;
        std::vector<art::Ptr<FUZZYKVERTEXOUTPUT>> SelectedProng2D;
        //check if instance label is defined to handle sorting of prongs
        if ((fInstLabel3D != "")&&(fInstLabel2D != "")){
          SelectedProng = fmProng3D.at(0);
          SelectedProng2D = fmProng2D.at(0);
        }
        else{
          for( unsigned int ip=0; ip<prongs.size(); ++ip ){
            art::Ptr<FUZZYKVERTEXOUTPUT> current_prong = prongs[ip];
            if(current_prong->IsNoise() )           continue;
            if(current_prong->View() == geo::kXorY)//3D prongs
              SelectedProng.push_back(current_prong);
            else
              SelectedProng2D.push_back(current_prong);
          }
        }
        
        int    nprongs  = SelectedProng.size();
        double Ebalance = 1.0;//energy balance between 2 most energetic prongs
        std::vector<double> ProngEnergy;
        std::vector<double> ProngEnergyXView;
        std::vector<double> ProngEnergyYView;
        for(int ip=0; ip<nprongs; ++ip ){
          double Eprong  = 0.0;
          double EprongX = 0.0;
          double EprongY = 0.0;
          const art::Ptr<FUZZYKVERTEXOUTPUT> current_prong = SelectedProng[ip];
          
          const unsigned int nselectedprong_cells = current_prong->NCell();

          for( unsigned int icell=0; icell<nselectedprong_cells; ++icell){
            const art::Ptr<rb::CellHit>& chit = current_prong->Cell(icell);
            const rb::RecoHit rhit = current_prong->RecoHit(icell);

            if( !rhit.IsCalibrated() ) continue;
            Eprong += rhit.GeV();
            if( chit->View() == geo::kX ) EprongX += rhit.GeV();
            if( chit->View() == geo::kY ) EprongY += rhit.GeV();
          }
          ProngEnergy.push_back(Eprong);
          ProngEnergyXView.push_back(EprongX);
          ProngEnergyYView.push_back(EprongY);
        }//loop all prongs
        
        if (nprongs > 1){  
          int Index1 = -1;
          int Index2 = -1;
          double MaxEnergy1 = -999.0;
          double MaxEnergy2 = -999.0;
          for( int i=0;i<nprongs; ++i ){
            if( MaxEnergy1<ProngEnergy[i] ){
              MaxEnergy1 = ProngEnergy[i];
              Index1     = i;
            }
          }//leading
          for( int i=0;i<nprongs; ++i ){
            if( i==Index1 ) continue;
            if( MaxEnergy2<ProngEnergy[i] ){
              MaxEnergy2=ProngEnergy[i];
              Index2=i;
            }
          }//sub-leading
          double Eden = ProngEnergy[Index1]+ProngEnergy[Index2];
          double Enum = ProngEnergy[Index1]-ProngEnergy[Index2];
          if( Eden>0.0 ) Ebalance=Enum/Eden;
        }//has at least 2 prongs
        //study 2D prongs
        
        std::vector<int>    Prong2DInXView;
        std::vector<double> ProngEnergy2D;
        int nprongs2D=SelectedProng2D.size();
        for(int ip=0; ip<nprongs2D; ++ip ){
          const art::Ptr<FUZZYKVERTEXOUTPUT> current_prong = SelectedProng2D[ip];

          int prongview=SelectedProng2D[ip]->View();
          Prong2DInXView.push_back(prongview);
          
          //reco momentum
          double Eprong = 0.0;
          
          const unsigned int nselectedprong2d_cells = current_prong->NCell();
          for( unsigned int icell=0; icell<nselectedprong2d_cells; ++icell){
            const rb::RecoHit rhit = current_prong->RecoHit(icell);

            if( !rhit.IsCalibrated() ) continue;
            Eprong += rhit.GeV();
          }
          ProngEnergy2D.push_back(Eprong);
        }

        int    index3D = -1;
        double max3D   = -99.0;
        for( int i=0; i<nprongs; ++i ){
          if( max3D<ProngEnergy[i] ){
            max3D   = ProngEnergy[i];
            index3D = i;
          }
        }//loop 3D prongs
        int    index2D = -1;
        double max2D   = -99.;
        for( int i=0;i<nprongs2D; ++i ){
          if( max2D<ProngEnergy2D[i] ){
            max2D   = ProngEnergy2D[i];
            index2D = i;
          }
        }//loop 2D prongs
        double ebalance2D = 1.0;
        if( nprongs>0 && nprongs2D>0 ){
          double Eden = 0.0;
          double Enum = 0.0;
          if( Prong2DInXView[index2D]==0 ){//XZ view
            Eden = ProngEnergy2D[index2D]+ProngEnergyXView[index3D];
            Enum = ProngEnergy2D[index2D]-ProngEnergyXView[index3D];
          }
          else {//YZ view
            Eden = ProngEnergy2D[index2D]+ProngEnergyYView[index3D];
            Enum = ProngEnergy2D[index2D]-ProngEnergyYView[index3D];
          }
          if( Eden > 0.0 ) ebalance2D=fabs(Enum)/Eden;
        }
        
        fInputVarsRVP[0]  = (Float_t)nCells;
        fInputVarsRVP[1]  = (Float_t)recoEnergy;
        fInputVarsRVP[2]  = (Float_t)longestTrack;
        if (nCells > 0) fInputVarsRVP[3]  = (Float_t)longestTrack_ncells/nCells;
        else fInputVarsRVP[3] = -1.0;
        if (Ntotal > 0) fInputVarsRVP[4]  = (Float_t)nmip/Ntotal;
        else fInputVarsRVP[4] = -1.0;
        fInputVarsRVP[5]  = (Float_t)nmip;
        //make sure total rec energy is nonzero
        if (recoEnergy > 0.0) fInputVarsRVP[6] = (Float_t)Eplane20/recoEnergy;
        else                  fInputVarsRVP[6] = -1.0;
        fInputVarsRVP[7]  = (Float_t)efrac_2slides;
        fInputVarsRVP[8]  = (Float_t)efrac_6slides;
        if (recoEnergy > 0.0) fInputVarsRVP[9] = (Float_t)totalE_road2sigma/recoEnergy;
        else                  fInputVarsRVP[9] = -1.0;
        fInputVarsRVP[10] = (Float_t)nprongs;
        fInputVarsRVP[11] = (Float_t)Ebalance;
        fInputVarsRVP[12] = (Float_t)nprongs2D;
        fInputVarsRVP[13] = (Float_t)ebalance2D;
        if (Esig_3     > 0.0) fInputVarsRVP[14] = (Float_t)(recoEnergy - Esig_3)/Esig_3;
        else                  fInputVarsRVP[14] = -1.0;

        if(fRunTMVA){
          BDT15 = fbdtRVP->EvaluateMVA(fbdtRVPMethodName);
          BDT12 = fbdtRVP12->EvaluateMVA(fbdtRVP12MethodName);
        }

        /// Slice for the PID
        const art::Ptr<rb::Cluster> bestSlice(slices,sliceIdx);

        RVP pid;
        pid.SetVal         (BDT15);
        pid.SetPdg         (12);
        pid.SetNCells      (fInputVarsRVP[0]);
        pid.SetRecoE       (fInputVarsRVP[1]);
        pid.SetLongTr      (fInputVarsRVP[2]);
        pid.SetLongTrFrac  (fInputVarsRVP[3]);
        pid.SetMipFrac     (fInputVarsRVP[4]);
        pid.SetMipHits     (fInputVarsRVP[5]);
        pid.SetEPl20Frac   (fInputVarsRVP[6]);
        pid.SetEFrac2PlWin (fInputVarsRVP[7]);
        pid.SetEFrac6PlWin (fInputVarsRVP[8]);
        pid.SetEFrac2SigRd (fInputVarsRVP[9]);
        pid.SetProngs3D    (fInputVarsRVP[10]);
        pid.SetProngEBal3D (fInputVarsRVP[11]);
        pid.SetProngs2D    (fInputVarsRVP[12]);
        pid.SetProngEBal2D (fInputVarsRVP[13]);
        pid.SetEIso3Sig    (fInputVarsRVP[14]);
        pid.SetRVP12       (BDT12);
          
        pidcol->push_back(pid);
        util::CreateAssn(*this,evt,*(pidcol.get()),bestSlice,*(assncol.get()));
        
        
        
      }//check Assn between Vertex and Slice exists
    }//has reco slicer
    
    /// If producing RVP, move products
    evt.put(std::move(pidcol));
    evt.put(std::move(assncol));
    
  }// end of produce

  
  //...........................................................................
  bool RecVarPID::bookVariables(const std::string& weight_file, TMVA::Reader& rvp, std::string& method_name)
  {
    /// Get all variable names from the weight XML
    const std::vector<std::string> varnames_from_xml = getVarNamesFromXML(weight_file);
    
    /// Number of variables
    const unsigned int nvars = varnames_from_xml.size();
    
    /// Add variables to TMVA
    for(unsigned int i=0; i<nvars; ++i){
      rvp.AddVariable(varnames_from_xml[i], &fInputVarsRVP[i]);
    }// end of adding variables to TMVA
    
    /// Extract method name and add it to the current name
    method_name += getMethodNameFromXML(weight_file);
    
    /// Book TMVA
    rvp.BookMVA(method_name, weight_file);
    
    return true;
  }// end of RecVarPID::bookVariables(TMVA::Reader& rvp)
  
  
  //---------------------------------------------------------------------------
  // Parse XML file to extract TMVA method name
  std::string RecVarPID::getMethodNameFromXML(const std::string weights_file_name,
                                              const bool        get_full_name
                                              ) const
  {
    
    /// System command to read out variables
    std::string cmd = "cat " + weights_file_name;
    cmd += "| grep \"<MethodSetup\"";
    
    /// d
    std::string data; ///< stream from system
    
    /// Perform a system call of cmd with the output stream going to data
    {
      FILE *stream;
      char buffer[MAX_BUFFER];
      
      stream = popen(cmd.c_str(), "r");
      
      while ( fgets(buffer, MAX_BUFFER, stream) != NULL ){
        data.append(buffer);
      }
      pclose(stream);
    }// end of the system call
    
    /// Replace stuff
    boost::replace_all(data, "<MethodSetup Method=", "");
    boost::replace_all(data, ">"                   , "");
    boost::replace_all(data, "\""                  , "");
    boost::replace_all(data, " "                   , "");
    boost::replace_all(data, "\n"                  , "");
    
    /// If getting a full name, return now
    ///     Full name is, for instance, "BDT::BDTG"
    ///     Short name would be "BDTG"
    if(get_full_name) return data;
    
    
    /// Getting the short name from the full name
    
    std::vector<std::string> names;
    
    /// Tokenize the string, with a separator ";"
    const boost::char_separator<char> sep(":");
    boost::tokenizer< boost::char_separator<char> > tokens(data, sep);
    BOOST_FOREACH (const std::string& t, tokens) {
      
      /// If the size is non-zero, add the string to the vector of arguments
      if(t.size() > 0)
        names.push_back(t);
    }// end of loop over tokens

    std::string out_string;
    
    if(names.size() > 0)
      out_string = names[names.size() - 1];

    return out_string;
  }
  
  //---------------------------------------------------------------------------
  // Parse XML file to extract variable names
  std::vector<std::string> RecVarPID::getVarNamesFromXML(const std::string file_name) const
  {
    /// System command to read out variables
    std::string cmd = "cat " + file_name;
    cmd += "| grep \"Variable VarIndex=\" | awk \'{print $3}\'";
    
    /// d
    std::string data; ///< stream from system
    
    /// Perform a system call of cmd with the output stream going to data
    {
      FILE *stream;
      char buffer[MAX_BUFFER];
      
      stream = popen(cmd.c_str(), "r");
      
      while ( fgets(buffer, MAX_BUFFER, stream) != NULL ){
        data.append(buffer);
      }
      pclose(stream);
    }// end of the system call

    /// Replace stuff
    boost::replace_all(data, "Expression=", ";");
    boost::replace_all(data, "\""         , "");
    boost::replace_all(data, "\n"         , "");
    boost::replace_all(data, ","          , "");
    boost::replace_all(data, " "          , "");
    /// After the replacement the string of variables have variables separated by ";"
    
    /// Output vector of variables in the XML file
    std::vector<std::string> vars;
    
    /// Tokenize the string, with a separator ";"
    const boost::char_separator<char> sep(";");
    boost::tokenizer< boost::char_separator<char> > tokens(data, sep);
    BOOST_FOREACH (const std::string& t, tokens) {
      
      /// If the size is non-zero, add the string to the vector of arguments
      if(t.size() > 0)
        vars.push_back(t);
    }// end of loop over tokens
        
    return vars;
  }
  
  DEFINE_ART_MODULE(RecVarPID)
  
}//end namespace