TCTrack.cxx

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
/// \brief   A TCHit contains pchits and information related to performing a timing calibration
/// \author  edniner@indiana.edu
/// \date    Feb. 2013
///////////////////////////////////////////////////////////////////////////// 
#include "CalibrationDataProducts/TCTrack.h"

#include "CalibrationDataProducts/PCHit.h"
#include "Utilities/func/MathUtil.h"

#include <cfloat>
#include <climits>
 
namespace caldp 
{
  TCTrack::TCTrack() :
  fTotalLength(0.0),
  fFiberVelocity(0.0)
  {  }
  
  
  TCTrack::~TCTrack(){  }

  //--------------------------------------------------------------------------
  void TCTrack::Add(const art::Ptr<caldp::PCHit>& pchit)
  {
    //if the nit view was never we want this to fail
    assert(pchit->View() != geo::kXorY);

    if (pchit->View() == geo::kX) fXPCHits.push_back(pchit);
    if (pchit->View() == geo::kY) fYPCHits.push_back(pchit);
  }

  //--------------------------------------------------------------------------
  void TCTrack::Add(const std::vector<art::Ptr<caldp::PCHit> >& pchits)
  {
    for (unsigned int n=0; n<pchits.size(); ++n){
      Add(pchits[n]);
    }
  }

  //--------------------------------------------------------------------------
  unsigned int TCTrack::NPCHit(geo::View_t view) const
  {
    switch(view){
      case geo::kX: return NXPCHit();
      case geo::kY: return NYPCHit();
      case geo::kXorY: return NPCHit();
      default: assert(0 && "Unknown view");
    }
  }

  //--------------------------------------------------------------------------
  art::Ptr<caldp::PCHit> TCTrack::PCHit(geo::View_t view, unsigned int viewIdx) const
  {
    switch(view){
      case geo::kX: return XPCHit(viewIdx);
      case geo::kY: return YPCHit(viewIdx);
      case geo::kXorY: return PCHit(viewIdx);
      default: assert(0 && "Unknown view");
    }
  }

  //--------------------------------------------------------------------------
  art::Ptr<caldp::PCHit> TCTrack::XPCHit(unsigned int xIdx) const
  {
    assert(xIdx < fXPCHits.size());
    return fXPCHits[xIdx];
  }

  //--------------------------------------------------------------------------
  art::Ptr<caldp::PCHit> TCTrack::YPCHit(unsigned int yIdx) const
  {
    assert(yIdx < fYPCHits.size());
    return fYPCHits[yIdx];
  }

  //--------------------------------------------------------------------------
  art::Ptr<caldp::PCHit> TCTrack::PCHit(unsigned int globalIdx) const
  {
    if (globalIdx < fXPCHits.size()) return XPCHit(globalIdx);
    return YPCHit(globalIdx - NXPCHit());
  }

  //--------------------------------------------------------------------------
  std::vector<art::Ptr<caldp::PCHit> > TCTrack::AllPCHits() const
  {
    std::vector<art::Ptr<caldp::PCHit> > all;
    for(size_t c=0; c<fXPCHits.size(); ++c) all.push_back(fXPCHits[c]);
    for(size_t c=0; c<fYPCHits.size(); ++c) all.push_back(fYPCHits[c]);
    return all;
  }

  //--------------------------------------------------------------------------
  void TCTrack::RemovePCHit(const art::Ptr<caldp::PCHit> pch)
  {
    if(pch->View() == geo::kX){
      for(unsigned int i=0; i<NXPCHit(); ++i){
        if(fXPCHits[i] == pch) fXPCHits.erase(fXPCHits.begin()+i);
        break;
      }
    }
    if(pch->View() == geo::kY){
      for(unsigned int i=0; i<NYPCHit(); ++i){
        if(fYPCHits[i] == pch) fYPCHits.erase(fYPCHits.begin()+i);
        break;
      }
    }
  }

  //--------------------------------------------------------------------------
  void TCTrack::RemovePCHit(unsigned int globalIdx)
  {
    assert(globalIdx < NPCHit());
    if(globalIdx < NXPCHit()) fXPCHits.erase(fXPCHits.begin()+globalIdx);
    else fYPCHits.erase(fYPCHits.begin()+globalIdx-NXPCHit());
  }

  //--------------------------------------------------------------------------
  int TCTrack::MinPlane(geo::View_t view) const
  {
    int ret = INT_MAX;
    for(unsigned int i=0; i<NPCHit(view); ++i){
      if(PCHit(view, i)->Plane() < ret) ret = PCHit(view, i)->Plane();
    }
    return ret;
  }

  //--------------------------------------------------------------------------
  int TCTrack::MinCell(geo::View_t view) const
  {
    int ret = INT_MAX;
    for(unsigned int i=0; i<NPCHit(view); ++i){
      if(PCHit(view, i)->Cell() < ret) ret = PCHit(view, i)->Cell();
    }
    return ret;
  }

  //--------------------------------------------------------------------------
  int TCTrack::MaxPlane(geo::View_t view) const
  {
    int ret = -INT_MAX;
    for(unsigned int i=0; i<NPCHit(view); ++i){
      if(PCHit(view, i)->Plane() > ret) ret = PCHit(view, i)->Plane();
    }
    return ret;
  }

  //--------------------------------------------------------------------------
  int TCTrack::MaxCell(geo::View_t view) const
  {
    int ret = -INT_MAX;
    for(unsigned int i=0; i<NPCHit(view); ++i){
      if(PCHit(view, i)->Cell() > ret) ret = PCHit(view, i)->Cell();
    }
    return ret;
  }

  //--------------------------------------------------------------------------
  double TCTrack::MinTNS(geo::View_t view) const
  {
    double ret = DBL_MAX;
    for(unsigned int i=0; i<NPCHit(view); ++i){
      if(PCHit(view, i)->TNS() < ret) ret = PCHit(view, i)->TNS();
    }
    return ret;
  }

  //--------------------------------------------------------------------------
  double TCTrack::MaxTNS(geo::View_t view) const
  {
    double ret = -DBL_MAX;
    for(unsigned int i=0; i<NPCHit(view); ++i){
      if(PCHit(view, i)->TNS() > ret) ret = PCHit(view, i)->TNS();
    }
    return ret;
  }

  //--------------------------------------------------------------------------
  double TCTrack::MinPath(geo::View_t view) const
  {
    double ret = DBL_MAX;
    for(unsigned int i=0; i<NPCHit(view); ++i){
      if(PCHit(view, i)->FlightLen() < ret) ret = PCHit(view, i)->FlightLen();
    }
    return ret;
  }

  //--------------------------------------------------------------------------
  double TCTrack::MaxPath(geo::View_t view) const
  {
    double ret = -DBL_MAX;
    for(unsigned int i=0; i<NPCHit(view); ++i){
      if(PCHit(view, i)->FlightLen() > ret) ret = PCHit(view, i)->FlightLen();
    }
    return ret;
  }

  //--------------------------------------------------------------------------
  double TCTrack::TNSUncertainty(float pe, bool goodTime, novadaq::cnv::DetId detectorID) const
  {

    //This is a fit for the timing resolution as a function of PE.  
    //Current values of the fit come from S14-04-14 data-like MC.  This fit can be found in
    //slide 4 of docdb-11252.
    //TODO: The constants in this function should be moved into the DB for ease of retrival

    // uncertainty derived separately for single and multipoint readout.  If timing is good, use multip, otherwise 
    // stay with the single point value.  These values are for the far detector, not checking detector type yet.

    double stdev = 144.0; //simple assumption of 500ns/sqrt(12)
    if (detectorID == novadaq::cnv::kFARDET) {
      if (goodTime){
        stdev = 231594.0/(2267.16+pow(pe,2.01011)) + 9.55689;   //2014-10-13 fit, multipoint readout from data, docdb-12122, with 460 ns rise time
      }
      else {
        stdev = 848357.0/(5734.75+pow(pe,2.20060)) + 142.221; //2014-06-13 fit, from data
      }
    }
    else if (detectorID == novadaq::cnv::kNEARDET) {
      if (goodTime){
        stdev = 2229.53/(77.0043+pow(pe,1.22743)) + 4.89355; //2014-09-15 fit from data with 140 ns rise time, 4.5 us fall time
      }
      else {
        stdev = 6136470.0/(185518.0+pow(pe,2.92360)) + 31.6071; //2014-09-15 singe point fit of data
      }
    }

    return stdev;
  }

  //--------------------------------------------------------------------------
  double TCTrack::ReadoutDistCorrection(float x, novadaq::cnv::DetId detectorID) const
  {

    //This is a fit for the timewalk as a function of distance of the hit to the readout electronics. 
    //This effect comes from the fiber being looped, causing early and late arrival of light which can distort the timing fit
    //The timewalk has been parameterized as an 8 degree polynomial which is ok since extrapolation is never needed.
    //Current values of this fit come from MC with multipoint readout looking at the difference between expected and reconstructed photon times.
    //This fit was made simultaneously with a fit of fiberspeed as a function of distance to readout.  More details can be found on slide 3 of docdb-11252. 
    //TODO: The constants in this function should be moved into the DB for ease of retrival
    //TODO: Incorporate timewalk for NearDet, where the effect will be different.  For now apply FD values for both detectors.

    double tw = 0.0;  //default to a value of 0 for NDOS
    if (detectorID == novadaq::cnv::DetId::kFARDET) {
      tw = -9.63368 
           + x*(0.0185403 
           + x*(-(1.29139e-4)
           + x*((3.98478e-7)
           + x*(-(7.66614e-10)
           + x*((8.9566e-13)
           + x*(-(6.2088e-16)
           + x*((2.37032e-19)
           + x*(-3.81279e-23))))))));
    }
    else if (detectorID == novadaq::cnv::DetId::kNEARDET) {
      /*tw = -9.63368 
           + x*(0.0185403 
           + x*(-(1.29139e-4)
           + x*((3.98478e-7)
           + x*(-(7.66614e-10)
           + x*((8.9566e-13)
           + x*(-(6.2088e-16)
           + x*((2.37032e-19)
           + x*(-3.81279e-23))))))));*/
    }
    return tw;

  }

  //--------------------------------------------------------------------------
  double TCTrack::CalcFiberVelocity(float d, novadaq::cnv::DetId detectorID) const
  {

    //The effective fiber speed has some dependence on the location of light within the cell.  Near the 
    //top of the cell the speed will be slower since light cannot travel both ways in scintillator before being absorbed
    //in the fiber.  This fit is done for now from multipoint MC and is documented in docdb-10863 and the final fit result is on slide 3 of docdb-11252
    //TODO: This fit is specific to the FD, will have to be redone for ND.  For now assume a flat fiber speed for ND
    double fspeed = 15.7446; //cm/ns
    if (detectorID == novadaq::cnv::DetId::kFARDET) {
      fspeed = 15.7446/(((15.7446*1.63517)/d)+1.0);
    }
    else if (detectorID == novadaq::cnv::DetId::kNEARDET) {
      //fspeed = 15.7446/(((15.7446*1.63517)/d)+1.0);
    }

    return fspeed;

  }

   //--------------------------------------------------------------------------
   double TCTrack::CalculateFiberVelocity(double timeResid, 
                                          std::vector<double> *resids,
                                          double *icept,
                                          double *chisqr,
                                          novadaq::cnv::DetId detectorID,
                                          bool removeHits,
                                          bool overwrite)
   {
     //define the speed of light in cm/ns
     const double sCmPerNsec = 29.9792458;

     //Containers to store pieces needed for fiber speed fit.
     //Times are the raw TNS for a cell - TOF correction
     //Where TOF is the time-of-flight from the start of the track to a hit
     std::vector<double> time;
     std::vector<double> rDist;
     std::vector<double> sterror;
     std::vector<double> weight;

     for(unsigned int i=0; i<NPCHit(); ++i){
       art::Ptr<caldp::PCHit> pch = PCHit(i);
       double tns;
       double pe = pch->PE();
       double x = pch->ReadDist();
       double tw = ReadoutDistCorrection(x, detectorID);

       tns = pch->TNS() + tw;
       time.push_back(tns - pch->FlightLen()/sCmPerNsec);
       rDist.push_back(pch->ReadDist());
       double stdev = TNSUncertainty(pe, pch->GoodTime(), detectorID);
       sterror.push_back(stdev);
       weight.push_back(1.0/(stdev*stdev));
     }
     double slope, intercept, x1, y1;
     double fspeed = -99.0;
     //do iterative fit for the speed light in the optical fiber
     while (true){
       //don't let fit go below two points
       if (rDist.size() < 2) break;
       resids->clear();
       *chisqr = util::LinFit(rDist, time, weight, slope, intercept);
       fspeed = slope;
       x1 = rDist[0];
       y1 = slope*x1 + intercept;
       double maxresid = 0;
       unsigned int maxpos = 0;
       for(unsigned int j=0; j<rDist.size(); ++j){
         double recresid = time[j] - (fspeed*(rDist[j]-x1)+y1);
         resids->push_back(recresid);
         if (fabs(recresid/sterror[j]) > fabs(maxresid)){
           maxresid = fabs(recresid/sterror[j]);
           maxpos = j;
         }
       }
       //if fit only has two points, dont remove any
       if (rDist.size() == 2) break;
       //if there was a hit with a residual above the maximum
       //remove that hit and recalculate
       if (fabs(maxresid) > timeResid){
         weight.erase(weight.begin()+maxpos);
         rDist.erase(rDist.begin()+maxpos);
         time.erase(time.begin()+maxpos);
         sterror.erase(sterror.begin()+maxpos);
         resids->erase(resids->begin()+maxpos);
         if (removeHits == true){
           fDropPCHits.push_back(PCHit(maxpos));
           RemovePCHit(maxpos);
         }
         continue;
       }
       //stop iterating if no timing residuals were too large
       break;
     }
     *icept = intercept;
     if (overwrite == true) fFiberVelocity = 1.0/fspeed;
     return 1.0/fspeed;
   }

   //--------------------------------------------------------------------------
   double TCTrack::CalculateMuonVelocity(std::vector<double> *resids,
                                         double *icept,
                                         double *chisqr,
                                         novadaq::cnv::DetId detectorID,
                                         double fiberVel) const
   {
     //define the speed of light in cm/ns
     const double sCmPerNsec = 29.9792458;



     //Containers to store pieces needed for fiber speed fit.
     //Times are the raw TNS for a cell  - Fiber speed correction
     //Where this correction is the time for light to travel up the optical fiber
     std::vector<double> time;
     std::vector<double> rDist;
     std::vector<double> weight;

     for(unsigned int i=0; i<NPCHit(); ++i){
       art::Ptr<caldp::PCHit> pch = PCHit(i);
       double tns;
       double pe = pch->PE();
       double x = pch->ReadDist();
       double tw = ReadoutDistCorrection(x, detectorID);

       //if a specific fibr speed is not specified, use parameterization as a function of distance to readout
       if (fiberVel == 9999.0) fiberVel = CalcFiberVelocity(x, detectorID);

       tns = pch->TNS() + tw;
       time.push_back(tns - pch->ReadDist()/fiberVel);
       rDist.push_back(pch->FlightLen());
       double stdev = TNSUncertainty(pe, pch->GoodTime(), detectorID);

       weight.push_back(1.0/(stdev*stdev));
     }

     double slope, intercept, x1, y1;
     double mspeed = -5.0;
     
     //do the fit one time
     resids->clear();
     //make sure there are at least two points before doing fit;
     if (rDist.size() < 2) return -5.0;
     *chisqr = util::LinFit(rDist, time, weight, slope, intercept);
     mspeed = slope;
     x1 = rDist[0];
     y1 = slope*x1 + intercept;
     for(unsigned int j=0; j<rDist.size(); ++j){
       double recresid = time[j] - (mspeed*(rDist[j]-x1)+y1);
       resids->push_back(recresid);
     }

     *icept = intercept;
     //return muon velocity in terms of 1/Beta
     return mspeed*sCmPerNsec;

   }

    
  
} // end namespace calib