---

SiDetLadder.cxx

Go to the documentation of this file.

// $Id: SiDetLadder.cxx,v 1.6 2000/11/30 01:42:34 jrb Exp $
// 
//  Original author: Sawyer Gillespie
//                   hgillesp@u.washington.edu
//
//  Clipping code, etc. provided by Marios Nikolaou, SLAC
//                   marios@slac.stanford.edu
//

#include "instrument/SiDetLadder.h"
#include <cfloat>
#include <cassert>
#include <algorithm>
#include "dom/DOM_Element.hpp"  // or maybe DOM_Document??
#include "xml/Dom.h"

// local namespace
namespace {

  // utility function declarations
#ifdef _MSC_VER
  inline double copysign(double a, double b){return _copysign(a,b);}
#endif


#define EPS 1e-8 // 0.1nm. Have to find a better way. This because of 
  // floating point errors.

  // ------------- 2D Liang-Barsky clipping algorithm ------------------------
  // Subroutine used by clipping algoritm. (2D Liang-Barsky algorithm)
  static inline bool clipT(double p, double q, double *tE, double *tL)
    {
      if (p < 0.0) {
        double t = q / p;
        if (t > *tL) return false;
        else  if (t > *tE)    *tE = t;
      }
      else  if (p > 0.0) {
        double t = q / p;
        if (t < *tE) return false;
        else  if (t < *tL)    *tL = t;
      }
      else       /* p == 0 */
        if (q < 0.0) return false;
            
      return true;
    }
    
  // Clip a line to a specified rectangle.
  // ..first & second point
  static bool clipDie(double *x0, double *y0, double *x1, double *y1, 
                       double xmin, double xmax, double ymin, double ymax)
    {
      double tE = 0.0;
      double tL = 1.0;
      double dX = *x1 - *x0;
      if (clipT(-dX, *x0 - xmin, &tE, &tL) && clipT(dX, xmax - *x0, &tE, &tL))
      {
        double dY = *y1 - *y0;
        if (clipT(-dY, *y0 - ymin, &tE, &tL) && 
            clipT(dY, ymax - *y0, &tE, &tL))   
        {
          if (tL < 1.0)  {
            *x1 = *x0 + (tL * dX);
            *y1 = *y0 + (tL * dY);
          }
          if (tE > 0.0)  {
            *x0 += tE * dX;
            *y0 += tE * dY;
          }
        }
        else return false;
      }
      else return false;
        
      // patch for floating point errors
      if (fabs(*x0-xmin) < 10*DBL_EPSILON) *x0=xmin;
      if (fabs(*x0-xmax) < 10*DBL_EPSILON) *x0=xmax;
      if (fabs(*x1-xmin) < 10*DBL_EPSILON) *x1=xmin;
      if (fabs(*x1-xmax) < 10*DBL_EPSILON) *x1=xmax;
      if (fabs(*y0-ymin) < 10*DBL_EPSILON) *y0=ymin;
      if (fabs(*y0-ymax) < 10*DBL_EPSILON) *y0=ymax;
      if (fabs(*y1-ymin) < 10*DBL_EPSILON) *y1=ymin;
      if (fabs(*y1-ymax) < 10*DBL_EPSILON) *y1=ymax;
        
      // Should not happen!
      assert((*x0>=xmin) && (*x0<=xmax));
      assert((*x1>=xmin) && (*x1<=xmax));
      assert((*y0>=ymin) && (*y0<=ymax));
      assert((*y1>=ymin) && (*y1<=ymax));
        
      return true;
    }
}

SiDetLadder::SiDetLadder(const DOM_Element& spec)
{
  // get all the parameters from the XML spec.
  m_diegap = atof(xml::Dom::transToChar(spec.getAttribute("diegap")));
  m_dieheight = atof(xml::Dom::transToChar(spec.getAttribute("dieheight")));
  m_guardring = atof(xml::Dom::transToChar(spec.getAttribute("guardring")));
  m_ndies = atof(xml::Dom::transToChar(spec.getAttribute("ndies")));
  m_strippitch = atof(xml::Dom::transToChar(spec.getAttribute("strippitch")));
  m_width = atof(xml::Dom::transToChar(spec.getAttribute("diewidth")));
  m_stripzero = 0;
}

SiDetLadder::~SiDetLadder ()
{
}

void    SiDetLadder::score(const Point& o, const Point& p, float eloss,
                           SiStripList& l, int stripoffset) const
{
  Vector dir = (p - o);
  double xDir = dir.x();
    
  double in = o.x();                                // entry point -- x
  double ex = p.x();                                // exit point -- x
 // enter/exit strips (if valid)
  unsigned ins = stripId(in), exs = stripId(ex); 
    
  // if these assertions fire, it probably is because the clipping algorithm
  // may return an intersection, which lies just outside the active area.
  // This due to floting point errors.
  // Has temporarily been fixed by giving a slightly smaller clipping region.
  //    assert(ins != SiStrip::undef_strip());
  //    assert(exs != SiStrip::undef_strip());
    
  if (ins == exs) {
    if (ins != SiStrip::undef_strip()) l.push_back(SiStrip(ins + stripoffset,
                                                            eloss));
  }
  else {
    double sx = localX (ins);
    double dx = m_strippitch/2. - (in-sx)*copysign(1.,xDir);
    l.push_back(SiStrip(ins + stripoffset, eloss*dx/fabs(xDir)));// entry strip
    sx = localX (exs);
    dx = m_strippitch/2. + (ex-sx)*copysign(1.,xDir);
    l.push_back(SiStrip(exs + stripoffset, eloss*dx/fabs(xDir))); // exit strip
        
    // add energy to strips between entry and exit
    short sinc = (ins < exs) ? 1 : -1;
    if (fabs(xDir) > m_strippitch)
      for (unsigned sid = ins + sinc; sid != exs; sid += sinc)   
        l.push_back(SiStrip(sid + stripoffset, eloss*m_strippitch/fabs(xDir)));
    else
      for (unsigned sid = ins + sinc; sid != exs; sid += sinc) { 
        //THB fprintf(std::cerr,"Should not have printed this!\n");
        l.push_back(SiStrip(sid + stripoffset, eloss));
      }
  }
}

SiStripList SiDetLadder::score (const Point& entry, const Point& exit
                                 , float eloss, int stripoffset) const
{
  SiStripList retval; // return value

  // call the other method
  score (entry, exit, eloss, retval, stripoffset);

  // this calls SiStripList copy constructor (deep copy)
  return retval;
}
double   SiDetLadder::insideActiveArea (double x) const {
  double  xin = 0., yin = 0.;
    
  // compute the distance to the nearest x-edge
  xin = x - m_guardring;
  xin = std::min (xin, m_width - m_guardring*2 - xin);
    
  return xin;
}


double   SiDetLadder::insideActiveArea (double x, double y) const {
  double  xin = 0., yin = 0.;
    
  // compute the distance to the nearest x-edge
  xin = x - m_guardring;
  xin = std::min (xin, m_width - m_guardring*2 - xin);
  if (xin < 0) return xin;

  // compute the distance to the nearest y edge
  double diefull = m_dieheight + m_diegap;
  double die = floor (y / diefull);
  yin = y - diefull * die - m_guardring;
  yin = std::min(yin, diefull - m_diegap - m_guardring*2 - yin);
  if (yin < 0) return yin;

  // finally, return the minimum distance to the edge
  return std::min(xin, yin);
}

---