00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "GaudiKernel/MsgStream.h"
00026 #include "GaudiKernel/AlgFactory.h"
00027 #include "GaudiKernel/ISvcLocator.h"
00028 #include "GaudiKernel/SmartDataPtr.h"
00029 #include "GaudiKernel/IDataProviderSvc.h"
00030 #include "GaudiKernel/PropertyMgr.h"
00031 #include "EventModel/Event.h"
00032 #include "MdcRawEvent/MdcDigi.h"
00033 #include "EmcRawEvent/EmcDigi.h"
00034 #include "McTruth/McParticle.h"
00035
00036 #include "MdcGeomSvc/IMdcGeomSvc.h"
00037 #include "MdcGeomSvc/MdcGeoWire.h"
00038 #include "MdcGeomSvc/MdcGeoLayer.h"
00039 #include <vector>
00040 #include <iostream>
00041 #include <cstdio>
00042
00043
00044 #include "EsTimeAlg/Emc_helix.h"
00045 #include "CLHEP/Vector/ThreeVector.h"
00046 #include "MdcFastTrkAlg/FTFinder.h"
00047
00048 #include "CLHEP/Geometry/Point3D.h"
00049
00050 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00051 typedef HepGeom::Point3D<double> HepPoint3D;
00052 #endif
00053
00054 using CLHEP::HepVector;
00055 using CLHEP::Hep3Vector;
00056
00057
00058 Emc_helix::Emc_helix(void) {
00059 piby1=3.141593;
00060 pi2=2.0*piby1;
00061 _debug=0;
00062 _pathl_cut= 1000.0;
00063
00064 }
00065
00066
00067 int Emc_helix::Emc_Get(double rEmc, int id, double fzisan[]) {
00068
00069 double xc,yc;
00070 double rho;
00071 double xx[2],yy[2];
00072 double emcx,emcy;
00073 double nx,ny;
00074 double cosx,sinx;
00075 double aa,ca,cb,cc,detm;
00076
00077
00078
00079 int Id_cdfztrk = id; R_emc = rEmc;
00080
00081 if( Id_cdfztrk<0 ) {
00082 if (_debug) std::cout << " Emc_helix *** wrong id_cdfztrk ="
00083 << Id_cdfztrk << std::endl;
00084 return(-1);
00085 }
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106 const HepPoint3D pivot1(0.,0.,0.);
00107 HepVector a(5,0);
00108 a[0] = fzisan[0];
00109 a[1] = fzisan[1];
00110 a[2] = fzisan[2];
00111 a[3] = fzisan[3];
00112 a[4] = fzisan[4];
00113
00114
00115 if (abs(a[3])>50.0 || abs(a[4])>500.0) return (-5);
00116
00117
00118 Helix helix1(pivot1,a);
00119 Dr = helix1.a()[0];
00120 Phi0 = helix1.a()[1];
00121 Kappa = helix1.a()[2];
00122 Dz = helix1.a()[3];
00123 Tanl = helix1.a()[4];
00124
00125
00126
00127
00128
00129
00130 rho = helix1.radius();
00131
00132 HepPoint3D xyz = helix1.center();
00133 xc = xyz(0);
00134 yc = xyz(1);
00135
00136 Hep3Vector vxyz = helix1.direction();
00137 nx = vxyz(0);
00138 ny = vxyz(1);
00139
00140
00141 if( xc==0.0 && yc==0.0 ) return(-3);
00142
00143
00144 ca = xc*xc + yc*yc ;
00145 aa = 0.5 * ( ca - rho*rho + R_emc*R_emc );
00146
00147 if ( xc != 0.0 ) {
00148 cb = aa * yc;
00149 cc = aa*aa - R_emc*R_emc*xc*xc;
00150
00151 detm = cb*cb - ca*cc;
00152 if ( detm > 0.0 ) {
00153 yy[0] = ( cb + sqrt(detm) )/ ca;
00154 xx[0] = ( aa - yy[0]*yc )/xc;
00155 yy[1] = ( cb - sqrt(detm) )/ ca;
00156 xx[1] = ( aa - yy[1]*yc )/xc;
00157 } else return(-1);
00158 }
00159 else{
00160 cb = aa * xc;
00161 cc = aa*aa - R_emc*R_emc*yc*yc;
00162
00163 detm = cb*cb - ca*cc;
00164 if ( detm > 0.0 ) {
00165 xx[0] = ( cb + sqrt(detm) )/ca;
00166 yy[0] = ( aa - xx[0]*xc )/yc;
00167 xx[1] = ( cb - sqrt(detm) )/ca;
00168 yy[1] = ( aa - xx[1]*xc )/yc;
00169 } else return(-2);
00170 }
00171
00172
00173
00174
00175 if( xx[0]*nx + yy[0]*ny > 0.0 ) { emcx = xx[0]; emcy = yy[0]; }
00176 else { emcx = xx[1]; emcy = yy[1]; }
00177
00178 double fi = atan2(emcy ,emcx );
00179
00180 if( fi < 0.0 ) fi=pi2+fi;
00181 Fi_emc = fi;
00182
00183
00184 const HepPoint3D pivot2( emcx, emcy, 0.0);
00185 helix1.pivot(pivot2);
00186
00187
00188 Phi1=helix1.a()[1];
00189 Z_emc=helix1.a()[3];
00190
00191
00192 Hep3Vector intersec(emcx, emcy, Z_emc);
00193 theta_emc = intersec.getTheta();
00194 phi_emc = intersec.getPhi();
00195
00196 return(1);
00197
00198 }
00199
00200