TMDCUtil.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TMDCUtil.cxx,v 1.7 2010/03/31 09:58:59 liucy Exp $
//-----------------------------------------------------------------------------
// Filename : TMDCUtil.cc
// Section  : Tracking MDC
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : Utilities for MDC analyses and tracking.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
/* for erfc */
#if defined( __sparc )
#  if defined( __EXTENSIONS__ )
#    include <cmath>
#  else
#    define __EXTENSIONS__
#    include <cmath>
#    undef __EXTENSIONS__
#  endif
#elif defined( __GNUC__ )
#  if defined( _XOPEN_SOURCE )
#    include <cmath>
#  else
#    define _XOPEN_SOURCE
#    include <cmath>
#    undef _XOPEN_SOURCE
#  endif
#endif
 
#if defined( __SUNPRO_CC )
// for erfc and other functions (lgamma and cbrt
#  include <math.h>
#endif
 
// #include "panther/panther.h"
// #include MDC_H
// #include "MdcRecGeo/MdcRecGeo.h"
// #include "MdcGeomSvc/MdcGeomSvc.h"
#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TMLink.h"
 
const HepPoint3D ORIGIN = HepPoint3D( 0., 0., 0. );
 
// added by matsu
 
// CathodeSectorId
//---------------------------------------
// input    : wire id ( 0 - 189 )
// output   : cathode sector id ( 0 - 23 )
//      ( layer 0 : 0,1,...,7 )
//      ( layer 1 : 8,8.5,9,...,15,15.5)
//      ( layer 2 : 16,17,...,23)
 
float CathodeSectorId( unsigned id ) {
 
  unsigned layer = id / 64;
 
  if ( layer == 0 ) { return int( id / 8 ); }
 
  if ( layer == 1 )
  {
    if ( id >= 127 ) id -= 64;
    if ( ( id - 6 ) % 8 == 0 ) return ( id - 6 ) / 8 + 0.5;
    else return int( ( id + 1 ) / 8 );
  }
 
  if ( layer == 2 )
  {
    if ( id <= 129 ) id += 64;
    return int( ( id - 2 ) / 8 );
  }
 
  return 9999;
}
// end
 
void bitDisplay( unsigned val ) { bitDisplay( val, 31, 0 ); }
 
void bitDisplay( unsigned val, unsigned f, unsigned l ) {
  unsigned i;
  for ( i = 0; i < f - l; i++ )
  {
    if ( ( i % 8 ) == 0 ) std::cout << " ";
    std::cout << ( val >> ( f - i ) ) % 2;
  }
}
 
int intersection( const HepPoint3D& c1, double r1, const HepPoint3D& c2, double r2, double eps,
                  HepPoint3D& x1, HepPoint3D& x2 ) {
 
  double c0x  = c2.x() - c1.x();
  double c0y  = c2.y() - c1.y();
  double c0   = sqrt( c0x * c0x + c0y * c0y );
  double rr1  = abs( r1 );
  double rr2  = abs( r2 );
  double Radd = rr1 + rr2;
  double Rsub = abs( rr1 - rr2 );
 
  // no intersections
 
  if ( c0 > Radd + eps || c0 < 0.001 || c0 < Rsub - eps )
  {
    //-- debug
    // std::cout << "Int2Cir return 0 " << std::endl;
    //-- debug end
    return 0;
  }
 
  // single intersection
 
  else
  {
    if ( c0 > Radd - eps )
    {
      x1.setX( c1.x() + rr1 * c0x / c0 );
      x1.setY( c1.y() + rr1 * c0y / c0 );
      x2.setX( 0.0 );
      x2.setY( 0.0 );
      //--debug
      // std::cout << "Int2Cir return 1" << std::endl;
      //--debug end
      return 1;
    }
  }
 
  // two intersections
 
  double chg    = abs( r1 ) / r1;
  double cosPsi = ( c0 * c0 + rr1 * rr1 - rr2 * rr2 ) / ( 2. * c0 * rr1 );
  double sinPsi = -( chg / abs( chg ) ) * sqrt( 1.0 - cosPsi * cosPsi );
  x1.setX( c1.x() + ( rr1 / c0 ) * ( cosPsi * c0x - sinPsi * c0y ) );
  x1.setY( c1.y() + ( rr1 / c0 ) * ( cosPsi * c0y + sinPsi * c0x ) );
  x2.setX( c1.x() + ( rr1 / c0 ) * ( cosPsi * c0x + sinPsi * c0y ) );
  x2.setY( c1.y() + ( rr1 / c0 ) * ( cosPsi * c0y - sinPsi * c0x ) );
  //-- debug
  // std::cout << "Int2Cir return 2" << std::endl;
  //-- debug end
  return 2;
}
 
// from by jtanaka
double chisq2confLevel( int n, double chi2 ) {
#define SRTOPI 0.7978846
#define UPL 340.0
#define ROOT2I 0.70710678
 
  double prob = 0.0;
  double sum, term;
  int    m;
  int    i, k;
  double temp_i, temp_n;
  double srty;
 
  if ( ( n <= 0 ) || ( chi2 < 0.0 ) ) { return prob; }
  if ( n > 60 )
  {
    temp_n = (double)n;
    srty   = sqrt( chi2 ) - sqrt( temp_n - 0.5 );
    if ( srty < 12.0 )
    {
      prob = 0.5 * erfc( srty );
      return prob;
    }
    return prob;
  }
  if ( chi2 > UPL ) { return prob; }
  sum  = exp( -0.5 * chi2 );
  term = sum;
  m    = (int)floor( n / 2. );
 
  if ( 2 * m == n )
  {
    if ( m == 1 )
    {
      prob = sum;
      return prob;
    }
    else
    {
      for ( i = 2; i < m + 1; i++ )
      {
        temp_i = (double)i;
        term   = 0.5 * chi2 * term / ( temp_i - 1.0 );
        sum    = sum + term;
      }
      prob = sum;
      return prob;
    }
  }
  else
  {
    srty = sqrt( chi2 );
    prob = erfc( ROOT2I * srty );
    if ( n == 1 ) { return prob; }
    if ( n == 3 )
    {
      prob = SRTOPI * srty * sum + prob;
      return prob;
    }
    else
    {
      k = m - 1;
      for ( i = 1; i < k + 1; i++ )
      {
        temp_i = (double)i;
        term   = term * chi2 / ( 2.0 * temp_i + 1.0 );
        sum    = sum + term;
      }
      prob = SRTOPI * srty * sum + prob;
      return prob;
    }
  }
}