TrackCleanUpAlg.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file  TrackCleanUpAlg.cxx                                                                              
/// \brief Remove raw digits corresponding to muon hits.                                                   
///                                                                                                         
/// \version $Id: TrackCleanUpAlg.cxx,v 1.11 2012-10-26 22:19:03 bckhouse Exp $                             
/// \author  Kanika Sachdev (ksachdev@physics.umn.edu)                  
////////////////////////////////////////////////////////////////////////

#include "messagefacility/MessageLogger/MessageLogger.h"

#include "MuonRemove/art/TrackCleanUpAlg.h"

#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/WeightedHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "Utilities/func/MathUtil.h"

#include "TGeoManager.h"
#include "TGeoNode.h"

namespace murem
{

  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                                   
  
  TrackCleanUpAlg::TrackCleanUpAlg(fhicl::ParameterSet const& pset)
  {
    this->reconfigure(pset);
  }
  
  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////               


  TrackCleanUpAlg::~TrackCleanUpAlg()
  {
  
  }
  
  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                                   
  
  void TrackCleanUpAlg::reconfigure(const fhicl::ParameterSet& pset)
  {
    fMipRangeHigh                = pset.get< float >("MipRangeHigh");
    fMipRangeLow                 = pset.get< float >("MipRangeLow");
    fMipValue                    = pset.get< float >("MipValue");
    fVertexRegionDeDxCutOff      = pset.get< float >("VertexRegionDeDxCutOff");

    fMipRangeHighInGev           = pset.get< float >("MipRangeHighInGev");
    fMipRangeLowInGev            = pset.get< float >("MipRangeLowInGev");
    fMipValueInGev               = pset.get< float >("MipValueInGev");
    fVertexRegionDeDxCutOffInGev = pset.get< float >("VertexRegionDeDxCutOffInGev");
  }
  
  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                                   
  
  /// Track Clean up close to the vertex.                                                                  
  /// For hits on the track, checks if the energy per length                                               
  /// in the first few planes is consistent with a MIP.                                                     
  /// If it is not, modifies the hits and or the number of hits                                            
  /// associated with the track in a given plane.                                                          
  
  std::vector<rb::WeightedHit> TrackCleanUpAlg::CleanUpTrack(const rb::Track &muonTrack, 
                                                             const rb::Cluster &slice,
                                                             double*  mipRangeHigh,
                                                             double*  mipRangeLow,
                                                             double*  mipValue,
                                                             double*  vertexRegionDeDxCutOff) 
  {
                
    // Check if the mip numbers are provided to us. If not, use 
    // defaults which are good only for muons.

    if(mipRangeHigh)
      fMipRangeHigh = *mipRangeHigh;
    if(mipRangeLow)
      fMipRangeLow  = *mipRangeLow;
    if(mipValue)
      fMipValue     = *mipValue;
    
    // NB: everywhere in this script energy is calculated in units of MIP
    
    fExtraEOnTrack       = 0;
    fExtraEOnTrackPlanes.clear();

    double mipDeDx       = fMipValue;
    bool xViewCleaned    = false;
    bool yViewCleaned    = false;
    int vertexRegion[2];
    
    //if(vertexRegionDeDxCutOff)
    
    ComputeVertexRegion(muonTrack, vertexRegion, vertexRegionDeDxCutOff);
    
    int nPlanes = fPlanes.size();

    std::vector<rb::WeightedHit > trackHits;
    
    for(int i = 0; i < (int)muonTrack.NCell(); i++){
      rb::WeightedHit newHit( muonTrack.Cell(i), 1 );
      trackHits.push_back( newHit );
    }

    for( int plidx = 0; plidx < nPlanes; plidx++ ){      
      if ( xViewCleaned && yViewCleaned )
        break;

      int iPlane       = fPlanes[plidx];
      geo::View_t view = fgeom->Plane( iPlane )->View();

      if ( view == geo::kX){
        if (iPlane > vertexRegion[0] ){
          xViewCleaned = true;
          continue;
        }
      }
      else{
        if ( iPlane > vertexRegion[1] ){
          yViewCleaned = true;
          continue;
        }
      }

      // Track clean is only done for the part of the track impacted                                                     
      // by vertex energy. There can be variations in the dEdX beyond                                                    
      // the MIP RMS due to delta rays escaping the cell etc. We want                                                    
      // to make sure we don't have to deal with that by restricting                                                     
      // ourselves to the first half.                                                                                    

      // This is the energy a MIP should deposit in this plane.                                                         
      double eDep = (mipDeDx * fPlaneTrkLength[iPlane]);
  
      int flag = IsMIP( fPlaneDeDx[iPlane] );

      if(flag == 2){ //More than MIP                                                                       

        SortByDistFromTrack( fPlaneHits[iPlane], muonTrack, iPlane );
        
        float trialE = 0;
        int nPlaneHits = fPlaneHits[iPlane].size();
        for( int iHit = 0; iHit< nPlaneHits; iHit++){
          geo::OfflineChan p( iPlane, fPlaneHits[iPlane].at(iHit)->Cell() ) ;
          trialE += fCellEnergy[p];

          int cellFlag = IsMIP( trialE/fPlaneTrkLength[iPlane]) ;

          if (cellFlag == 0){ // Not enough energy to make a MIP yet.
            
            rb::WeightedHit newHit( fPlaneHits[iPlane].at(iHit), 1);
            fCleanUpWeightedHits.push_back(newHit);

            continue;
          }
          else{
            // If we are not less than MIP, erase all
            // other plane hits from trackHits
              
            for( int k = iHit+1; k<nPlaneHits; k++){
              float weight = 0;
              trackHits = ResetHitWeight( trackHits, fPlaneHits[iPlane].at(k), weight );
              geo::OfflineChan pa( iPlane, fPlaneHits[iPlane].at(k)->Cell() ) ;
              rb::WeightedHit newHit( fPlaneHits[iPlane].at(k), weight);
              fCleanUpWeightedHits.push_back(newHit);
              fExtraEOnTrack               += fCellEnergy[pa];
              fExtraEOnTrackPlanes[iPlane] += fCellEnergy[pa];

            }
                      
            if (cellFlag == 2){
              double muEnergyInCell = eDep + fCellEnergy[p] - trialE;
              float weight = muEnergyInCell/fCellEnergy[p];
              trackHits = ResetHitWeight( trackHits, fPlaneHits[iPlane].at(iHit), weight );
              rb::WeightedHit newHit(fPlaneHits[iPlane].at(iHit), weight);
              fCleanUpWeightedHits.push_back(newHit);
              fExtraEOnTrack               += fCellEnergy[p] - muEnergyInCell;
              fExtraEOnTrackPlanes[iPlane] += fCellEnergy[p] - muEnergyInCell;

            }

            break;
          }

        }// end loop over fPlaneHits 
      }// end if on energy in plane is greater than MIP                                                    
      
      

      if( flag == 0){ //Less than MIP

        int nPlaneHits = fPlaneHits[iPlane].size();     
        unsigned int cell[ nPlaneHits ];
        std::vector< art::Ptr<rb::CellHit> > slicePlaneHits;

        int nCells     = slice.NCell();

        for( int h = 0; h < nPlaneHits; h++)
          cell[h] = fPlaneHits[iPlane].at(h)->Cell();

        for( int hit = 0; hit < nCells; hit++){
          for(int iHit = 0; iHit < nPlaneHits; iHit++){
            if (slice.Cell(hit)->Plane() == iPlane && 
                !(*(slice.Cell(hit)) == *( fPlaneHits[iPlane].at(iHit) )) ) 
              slicePlaneHits.push_back(slice.Cell(hit));
          }
        }// end loop over slice hits                                                                       

        SortByDistFromTrack( slicePlaneHits, muonTrack, iPlane );
        float trialE = fPlaneEnergy[iPlane];
        
        int nSlicePlaneHits = slicePlaneHits.size();
        for(int hit = 0; hit< nSlicePlaneHits; hit++){

          bool closeEnough = false;
          for(int c = 0; c < nPlaneHits; c++ ){
            int diff = slicePlaneHits[hit]->Cell()  - cell[c];
            if ( diff < 5 && diff > -5 )
              closeEnough = true;
          }

          if( closeEnough){

            // Need to add each slicePlaneHit to muonTrack so that it is
            // calibrated accordingly

            rb::RecoHit rHit = muonTrack.RecoHit(slicePlaneHits[hit] );

            if( rHit.IsCalibrated() ){
              trialE += rHit.MIP();

              int cellFlag = IsMIP( trialE/fPlaneTrkLength[iPlane] ) ;
              geo::OfflineChan p( iPlane, slicePlaneHits[hit]->Cell() ) ;

              if (cellFlag == 0){
                rb::WeightedHit newHit(slicePlaneHits[hit], 1);

                trackHits.push_back(newHit);
                fCleanUpWeightedHits.push_back(newHit);
                fExtraEOnTrack               -= fCellEnergy[p];
                fExtraEOnTrackPlanes[iPlane] -= fCellEnergy[p];
              }   

              if (cellFlag == 1) {
                rb::WeightedHit newHit(slicePlaneHits[hit], 1);

                trackHits.push_back(newHit);
                fCleanUpWeightedHits.push_back(newHit);
                fExtraEOnTrack               -= fCellEnergy[p];
                fExtraEOnTrackPlanes[iPlane] -= fCellEnergy[p];
                break;
              }

              else if (cellFlag == 2){

                float tempEnergy = eDep + rHit.MIP() - trialE;

                rb::WeightedHit tempHit(slicePlaneHits[hit], tempEnergy/rHit.MIP() );
                fExtraEOnTrack               -= tempEnergy;
                fExtraEOnTrackPlanes[iPlane] -= tempEnergy;
                trackHits.push_back( tempHit );
                break;
              }

            
            }
          }// end of closeEnough                                                                           
        }// end loop over slice hits in iPlane                                                             
      }// end if energy in plane is less than MIP                                                          
    
  
    }// end loop over iPlanes.                                                                               
    
    return trackHits;
  }// end function CleanUpTrack                                                                              



  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                                   

  /// Checks if the given dE/dx is consistent with a MIP.                                                    
  /// Returns 1 if dE/dx is consistent with MIP.                                                             
  /// Returns 0, if dE/dx is less than MIP.                                                                  
  /// Returns 2, if dE/dx is more than MIP.                                                                  

  int TrackCleanUpAlg::IsMIP( double dEdX ) const
  {
    // This will probably be changed to compare to MIP scale as                                              
    // returned by recohit, please put up with the hard coded           
    // numbers for now! Refer NOVA-doc-7658-v2, slide 11 for origin.                                         
    if ( dEdX < fMipRangeLow ) return 0;     
    // Less than MIP
    if ( dEdX >= fMipRangeLow && dEdX <= fMipRangeHigh ) return 1;
    // Consistent with MIP
    if ( dEdX > fMipRangeHigh ) return 2;
    // More than MIP
    else return -1;
  }//end function IsMIP                                                                                      


  /// Finds and stores information plane by plane for a track.                                               
  /// The information it stores is energy per plane,                                                         
  /// track length in plane, dE/dx in plane, energy per cell                                                 
  /// and calibration constant for every cell on track.                                                      


  std::map< unsigned int, double> TrackCleanUpAlg::DeDxInPlane(const rb::Track &muonTrack)
  {

    fPlaneEnergy.clear();
    fPlaneTrkLength.clear();
    fPlaneDeDx.clear();
    fPlanes.clear();
    fPlanesX.clear();
    fPlanesY.clear();
    fPlaneHits.clear();

    int minPlane = muonTrack.MinPlane();
    int maxPlane = muonTrack.MaxPlane();

    art::PtrVector<rb::CellHit> trkHits = muonTrack.AllCells();
    SortByPlaneAndCell(trkHits);
    

    // Check if the track is vertical. It's all very easy in that case                                                                           
    if(minPlane == maxPlane){
      fPlanes.push_back(minPlane);
      if( fgeom->Plane( minPlane)->View() == geo::kX)
        fPlanesX.push_back( minPlane );
      else
        fPlanesY.push_back( maxPlane );
      
      double planeWidth = 2.0*fgeom->Plane(minPlane)->Cell(0)->HalfD();
      // MIP() is calculated in one cell in z, calculate length in terms of cell width
      fPlaneTrkLength[minPlane] = muonTrack.TotalLength()/planeWidth; 
      fPlaneEnergy[minPlane] = TrackEinMIP(muonTrack);
      fPlaneDeDx[minPlane] = fPlaneEnergy[minPlane]/fPlaneTrkLength[minPlane];
      for( auto const& iHit : trkHits){
        fPlaneHits[ minPlane ].push_back( iHit );
        rb::RecoHit rhit = muonTrack.RecoHit(iHit);
        if( rhit.IsCalibrated() )
          fCellEnergy[ iHit->OfflineChan()] = rhit.MIP();
      }
    }
    
    else{

      

      // Find all the unique planes that the track passes through
      for( auto const& iHit : trkHits){
        int iPlane = iHit->Plane();     
        fPlaneHits[ iPlane ].push_back( iHit );
        rb::RecoHit rhit = muonTrack.RecoHit(iHit);
        if( rhit.IsCalibrated() )
          fCellEnergy[ iHit->OfflineChan()] = rhit.MIP();
  
        if( std::find( fPlanes.begin(), fPlanes.end(), iPlane) == fPlanes.end() )
          // This plane is a new one! 
          fPlanes.push_back( iPlane );
        if( fgeom->Plane( iPlane)->View() == geo::kX)
          fPlanesX.push_back( iPlane );
        else
          fPlanesY.push_back( iPlane );
        
        if( rhit.IsCalibrated() )
          fPlaneEnergy[ iPlane ] += rhit.MIP(); 
      }// end loop over track hits
     
      // Now loop through the planes, find path-length through each and
      // store the dE/dx. 


      for(int iPl = 0, nPlanes = fPlanes.size(); iPl < nPlanes; ++iPl){

        int iPlane = fPlanes[iPl];

        // We can give exaggerated results for vertical or close to vertical
        // tracks, so just check the z direction cosine 

        if(  muonTrack.Dir().Z()  < 0.2 && muonTrack.Dir().Z() > -0.2 ){
          // Looks like we have a vertical-ish track on our hands.
          // We will find the max and min cells in 
          // double cellWidth = 2*fgeom->Plane(iPlane)->Cell(0)->HalfD();
          int ncell = fPlaneHits[ iPlane ].back()->Cell() - fPlaneHits[ iPlane ].front()->Cell();
          // MIP() is calculated in one cell in z, calculate length in terms of cell width
          fPlaneTrkLength[iPlane] = ncell - 1; // was cellWidth*( ncell - 1) previously
          // the minus 1 here is assuming that on average only
          // the half os the first and last cells in the plane
          // are traversed.
        }
        else{
          double planeHalfWidth = fgeom->Plane(iPlane)->Cell(0)->HalfD();
          double cellCenter[3];
          fgeom->Plane( iPlane )->Cell( 10 )->GetCenter(cellCenter);
          double xbound1, ybound1, xbound2, ybound2;
          muonTrack.InterpolateXY( cellCenter[2]-planeHalfWidth, xbound1, ybound1);
          muonTrack.InterpolateXY( cellCenter[2]+planeHalfWidth, xbound2, ybound2);
          TVector3 bound1( xbound1, ybound1, cellCenter[2]-planeHalfWidth);
          TVector3 bound2( xbound2, ybound2, cellCenter[2]+planeHalfWidth);

          // MIP() is calculated in one cell in z, calculate length in terms of cell width 
          fPlaneTrkLength[iPlane] = (bound2 - bound1).Mag()/planeHalfWidth/2.0;
        }
        
        fPlaneDeDx[iPlane] = fPlaneEnergy[iPlane]/fPlaneTrkLength[iPlane];
      }// end loop over planes

    }
    
    return fPlaneDeDx;
  }//End of DeDxInPlane                                                                                      
  

  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                                   

  /// Erases a specified hit from a std::vector of rb::CellHit.                                              

  std::vector<rb::WeightedHit > TrackCleanUpAlg::ResetHitWeight( std::vector< rb::WeightedHit> &trackHits,
                                                                 art::Ptr<rb::CellHit>         &hit,
                                                                 float                         &weight)
  {
    int nTrackHits = trackHits.size();
    for(int iHit = 0; iHit < nTrackHits ; iHit++){
      if(*(trackHits[iHit].hit) == *(hit))
        trackHits[iHit].weight = weight;
    }
    return(trackHits);
  }// End of EraseHit                                                                                        



  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////       

  // A function to compute the extent of hadronic energy impacted region                                     
  // along the length of the muon track.                                                                     

  void TrackCleanUpAlg::ComputeVertexRegion(const rb::Track& muonTrack, 
                                            int    *vertexRegion,
                                            double* vertexRegionDeDxCutOff) 
  {
    if(vertexRegionDeDxCutOff)
      fVertexRegionDeDxCutOff = *vertexRegionDeDxCutOff;
    
    if( !(muonTrack.NCell() > 0) )
      return;
    
    fAvgDeDx.clear();
    DeDxInPlane(muonTrack);

    vertexRegion[0]   = -1;
    vertexRegion[1]   = -1;


    for( int view = 0; view < 2 ; view++){
      std::vector<int> planes;

      if ( view == 0 )
        planes = fPlanesX;
      else 
        planes = fPlanesY;
      
      int nPlanes = planes.size();

      for(int ipl = 0; ipl < nPlanes - 2 ; ipl++ ){
        float avgDeDx = 0;
        for(int jpl = ipl; jpl < ipl+3; jpl++){
          avgDeDx +=  fPlaneDeDx[ planes[jpl] ];
        }
        
        fAvgDeDx[ planes[ipl] ] = avgDeDx/3;
      }
      
      if( nPlanes > 7){
        for(int ipl = 0; ipl < nPlanes - 3; ipl++ ){

          if( fAvgDeDx[ planes[ipl] + 1] < fVertexRegionDeDxCutOff && 
              fAvgDeDx[ planes[ipl] + 2] < fVertexRegionDeDxCutOff && 
              fAvgDeDx[ planes[ipl] + 3] < fVertexRegionDeDxCutOff && 
              fAvgDeDx[ planes[ipl] + 4] < fVertexRegionDeDxCutOff){
            vertexRegion[view] = planes[ipl];
            break;
          }
        }
        if( vertexRegion[view] == -1 && nPlanes>0 )
          vertexRegion[view] = planes[ nPlanes - 1 ];
      }
      if( nPlanes < 8){
        for(int ipl = 0; ipl < nPlanes - 3; ipl++ ){

          if( fPlaneDeDx[ planes[ipl] + 1] < fVertexRegionDeDxCutOff && 
              fPlaneDeDx[ planes[ipl] + 2] < fVertexRegionDeDxCutOff && 
              fPlaneDeDx[ planes[ipl] + 3] < fVertexRegionDeDxCutOff ){
            vertexRegion[view] = planes[ipl];
            break;
          }
        }
        if( vertexRegion[view] == -1 && nPlanes>0 )
          vertexRegion[view] = planes[ nPlanes - 1 ];
      }
    }
  }// End of ComputeVertexRegion       



  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                             
  
  
  void TrackCleanUpAlg::SortByDistFromTrack( std::vector< art::Ptr<rb::CellHit> >&hits, 
                                             const rb::Track &track, int plane)
  {
    
    int nHits = hits.size();

    for( int iHit = 0; iHit < nHits ; iHit++){

      TVector3 cellCenter, trkCenter;
      fgeom->Plane( hits[iHit]->Plane() )->Cell( hits[iHit]->Cell() )->GetCenter( cellCenter );

      // Find the z position of this plane, and use rb::Track::InterpolateXY to
      // figure the x,y position in this plane that the track passed through.
      fgeom->Plane( plane )->Cell( hits[iHit]->Cell() )->GetCenter( trkCenter );
      double x,y;
      track.InterpolateXY( trkCenter[2], x, y );
      TVector3 trkPoint( x, y, trkCenter[2]);

      float iCellDist = ( cellCenter - trkPoint ).Mag();

      for( int jHit = iHit+1; jHit < nHits; jHit++){
        fgeom->Plane( hits[jHit]->Plane() )->Cell( hits[jHit]->Cell() )->GetCenter( cellCenter );

        float jCellDist = ( cellCenter - trkPoint ).Mag();

        if( jCellDist < iCellDist ){
          art::Ptr<rb::CellHit> tmpCellHit = hits[iHit];
          hits[iHit] = hits[jHit];
          hits[jHit] = tmpCellHit;
          iCellDist  = jCellDist;
        }
      }//end inner loop to order hits in distance from track                                                      
    }//end outer loop to order hits in distance from track                                                        
    return;
  }


  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                             


  float TrackCleanUpAlg::ExtraEOnTrack(const rb::Track &track, 
                                       const rb::Cluster &slice,
                                       double*  mipRangeHigh,
                                       double*  mipRangeLow,
                                       double*  mipValue,
                                       double*  vertexRegionDeDxCutOff) 
  {
    CleanUpTrack( track, slice, mipRangeHigh, mipRangeLow, mipValue, vertexRegionDeDxCutOff);
    return fExtraEOnTrack;
  }

  float TrackCleanUpAlg::ExtraEOnTrackInGeV(const rb::Track &track,
                                            const rb::Cluster &slice,
                                            double*  mipRangeHigh,
                                            double*  mipRangeLow,
                                            double*  mipValue,
                                            double*  vertexRegionDeDxCutOff)
  {
    //std::cout<<"NB: entire TrackCleanUpAlg works in MIP units now, so your mipRange* values will be translated to the MIPs"<<std::endl;
    //std::cout<<"NB: if you haven't set mip arguments for this function default values from TrackCleanUpAlg.fcl will be called"<<std::endl;

    if(mipRangeHigh)
      fMipRangeHighInGev = *mipRangeHigh;
    if(mipRangeLow)
      fMipRangeLowInGev  = *mipRangeLow;
    if(mipValue)
      fMipValueInGev     = *mipValue;
    if(vertexRegionDeDxCutOff)
      fVertexRegionDeDxCutOffInGev     = *vertexRegionDeDxCutOff;
    
    double highInMIPs = fMipRangeHighInGev / fMipValueInGev; 
    double lowInMIPs = fMipRangeLowInGev / fMipValueInGev;
    double mipInMIPs = 1.0;
    double vertexRegionInMIPs = fVertexRegionDeDxCutOffInGev / fMipValueInGev;

    CleanUpTrack( track, slice, &highInMIPs, &lowInMIPs, &mipInMIPs, &vertexRegionInMIPs); // External packages will call this function with mipRange* in GeV/cm. Need to translate into MIP units
    return fExtraEOnTrack/94.62; // This a constant in Calibrator_service.cc, GetGeVToMIPScale function
  }

  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                             

  std::map<int, float> TrackCleanUpAlg::ExtraEOnTrackPlanes(const rb::Track &track,
                                                            const rb::Cluster &slice,
                                                            double*  mipRangeHigh,
                                                            double*  mipRangeLow,
                                                            double*  mipValue,
                                                            double*  vertexRegionDeDxCutOff) 
  {
    
    CleanUpTrack( track, slice, mipRangeHigh, mipRangeLow, mipValue, vertexRegionDeDxCutOff);
    return fExtraEOnTrackPlanes;
  }

  std::map<int, float> TrackCleanUpAlg::ExtraEOnTrackPlanesInGeV(const rb::Track &track,
                                                                 const rb::Cluster &slice,
                                                                 double*  mipRangeHigh,
                                                                 double*  mipRangeLow,
                                                                 double*  mipValue,
                                                                 double*  vertexRegionDeDxCutOff)
  {
    //std::cout<<"NB: entire TrackCleanUpAlg works in MIP units now, so your mipRange* values will be translated to the MIPs"<<std::endl;
    //std::cout<<"NB: if you haven't set mip arguments for this function default values from TrackCleanUpAlg.fcl will be called"<<std::endl;

    if(mipRangeHigh)
      fMipRangeHighInGev = *mipRangeHigh;
    if(mipRangeLow)
      fMipRangeLowInGev  = *mipRangeLow;
    if(mipValue)
      fMipValueInGev     = *mipValue;
    if(vertexRegionDeDxCutOff)
      fVertexRegionDeDxCutOffInGev     = *vertexRegionDeDxCutOff;

    double highInMIPs = fMipRangeHighInGev / fMipValueInGev;
    double lowInMIPs = fMipRangeLowInGev / fMipValueInGev;
    double mipInMIPs = 1.0;
    double vertexRegionInMIPs = fVertexRegionDeDxCutOffInGev / fMipValueInGev;

    CleanUpTrack( track, slice, &highInMIPs, &lowInMIPs, &mipInMIPs, &vertexRegionInMIPs); // External packages will call this function with mipRange* in GeV/cm. Need to translate into MIP units
    std::map<int, float> fExtraEOnTrackPlanesInGeV;
    for (auto itr = fExtraEOnTrackPlanes.begin(); itr != fExtraEOnTrackPlanes.end(); itr++){
         fExtraEOnTrackPlanesInGeV.insert ( std::pair<int, float>(itr->first, itr->second / 94.62) );  // This a constant in Calibrator_service.cc, GetGeVToMIPScale function
    }
    return fExtraEOnTrackPlanesInGeV; 
  }
  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                             

  void TrackCleanUpAlg::LastVertexPlanes(const rb::Track &track, 
                                         int    *vertexPlanes,
                                         double* vertexRegionDeDxCutOff) 
  {
    int vertexRegion[2];
    ComputeVertexRegion( track, vertexRegion, vertexRegionDeDxCutOff );
    vertexPlanes[0] = vertexRegion[0];
    vertexPlanes[1] = vertexRegion[1] ;
  }

  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                             
  std::vector<rb::WeightedHit> TrackCleanUpAlg::CleanUpWeightedHits(const rb::Track &track, 
                                                                    const rb::Cluster &slice,
                                                                    double*  mipRangeHigh,
                                                                    double*  mipRangeLow,
                                                                    double*  mipValue,
                                                                    double*  vertexRegionDeDxCutOff) 
  {
    
    CleanUpTrack( track, slice, mipRangeHigh, mipRangeLow, mipValue, vertexRegionDeDxCutOff);
    return fCleanUpWeightedHits;
  }

  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                             
  float TrackCleanUpAlg::TrackEinMIP (const rb::Track &track)
  {
    double ret = 0;
    for(unsigned int i = 0; i < track.NCell(); ++i){
        const rb::RecoHit rhit = track.RecoHit(i);
        if(rhit.IsCalibrated()) ret += rhit.MIP()*track.Weight(i);
    }
    return ret;

  }

}// End of namespace murem