dEdxCalculator.cxx

Go to the documentation of this file.

///
/// \file dEdxCalculator.cxx 
/// \brief Calculates dE/dx in cells passed through by a track.
///

#include "RecoBase/CellHit.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Track.h"

#include "BreakPointFitter/func/dEdxCalculator.h"
#include "BreakPointFitter/func/Path.h"

using namespace bpfit;

//static const double gsMmuon   = 105.658E-3;
//static const double gsMpion   = 139.570E-3;
//static const double gsMproton = 938.272E-3;

//......................................................................

dEdxCalculator::dEdxCalculator() { }

//......................................................................

void dEdxCalculator::computeDEDX(art::Ptr<rb::Track> &track, int pdg)
{

  // Reset.
  this->clearVectors();

  // get the cells in the track and get the trajectory points
  art::PtrVector<rb::CellHit> hits = track->AllCells();
  std::vector<TVector3>       traj = track->Trajectory();

  // keep track of trajectory length so far
  double stot = 0.0;
    
  //
  // For each track, we will follow along the trajectory segments and make
  // a list of all of the cells that the track passed through.  For each of
  // these cells, we will compute dx as the length between the entry point
  // and the exit point.  Since it is possible for a trajectory segment to
  // end inside a cell (with the next segment then continuing on to exit
  // the cell) we will sum the computed dx's for cells that appear on our
  // list more than once.
  //
    
  // Loop over trajectory points:
  for(unsigned int i = 0; i+1 < traj.size(); ++i) {

    // add length of previous segment to total path length
    if(i > 0) {
      double x0 = traj[i].X() - traj[i-1].X();
      double y0 = traj[i].Y() - traj[i-1].Y();
      double z0 = traj[i].Z() - traj[i-1].Z();
      double d0 = sqrt(x0*x0 + y0*y0 + z0*z0);
      stot += d0;
    }
        
    std::vector<geo::CellOnLine> Xhits;
    std::vector<geo::CellOnLine> Yhits;
        
    // make list of cells passed through by this segment of the trajectory
    fGeom->CountCellsOnLine(traj[i],traj[i+1],Xhits,Yhits);

    // Check if hits are already in my list.  If it is not on the list,
    // add it and its info.  If it is, then only add the additional
    // computed dx.
    bool onlist = false;
    unsigned int pos = 0;
    for(unsigned int j = 0; j < Xhits.size(); ++j) {

      onlist = false;
      pos    = 0;
      for(unsigned int k = 0; k < chan.size(); ++k) {
        if(Xhits[j].chan == chan[k]) {
          onlist = true;
          pos = k;
          break;
        }
      } // end loop over elements of chan
            
      // compute path length through the cell (for dx)
      double x = Xhits[j].entry.X() - Xhits[j].exit.X();
      double y = Xhits[j].entry.Y() - Xhits[j].exit.Y();
      double z = Xhits[j].entry.Z() - Xhits[j].exit.Z();
      double d = sqrt(x*x + y*y + z*z);
            
      // compute length from trajectory point i to entry point of cell
      // (for proper calculation of path length s for this cell)
      double x1 = Xhits[j].entry.X() - traj[i].X();
      double y1 = Xhits[j].entry.Y() - traj[i].Y();
      double z1 = Xhits[j].entry.Z() - traj[i].Z();
      double d1 = sqrt(x1*x1 + y1*y1 + z1*z1);

      if(onlist) {
        dx   [pos] += d;
        s    [pos] += d/2.0;
      }
      else {
        chan.push_back(Xhits[j].chan);
        dx.push_back(d);
        // The path length for this hit is taken to be the trajectory length
        // (up to the i-th trajectory point) plus the distance from the i-th
        // trajectory point to the entrance to the cell, plus half the distance
        // travelled through the cell.
        s.push_back(stot + d1 + d/2.0);
      }
            
    } // end loop over elements of Xhits
        
    onlist = false;
    pos = 0;
    for(unsigned int j = 0; j < Yhits.size(); ++j) {

      onlist = false;
      pos    = 0;
      for(unsigned int k = 0; k < chan.size(); ++k) {
        if(Yhits[j].chan == chan[k]) {
          onlist = true;
          pos = k;
          break;
        }
      } // end loop over elements of chan
            
      // compute path length through the cell
      double x = Yhits[j].entry.X() - Yhits[j].exit.X();
      double y = Yhits[j].entry.Y() - Yhits[j].exit.Y();
      double z = Yhits[j].entry.Z() - Yhits[j].exit.Z();
      double d = sqrt(x*x + y*y + z*z);
            
      // compute length from trajectory point i to entry point of cell
      // (for proper calculation of path length s for this cell)
      double x1 = Yhits[j].entry.X() - traj[i].X();
      double y1 = Yhits[j].entry.Y() - traj[i].Y();
      double z1 = Yhits[j].entry.Z() - traj[i].Z();
      double d1 = sqrt(x1*x1 + y1*y1 + z1*z1);

      if(onlist) {
        dx[pos]    += d;
        s [pos]    += d/2.0;
      }
      else {
        chan.push_back(Yhits[j].chan);
        dx.push_back(d);
        // The path length for this hit is taken to be the trajectory length
        // (up to the i-th trajectory point) plus the distance from the i-th
        // trajectory point to the entrance to the cell, plus half the distance
        // travelled through the cell.
        s.push_back(stot + d1 + d/2.0);
      }
            
    } // end loop over elements of Yhits
        
  } // end loop over trajectory points (i)
    
  // check that we haven't messed up
  if(chan.size() != dx.size() || chan.size() != s.size()) abort();

    
    
  // make the bpfit::path object to compute momentum along the trajectory
  bpfit::Path path(*fGeom, traj.size());
  for(unsigned int i = 0; i < traj.size(); ++i) {
    // set the points in path using the trajectory points
    double xyz[3];
    xyz[0] = traj[i].X();
    xyz[1] = traj[i].Y();
    xyz[2] = traj[i].Z();
    path.SetPoint(i,xyz);
  }
    
    
    
    
  //
  // Now populate the list of energies and momenta.  Since a trajectory can pass
  // through a cell without there being a cell hit, we will loop over all cells
  // that the trajectory passed through and only compute dE/dx for the ones that
  // have a corresponding cell hit on the track.
  //
    
  // loop over cells in chan and check for a cell hit in the list of hits in the track
  for(unsigned int i = 0; i < chan.size(); ++i) {

    art::Ptr<rb::CellHit> Ctemp = track->Cell(0);
    double dEtemp    = -1;
    double Wtemp     = -1e3;
    double ptemp     = -1;
    double Btemp     = -1;
    double Gtemp     =  1;
    double Ttemp     = -1;
    bool   match     = false;
    unsigned int pos = 0;
        
    // look for this cell in the list of hits in the track
    for(unsigned int j = 0; j < hits.size(); ++j) {
            
      // match by plane and cell number (can this assumption lead to trouble???)
      if(chan[i].Cell() == hits[j]->Cell() && chan[i].Plane() == hits[j]->Plane()) {
        match = true;
        pos = j;
        break;
      }
            
    } // end loop over hits on track (j)

    // A match was found, compute the real energy of the cell correcting for attenuation calibration.
    if(match) {
            
      // make reco hit from the cell hit
      const rb::CellHit *chit = track->Cell(pos).get();
      Ctemp = track->Cell(pos);
      Wtemp = track->W(chit);
      rb::RecoHit rhit(fCal->MakeRecoHit(*chit,Wtemp));
            
      // if reco hit is calibrated, get GeV to use as dE and compute the momentum using path
      if(rhit.IsCalibrated()) {
        dEtemp = rhit.GeV();
        
        if     (pdg == 13)   path.MuonParams(s[i], &Ttemp, &ptemp, &Btemp, &Gtemp);
        else if(pdg == 211)  path.PionParams(s[i], &Ttemp, &ptemp, &Btemp, &Gtemp);
        else if(pdg == 2212) path.ProtonParams(s[i], &Ttemp, &ptemp, &Btemp, &Gtemp);
        else {
          std::cerr << "\n\n\nERROR: PDG code " << pdg << " is not supported."
                    << "\tABORTING...";
          abort();
        }
        
      }

    }
        
    dE.push_back(dEtemp);
    W .push_back(Wtemp);
    bg.push_back(Btemp*Gtemp);
    cells.push_back(Ctemp);
        
  } // end loop over cells in chan list (i)

}

//......................................................................

void dEdxCalculator::getResults(std::vector<double> &dEtemp, std::vector<double> &dxtemp,
                                std::vector<double> &Wtemp,  std::vector<double> &bgtemp,
                                std::vector<double> &stemp,  double dxTol)
{

  // Construct the list of successfully calculated dE/dx values.
  for(unsigned int i = 0; i < chan.size(); ++i) {
    if(dE[i] > 0.0 && dx[i] > dxTol) {
      dEtemp.push_back(dE[i]);
      dxtemp.push_back(dx[i]);
      Wtemp .push_back(W[i]);
      bgtemp.push_back(bg[i]);
      stemp .push_back(s[i]);
    }
  }

}

//......................................................................

void dEdxCalculator::getDE(std::vector<double> &dEtemp, double dxTol)
{

  // Construct the list of successfully calculated values.
  for(unsigned int i = 0; i < chan.size(); ++i) {
    if(dE[i] > 0.0 && dx[i] > dxTol) {
      dEtemp.push_back(dE[i]);
    }
  }

}

//......................................................................

void dEdxCalculator::getDX(std::vector<double> &dxtemp, double dxTol)
{

  // Construct the list of successfully calculated values.
  for(unsigned int i = 0; i < chan.size(); ++i) {
    if(dE[i] > 0.0 && dx[i] > dxTol) {
      dxtemp.push_back(dx[i]);
    }
  }

}

//......................................................................

void dEdxCalculator::getW(std::vector<double> &Wtemp, double dxTol)
{

  // Construct the list of successfully calculated values.
  for(unsigned int i = 0; i < chan.size(); ++i) {
    if(dE[i] > 0.0 && dx[i] > dxTol) {
      Wtemp.push_back(W[i]);
    }
  }

}

//......................................................................

void dEdxCalculator::getBG(std::vector<double> &bgtemp, double dxTol)
{

  // Construct the list of successfully calculated values.
  for(unsigned int i = 0; i < chan.size(); ++i) {
    if(dE[i] > 0.0 && dx[i] > dxTol) {
      bgtemp.push_back(bg[i]);
    }
  }

}

//......................................................................

void dEdxCalculator::getS(std::vector<double> &stemp, double dxTol)
{

  // Construct the list of successfully calculated values.
  for(unsigned int i = 0; i < chan.size(); ++i) {
    if(dE[i] > 0.0 && dx[i] > dxTol) {
      stemp.push_back(s[i]);
    }
  }

}

//......................................................................

art::PtrVector<rb::CellHit> dEdxCalculator::getCellHits(double dxTol)
{
  
  art::PtrVector<rb::CellHit> returnedCells;

  // Construct the list of cells used to calculate the dE/dx values.
  for(unsigned int i = 0; i < chan.size(); ++i) {
    if(dE[i] > 0.0 && dx[i] > dxTol) {
      returnedCells.push_back(cells[i]);
    }
  }

  return returnedCells;

}

//......................................................................

void dEdxCalculator::clearVectors()
{

  // Clear all vectors used.
  chan .clear();
  cells.clear();
  dE   .clear();
  dx   .clear();
  W    .clear();
  s    .clear();
  bg   .clear();

}

////////////////////////////////////////////////////////////////////////