HitList3D.cxx

Go to the documentation of this file.

#include "BreakPointFitter/func/HitList3D.h"
#include <cstdlib>
#include <cmath>
#include <iostream>
#include <algorithm>
using namespace bpfit;

HitList3D::HitList3D() { }
HitList3D::HitList3D(unsigned int n) : std::vector<Hit3D>(n) { }

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

void HitList3D::AddHit(geo::View_t v,
                       double x,
                       double y,
                       double z, 
                       double t,
                       double sig_x,
                       double sig_y,
                       double sig_z,
                       double sig_t,
                       double w)
{
  this->push_back(Hit3D(v,x,y,z,t,sig_x,sig_y,sig_z,sig_t,w));
}

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

void HitList3D::SetHit(unsigned int i,
                       geo::View_t v,
                       double x,
                       double y,
                       double z, 
                       double t,
                       double sig_x,
                       double sig_y,
                       double sig_z,
                       double sig_t,
                       double w)
{
  Hit3D& h = (*this)[i];
  h.fView = v;
  h.fX    = x;
  h.fY    = y;
  h.fZ    = z;
  h.fT    = t;
  h.fSigX = sig_x;
  h.fSigY = sig_y;
  h.fSigZ = sig_z;
  h.fW    = w;
}

//......................................................................
///
/// Supply a default interpolation routine. This makes a staight line
/// fit to the (up to) 4 closest hits from the other view.
///
/// This function assumes the hits are sorted by z for optimum CPU
/// performance.
///
void gsDefaultInterp(Hit3D& h, HitList3D& hl) 
{
  unsigned int i, i0, i1, i2;
  unsigned int npts = 6;      // Minimum number of hits to use in interpolation
  double dz         = 100.0;  // Range in z to search for opposite view hits
  double dzmax      = 800.0;  // Largest z to search over
  
  for (i0=0; i0<hl.size(); ++i0) {
    if (&h==&hl[i0]) break;
  }

  unsigned int n = 0; // Number of hits used for interpolation
  double Sw=0, Swx=0, Swy=0, Swxy=0, Swx2=0, wz;
  while (dz<dzmax) {
    for (i1=i0; i1!=0; --i1) {
      if (fabs(h.fZ-hl[i1].fZ)>dz) { ++i1; break; }
    }
    for (i2=i0; i2<hl.size(); ++i2) {
      if (fabs(h.fZ-hl[i2].fZ)>dz) break;
    }
    n = 0;
    Sw = Swx = Swy = Swxy = Swx2 = 0;
    for (i=i1; i<i2; ++i) {
      //
      // Skip hits from the same view
      //
      if (h.fView==hl[i].fView) continue;
      if (fabs(h.fZ-hl[i].fZ)<dz) {
        ++n;
        Sw += hl[i].fW;
        wz = hl[i].fZ * hl[i].fW;
        Swx  += wz;
        Swx2 += wz*hl[i].fZ;
        if (hl[i].fView==geo::kX) {
          Swy += hl[i].fW*hl[i].fX;
          Swxy += wz*hl[i].fX;
        }
        else if (hl[i].fView==geo::kY) {
          Swy += hl[i].fW*hl[i].fY;
          Swxy += wz*hl[i].fY;
        }
        else abort();
      }
    }
    //
    // Test doneness
    //
    if (n>=npts) break; // Have enough points to interpolate
    //
    // Not done, open up the window and keep searching
    //
    dz *= 2;
  }
  //
  // Use what points we found to make a linear fit. In our coordinate
  // system "z" plays the role of x and "xy" plays the tole of y for a
  // fit y=m*x+b
  //
  double d = Sw*Swx2-Swx*Swx;
  if (d==0.0) {
    //
    // If d is zero we have no reliable information from the other
    // view. Put the hit at the detector center with large error.
    //
    if      (h.fView==geo::kX) { h.fY = 0.0; h.fSigY = 1000; }
    else if (h.fView==geo::kY) { h.fX = 0.0; h.fSigX = 1000; }
    return;
  }
  double m = (Sw*Swxy-Swx*Swy)/d;
  double b = -(Swx*Swxy-Swx2*Swy)/d;

  if (h.fView==geo::kX) {
    h.fY    = m*h.fZ+b;
    h.fSigY = sqrt(m*m*h.fSigZ*h.fSigZ + h.fSigX*h.fSigX);
  }
  else if (h.fView==geo::kY) {
    h.fX    = m*h.fZ+b;
    h.fSigX = sqrt(m*m*h.fSigZ*h.fSigZ + h.fSigY*h.fSigY);
  }
  else abort();
}

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

void HitList3D::CompleteOrtho(Interp3D interp)
{
  if (interp==0) {
    this->SortByZ();
    interp = gsDefaultInterp;
  }
  
  // ...for the sake of readability
  HitList3D& hitlist = (*this);
  
  unsigned int i;
  for (i=0; i<this->size(); ++i) {
    (*interp)(hitlist[i], hitlist);
  }
}

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

static bool gsSortHit3DByZ(const Hit3D& h1, const Hit3D& h2) 
{
  return h1.fZ<h2.fZ;
}

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

static bool gsSortHit3DByY(const Hit3D& h1, const Hit3D& h2) 
{
  return h1.fY>h2.fY;
}

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

void HitList3D::SortByZ() 
{
  sort(this->begin(), this->end(), gsSortHit3DByZ);
}

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

void HitList3D::SortByY() 
{
  sort(this->begin(), this->end(), gsSortHit3DByY);
}

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