00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #if defined(__sparc)
00015 # if defined(__EXTENSIONS__)
00016 # include <cmath>
00017 # else
00018 # define __EXTENSIONS__
00019 # include <cmath>
00020 # undef __EXTENSIONS__
00021 # endif
00022 #elif defined(__GNUC__)
00023 # if defined(_XOPEN_SOURCE)
00024 # include <cmath>
00025 # else
00026 # define _XOPEN_SOURCE
00027 # include <cmath>
00028 # undef _XOPEN_SOURCE
00029 # endif
00030 #endif
00031
00032 #if defined(__SUNPRO_CC)
00033
00034 # include <math.h>
00035 #endif
00036
00037
00038
00039
00040
00041 #include "TrkReco/TMDCUtil.h"
00042 #include "TrkReco/TMDCWire.h"
00043 #include "TrkReco/TMDCWireHit.h"
00044 #include "TrkReco/TMDCWireHitMC.h"
00045 #include "TrkReco/TMLink.h"
00046
00047 const HepPoint3D ORIGIN = HepPoint3D(0., 0., 0.);
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 float
00060 CathodeSectorId(unsigned id) {
00061
00062 unsigned layer = id/64;
00063
00064 if ( layer == 0 ){
00065 return int(id/8);
00066 }
00067
00068 if( layer == 1 ){
00069 if( id >= 127 ) id -= 64;
00070 if( (id-6)%8 == 0 ) return (id-6)/8 + 0.5;
00071 else return int((id+1)/8);
00072 }
00073
00074 if ( layer == 2 ) {
00075 if( id <= 129 ) id += 64;
00076 return int((id-2)/8);
00077 }
00078
00079 return 9999;
00080
00081 }
00082
00083
00084 void
00085 bitDisplay(unsigned val) {
00086 bitDisplay(val, 31, 0);
00087 }
00088
00089 void
00090 bitDisplay(unsigned val, unsigned f, unsigned l) {
00091 unsigned i;
00092 for (i = 0; i < f - l; i++) {
00093 if ((i % 8) == 0) std::cout << " ";
00094 std::cout << (val >> (f - i)) % 2;
00095 }
00096 }
00097
00098 int
00099 intersection(const HepPoint3D & c1,
00100 double r1,
00101 const HepPoint3D & c2,
00102 double r2,
00103 double eps,
00104 HepPoint3D & x1,
00105 HepPoint3D & x2) {
00106
00107 double c0x = c2.x() - c1.x() ;
00108 double c0y = c2.y() - c1.y() ;
00109 double c0 = sqrt ( c0x*c0x + c0y*c0y ) ;
00110 double rr1 = abs(r1) ;
00111 double rr2 = abs(r2) ;
00112 double Radd = rr1 + rr2 ;
00113 double Rsub = abs( rr1 - rr2 ) ;
00114
00115
00116
00117 if ( c0 > Radd + eps || c0 < 0.001 || c0 < Rsub - eps ) {
00118
00119
00120
00121 return 0 ;
00122 }
00123
00124
00125
00126 else {
00127 if ( c0 > Radd - eps ) {
00128 x1.setX(c1.x() + rr1*c0x/c0);
00129 x1.setY(c1.y() + rr1*c0y/c0);
00130 x2.setX(0.0);
00131 x2.setY(0.0);
00132
00133
00134
00135 return 1 ;
00136 }
00137 }
00138
00139
00140
00141 double chg = abs(r1) / r1 ;
00142 double cosPsi = ( c0*c0 + rr1*rr1 - rr2*rr2 ) / (2.*c0*rr1 ) ;
00143 double sinPsi = - ( chg/abs(chg) ) * sqrt(1.0 - cosPsi*cosPsi) ;
00144 x1.setX(c1.x() + ( rr1/c0 )*( cosPsi*c0x - sinPsi*c0y ));
00145 x1.setY(c1.y() + ( rr1/c0 )*( cosPsi*c0y + sinPsi*c0x ));
00146 x2.setX(c1.x() + ( rr1/c0 )*( cosPsi*c0x + sinPsi*c0y ));
00147 x2.setY(c1.y() + ( rr1/c0 )*( cosPsi*c0y - sinPsi*c0x ));
00148
00149
00150
00151 return 2 ;
00152
00153 }
00154
00155
00156 double chisq2confLevel(int n,double chi2){
00157 #define SRTOPI 0.7978846
00158 #define UPL 340.0
00159 #define ROOT2I 0.70710678
00160
00161 double prob = 0.0;
00162 double sum,term;
00163 int m;
00164 int i,k;
00165 double temp_i,temp_n;
00166 double srty;
00167
00168 if((n <= 0)||(chi2 < 0.0)){
00169 return prob;
00170 }
00171 if(n > 60){
00172 temp_n = (double)n;
00173 srty = sqrt(chi2) - sqrt(temp_n-0.5);
00174 if (srty < 12.0){
00175 prob = 0.5*erfc(srty);
00176 return prob;
00177 }
00178 return prob;
00179 }
00180 if(chi2 > UPL){
00181 return prob;
00182 }
00183 sum = exp( -0.5 * chi2 );
00184 term = sum;
00185 m = (int)floor(n/2.);
00186
00187 if( 2*m == n ){
00188 if( m == 1 ){
00189 prob = sum;
00190 return prob;
00191 }else{
00192 for(i=2;i<m+1;i++){
00193 temp_i = (double)i;
00194 term = 0.5*chi2*term/(temp_i-1.0);
00195 sum = sum + term;
00196
00197 }
00198 prob = sum;
00199 return prob;
00200 }
00201 }else{
00202 srty = sqrt(chi2);
00203 prob = erfc(ROOT2I*srty);
00204 if(n == 1){
00205 return prob;
00206 }
00207 if(n == 3){
00208 prob = SRTOPI*srty*sum + prob;
00209 return prob;
00210 }else{
00211 k = m - 1;
00212 for(i=1;i<k+1;i++){
00213 temp_i = (double)i;
00214 term = term*chi2/(2.0*temp_i + 1.0);
00215 sum = sum + term;
00216 }
00217 prob = SRTOPI*srty*sum + prob;
00218 return prob;
00219 }
00220 }
00221 }