TrackInfoFxs.cxx

Go to the documentation of this file.

#include "TrackInfo/TrackInfoFxs.h"

#include "Calibrator/Calibrator.h"
#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/RecoHit.h"

#include "art/Framework/Services/Registry/ServiceHandle.h"

#include "TMath.h"
#include "TVector3.h"

#include "MCCheater/BackTracker.h"
#include "CalibrationUtils/CalibUtil.h"

namespace trackinfo
{

  //-------------------------------------------------------
  TrackInfoFxs::TrackInfoFxs(){
  }

  //-------------------------------------------------------
  TrackInfoFxs::~TrackInfoFxs(){
  }

  //-----------------------------------------------------------------------------------

  float TrackInfoFxs::getAveTrackdEdxLastCells( const art::Ptr<rb::Track> track, unsigned int nLastCells){
    // This function returns the average dEdx in MeV/cm for the last nLastCells cells

    float avededxlast = -5; //default value

    if(track->NCell() < nLastCells){  //the number of cells in the track can not be less than the required number of last cells
      return avededxlast; 
    }
    else{
      art::ServiceHandle<cheat::BackTracker> btrk; 
      
      std::vector<std::vector<double> > vec_trkRecoHitInfo;
      std::pair<uint32_t, uint32_t> pc;
      TVector3 point;
      std::set<double> planeZBounds;
      calib::FindZBoundaries(planeZBounds);
      auto boundaryMap = calib::MakeZBoundaryMap(planeZBounds, track->Trajectory());
      if(boundaryMap.size() < 1) return avededxlast;
      //Loop over muon track hits                                                                                                 
      for(unsigned int jcell = 0; jcell < track->NCell();jcell++){
        const art::Ptr<rb::CellHit> cellhit = track->Cell(jcell);
        if(btrk->IsNoise(cellhit)) continue;
        pc.first = cellhit->Plane();                                                                                            
        pc.second = cellhit->Cell();                                                                                            
        
        rb::RecoHit rhit = track->RecoHit(cellhit);
        if(!rhit.IsCalibrated()) continue;
        
        std::vector<double> rhitInfo(5);
        rhitInfo[0] = rhit.X();
        rhitInfo[1] = rhit.Y();
        rhitInfo[2] = rhit.Z();
        point = TVector3(rhitInfo[0], rhitInfo[1], rhitInfo[2]);
        rhitInfo[3] = rhit.GeV();
        rhitInfo[4] = calib::PathLengthInCell(boundaryMap,point,pc);     

        vec_trkRecoHitInfo.push_back(rhitInfo);
      } //end loop over trk cell hits
      
        // using lambda function to sort it
      std::sort(vec_trkRecoHitInfo.begin(), vec_trkRecoHitInfo.end(),
                [](const std::vector<double>& a, const std::vector<double>& b) {
                  return a[2] < b[2];
                });
      
      float sumEnLast   = 0.0;
      float sumPathLast = 0.0;
      if (!vec_trkRecoHitInfo.empty()){
        if (nLastCells < 1 || vec_trkRecoHitInfo.size() < nLastCells){
          mf::LogWarning("TrackInfoFxs") << "nLastCells incompatible with reco hit info size";
          return avededxlast;
        }

        size_t min_index = vec_trkRecoHitInfo.size()-nLastCells;
        for(size_t a = min_index;a<vec_trkRecoHitInfo.size();a++) {
          sumEnLast   += (1000.*vec_trkRecoHitInfo[a][3]); // From GeV to MeV
          sumPathLast += vec_trkRecoHitInfo[a][4]; // Distance in cm 
        }
      }
      
      if(sumPathLast > 0){
        avededxlast = (sumEnLast)/(sumPathLast); // average dEdx in MeV/cm
      }      
    }
    return avededxlast;
  }//end of getAveTrackdEdxLastCells
    
  float TrackInfoFxs::getAveTrackdEdxEnd(std::vector<double> &dEtemp, std::vector<double> &dxtemp, 
                                         std::vector<double> &stemp,  double trk_len, double wanted_dist){
    float avededxend = -5; //default value
    if(trk_len <  wanted_dist){
      return avededxend;
    }
    
    double acc_dE = 0;
    double acc_dx = 0;
    //finding the minimum index:
    unsigned int min_index = 10000;
    for(unsigned int p = 0; p< stemp.size(); p++){
      double distcell_to_end = trk_len-stemp[p];
      if(distcell_to_end < 0){
        mf::LogWarning("TrackInfoFxs")<<"Something is wrong with this track, distance of a cell in the track to the end of the track is negative";
        return avededxend;
      }
      if(distcell_to_end < wanted_dist){
        min_index = p;
        break;
      }
    }
    
    for(unsigned int p = min_index; p < stemp.size(); p++){      
      if(dxtemp[p] <= 0){
        mf::LogWarning("TrackInfoFxs")<<"Something is wrong with dx... not positive value";
        return avededxend;
      }
      acc_dE += 1000. * dEtemp[p]; //in MeV
      acc_dx += dxtemp[p];         //in cm
    }  
    if(acc_dx>0)avededxend = acc_dE / acc_dx;

    return avededxend;
  }//end of getAveTrackdEdxEnd
  
} //end namespace