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    : [email protected]
//-----------------------------------------------------------------------------
// 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;
    }
  }
}