00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include "MdcTrkRecon/MdcTrack.h"
00017 #include <math.h>
00018 #include <stdlib.h>
00019 #include <assert.h>
00020 #include "MdcGeom/MdcLayer.h"
00021 #include "MdcData/MdcHitOnTrack.h"
00022 #include "MdcGeom/BesAngle.h"
00023 #include "MdcGeom/Constants.h"
00024 #include "CLHEP/Alist/AList.h"
00025 #include "TrkBase/TrkRep.h"
00026 #include "TrkFitter/TrkHelixRep.h"
00027 #include "TrkBase/TrkExchangePar.h"
00028 #include "TrkBase/TrkRecoTrk.h"
00029 #include "TrkFitter/TrkCircleMaker.h"
00030 #include "TrkBase/TrkContext.h"
00031
00032 #include "MdcData/MdcRecoHitOnTrack.h"
00033 #include "BField/BField.h"
00034 #include "MdcData/MdcHit.h"
00035 #include "MdcData/MdcHitOnTrack.h"
00036 #include "MdcRawEvent/MdcDigi.h"
00037 #include "Identifier/Identifier.h"
00038 #include "Identifier/MdcID.h"
00039
00040 MdcTrack::MdcTrack(TrkRecoTrk* aTrack){
00041
00042 _theTrack = aTrack;
00043 _firstLayer = _lastLayer = 0;
00044 _haveCurled = 0;
00045 }
00046
00047
00048 MdcTrack::MdcTrack(int nsupers, const TrkExchangePar& par, double chisq,
00049 TrkContext& context, double trackT0) {
00050
00051 TrkCircleMaker maker;
00052 _theTrack = maker.makeTrack(par, chisq, context, trackT0);
00053 _firstLayer = _lastLayer = 0;
00054 _haveCurled = 0;
00055 }
00056
00057
00058 MdcTrack::~MdcTrack() {
00059
00060 delete _theTrack;
00061 }
00062
00063 int MdcTrack::projectToR(double radius, BesAngle &phiIntersect, int lCurl)
00064 const {
00065
00066
00067
00068
00069
00070 const TrkFit* tkFit = track().fitResult();
00071 if (tkFit == 0) return -1;
00072 TrkExchangePar par = tkFit->helix(0.0);
00073 double d0 = par.d0();
00074 double omega = par.omega();
00075 double phi0 = par.phi0();
00076 double r2d2 = radius * radius - d0 * d0;
00077 if (r2d2 < 0) return -1;
00078 double rinv = 1./radius;
00079 double k2dinv = 1./(1 + omega * d0);
00080 if (!lCurl) {
00081 double arg = d0 * rinv + 0.5*omega * rinv * r2d2 * k2dinv;
00082 if (fabs(arg) > 1.0) return -1;
00083 phiIntersect.setRad( phi0 + asin(arg) );
00084 }
00085 else {
00086 double arg = -d0 * rinv - 0.5*omega * rinv * r2d2 * k2dinv;
00087 if (fabs(arg) > 1.0) return -1;
00088 phiIntersect.setRad( phi0 + Constants::pi + asin(arg) );
00089 }
00090
00091 return 0;
00092 }
00093
00094
00095 int MdcTrack::projectToR(double radius, BesAngle &phiIntersect,
00096 double &arcLength, int lCurl) const {
00097
00098
00099
00100
00101
00102
00103 const TrkFit* tkFit = track().fitResult();
00104 if (tkFit == 0) return -1;
00105 TrkExchangePar par = tkFit->helix(0.0);
00106 double d0 = par.d0();
00107 double omega = par.omega();
00108 double phi0 = par.phi0();
00109 double r2d2 = radius * radius - d0 * d0;
00110 if (r2d2 < 0) return -1;
00111 double rinv = 1./radius;
00112 double k2dinv = 1./(1 + omega * d0);
00113 if (k2dinv < 0.0) return -1;
00114 double arg;
00115
00116 if (!lCurl) {
00117 arg = d0 * rinv + 0.5*omega * rinv * r2d2 * k2dinv;
00118 if (fabs(arg) > 1.0) return -1;
00119 phiIntersect.setRad( phi0 + asin(arg) );
00120 }
00121 else {
00122 arg = -d0 * rinv - 0.5*omega * rinv * r2d2 * k2dinv;
00123 if (fabs(arg) > 1.0) return -1;
00124 phiIntersect.setRad( phi0 + Constants::pi + asin(arg));
00125 }
00126
00127 arg = 0.5*omega * sqrt(r2d2 * k2dinv);
00128 arcLength = 2. * asin(arg) / omega;
00129
00130 return 0;
00131 }
00132
00133
00134 bool
00135 MdcTrack::operator==(const MdcTrack& tk) const {
00136
00137
00138 return (track() == tk.track());
00139 }
00140
00141
00142 void
00143 MdcTrack::storeTrack(int trackId, RecMdcTrackCol* trackList, RecMdcHitCol* hitList,int tkStat){
00144
00145
00146
00147
00148 const TrkFit* fitresult = track().fitResult();
00149
00150 if (fitresult == 0) return;
00151
00152
00153 RecMdcTrack* recMdcTrack = new RecMdcTrack();
00154 const TrkHitList* aList = track().hits();
00155
00156 const BField& theField = track().bField();
00157 double Bz = theField.bFieldZ();
00158
00159 int hitId = 0;
00160 int nHits=0;
00161 int nSt=0;
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174 nHits = aList->nHit();
00175
00176
00177
00178
00179
00180 int q = fitresult->charge();
00181 double chisq = fitresult->chisq();
00182 int nDof = fitresult->nDof();
00183
00184 double fltLenPoca = 0.0;
00185 TrkExchangePar helix = fitresult->helix(fltLenPoca);
00186
00187 double phi0 = helix.phi0();
00188 double tanDip = helix.tanDip();
00189
00190 double d0 = helix.d0();
00191
00192
00193 HepPoint3D poca = fitresult->position(fltLenPoca);
00194
00195 double helixPar[5];
00196
00197 helixPar[0] = -d0;
00198
00199 double tphi0 = phi0 - Constants::pi/2.;
00200 if (tphi0 < 0. ) tphi0 += Constants::pi * 2.;
00201 helixPar[1] = tphi0;
00202
00203 double pxy = fitresult->pt();
00204 if (pxy == 0.) helixPar[2] = 9999.;
00205 else helixPar[2] = q/fabs(pxy);
00206 if(pxy>9999.) helixPar[2] = 0.00001;
00207
00208 helixPar[3] = helix.z0();
00209
00210 helixPar[4] = tanDip;
00211
00212
00213 HepSymMatrix mS(helix.params().num_row(),0);
00214 mS[0][0]=-1.;
00215 mS[1][1]=1.;
00216 mS[2][2]=-333.567/Bz;
00217 mS[3][3]=1.;
00218 mS[4][4]=1.;
00219 HepSymMatrix mVy = helix.covariance().similarity(mS);
00220 double errorMat[15];
00221 int k = 0;
00222 for (int ie = 0 ; ie < 5 ; ie ++){
00223 for (int je = ie ; je < 5 ; je ++) {
00224 errorMat[k] = mVy[ie][je];
00225 k++;
00226 }
00227 }
00228 double p,px,py,pz;
00229 px = pxy * (-sin(helixPar[1]));
00230 py = pxy * cos(helixPar[1]);
00231 pz = pxy * helixPar[4];
00232 p = sqrt(pxy*pxy + pz*pz);
00233
00234 double theta = acos(pz/p);
00235 double phi = atan2(py,px);
00236 recMdcTrack->setTrackId(trackId);
00237 recMdcTrack->setHelix(helixPar);
00238 recMdcTrack->setCharge(q);
00239 recMdcTrack->setPxy(pxy);
00240 recMdcTrack->setPx(px);
00241 recMdcTrack->setPy(py);
00242 recMdcTrack->setPz(pz);
00243 recMdcTrack->setP(p);
00244 recMdcTrack->setTheta(theta);
00245 recMdcTrack->setPhi(phi);
00246 recMdcTrack->setPoca(poca);
00247 recMdcTrack->setX(poca.x());
00248 recMdcTrack->setY(poca.y());
00249 recMdcTrack->setZ(poca.z());
00250 recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
00251 HepPoint3D pivot(0.,0.,0.);
00252 recMdcTrack->setPivot(pivot);
00253 recMdcTrack->setVX0(0.);
00254 recMdcTrack->setVY0(0.);
00255 recMdcTrack->setVZ0(0.);
00256 recMdcTrack->setError(errorMat);
00257 recMdcTrack->setError(mVy);
00258 recMdcTrack->setChi2(chisq);
00259 recMdcTrack->setNdof(nDof);
00260 recMdcTrack->setStat(tkStat);
00261 recMdcTrack->setNhits(nHits);
00262
00263 HitRefVec hitRefVec;
00264 vector<int> vecHits;
00265
00266 int maxLayerId = -1;
00267 int minLayerId = 43;
00268 double fiTerm = 999.;
00269 const MdcRecoHitOnTrack* fiTermHot = NULL;
00270 TrkHitList::hot_iterator hot = aList->begin();
00271 for (;hot!=aList->end();hot++){
00272 const MdcRecoHitOnTrack* recoHot;
00273 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
00274
00275 RecMdcHit* recMdcHit = new RecMdcHit;
00276
00277 recMdcHit->setId(hitId);
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289 double hotWireAmbig = recoHot->wireAmbig();
00290 double driftDist = fabs(recoHot->drift());
00291 double sigma = recoHot->hitRms();
00292 double doca = fabs(recoHot->dcaToWire());
00293 int flagLR=2;
00294 if ( hotWireAmbig == 1){
00295 flagLR = 0;
00296 doca *= -1.;
00297 }else if( hotWireAmbig == -1){
00298 flagLR = 1;
00299 }else if( hotWireAmbig == 0){
00300 flagLR = 2;
00301 }
00302 recMdcHit->setFlagLR(flagLR);
00303 recMdcHit->setDriftDistLeft((-1)*driftDist);
00304 recMdcHit->setErrDriftDistLeft(sigma);
00305 recMdcHit->setDriftDistRight(driftDist);
00306 recMdcHit->setErrDriftDistRight(sigma);
00307
00308 Identifier mdcId = recoHot->mdcHit()->digi()->identify();
00309 recMdcHit->setMdcId(mdcId);
00310
00311
00312
00313
00314 double res=999.,rese=999.;
00315 if (recoHot->resid(res,rese,false)){
00316 }else{}
00317 double deltaChi=0;
00318 recoHot->getFitStuff(deltaChi);
00319 recMdcHit->setChisqAdd( deltaChi * deltaChi );
00320
00321 recMdcHit->setTrkId(trackId);
00322
00323 recMdcHit->setDoca(doca);
00324
00325
00326 recMdcHit->setEntra(recoHot->entranceAngle());
00327
00328 HepPoint3D pos; Hep3Vector dir;
00329 double fltLen = recoHot->fltLen();
00330 recMdcHit->setFltLen(fltLen);
00331 recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
00332 recMdcHit->setZhit(pos.z());
00333 double tof = recoHot->getParentRep()->arrivalTime(fltLen);
00334 recMdcHit->setTdc(recoHot->mdcHit()->tdcIndex());
00335 recMdcHit->setAdc(recoHot->mdcHit()->adcIndex());
00336
00337 recMdcHit->setDriftT(recoHot->mdcHit()->driftTime(tof,pos.z()));
00338
00339 int layerId = recoHot->mdcHit()->layernumber();
00340 int wireId = recoHot->mdcHit()->wirenumber();
00341
00342
00343 if (layerId >= maxLayerId){
00344 maxLayerId = layerId;
00345 fiTermHot = recoHot;
00346 }
00347 if (layerId < minLayerId){
00348 minLayerId = layerId;
00349 }
00350
00351 if (recoHot->isActive()) {
00352 recMdcHit->setStat(1);
00353 }else{
00354 recMdcHit->setStat(0);
00355 }
00356
00357 if (recoHot->layer()->view()) {
00358 ++nSt;
00359 }
00360 hitList->push_back(recMdcHit);
00361 SmartRef<RecMdcHit> refHit(recMdcHit);
00362 hitRefVec.push_back(refHit);
00363 vecHits.push_back(mdcId.get_value());
00364 ++hitId;
00365 }
00366
00367 if (fiTermHot!=NULL){
00368 fiTerm=(1./sqrt(1.+tanDip*tanDip))*fiTermHot->fltLen()*helix.omega();
00369 }
00370 recMdcTrack->setFiTerm(fiTerm);
00371
00372 recMdcTrack->setNster(nSt);
00373 recMdcTrack->setFirstLayer(maxLayerId);
00374 recMdcTrack->setLastLayer(minLayerId);
00375 recMdcTrack->setVecHits(hitRefVec);
00376 trackList->push_back(recMdcTrack);
00377 }