00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 using namespace std;
00013
00014 #include <iostream>
00015 #include <vector>
00016 #include <CLHEP/Vector/ThreeVector.h>
00017 #include <CLHEP/Geometry/Point3D.h>
00018
00019 #include "MucRecEvent/RecMucTrack.h"
00020 #include "MucRecEvent/MucTrackParameter.h"
00021 #include "MucRecEvent/MucRec2DRoad.h"
00022 #include "MucRecEvent/MucRec3DRoad.h"
00023 #include "MucGeomSvc/MucConstant.h"
00024
00025
00026 RecMucTrack::RecMucTrack()
00027 :
00028 m_ExtTrackID(-1),
00029 m_MdcPos(0.0, 0.0, 0.0),
00030 m_MdcMomentum(0.0, 0.0, 0.0),
00031 m_MucPos(0.0, 0.0, 0.0),
00032 m_MucPosSigma(0.0, 0.0, 0.0),
00033 m_MucMomentum(0.0, 0.0, 0.0),
00034 m_CurrentPos(0.0, 0.0, 0.0),
00035 m_CurrentDir(0.0, 0.0, 0.0),
00036 m_CurrentInsct(0.0, 0.0, 0.0),
00037 m_Good3DLine(0),
00038 m_pHits(0),
00039 m_pExpectedHits(0),
00040 m_Intersections(0),
00041 m_Directions(0)
00042 {
00043
00044 for(int igap = 0; igap < 9; igap++){
00045 m_IntersectionInner[igap].set(-9999,-9999,-9999);
00046 m_IntersectionOuter[igap].set(-9999,-9999,-9999);
00047 }
00048 m_id = 0;
00049 m_status = -1;
00050 m_type = -1;
00051
00052 m_numHits = 0;
00053 m_startPart = -1;
00054 m_endPart = -1;
00055 m_brLastLayer = -1;
00056 m_ecLastLayer = -1;
00057 m_brFirstLayer = -1;
00058 m_ecFirstLayer = -1;
00059 m_ecPart = -1;
00060 m_numLayers = 0;
00061 m_maxHitsInLayer = 0;
00062 m_depth = -99;
00063 m_dof = 0;
00064 m_chi2 = 0.0;
00065 m_rms = 0.0;
00066 m_deltaPhi = 0.0;
00067
00068 m_xPosSigma = 0.0;
00069 m_yPosSigma = 0.0;
00070 m_zPosSigma = 0.0;
00071
00072 m_changeUnit = false;
00073 m_recmode = 0;
00074
00075 m_kalrechi2 =0. ;
00076 m_kaldof =0;
00077 m_kaldepth = -99;
00078 m_kalbrLastLayer = -1;
00079 m_kalecLastLayer = -1;
00080 }
00081
00082
00083 RecMucTrack&
00084 RecMucTrack::operator=(const RecMucTrack& orig)
00085 {
00086
00087 if ( this != &orig ) {
00088
00089 m_ExtTrackID = orig.m_ExtTrackID;
00090 m_MdcPos = orig.m_MdcPos;
00091 m_MdcMomentum = orig.m_MdcMomentum;
00092 m_MucPos = orig.m_MucPos;
00093 m_MucPosSigma = orig.m_MucPosSigma;
00094 m_MucMomentum = orig.m_MucMomentum;
00095 m_CurrentPos = orig.m_CurrentPos;
00096 m_CurrentDir = orig.m_CurrentDir;
00097 m_CurrentInsct = orig.m_CurrentInsct;
00098 m_id = orig.m_id;
00099 m_status = orig.m_status;
00100 m_type = orig.m_type;
00101 m_numHits = orig.m_numHits;
00102 m_startPart = orig.m_startPart;
00103 m_endPart = orig.m_endPart;
00104 m_brLastLayer = orig.m_brLastLayer;
00105 m_ecLastLayer = orig.m_ecLastLayer;
00106 m_numLayers = orig.m_numLayers;
00107 m_maxHitsInLayer= orig.m_maxHitsInLayer;
00108 m_depth = orig.m_depth;
00109 m_dof = orig.m_dof;
00110 m_chi2 = orig.m_chi2;
00111 m_rms = orig.m_rms;
00112 m_deltaPhi = orig.m_deltaPhi;
00113 m_pHits = orig.m_pHits;
00114 m_pExpectedHits = orig.m_pExpectedHits;
00115 m_Intersections = orig.m_Intersections;
00116 m_Directions = orig.m_Directions;
00117 m_xPosSigma = orig.m_xPosSigma;
00118 m_yPosSigma = orig.m_yPosSigma;
00119 m_zPosSigma = orig.m_zPosSigma;
00120 m_changeUnit = orig.m_changeUnit;
00121 m_recmode = orig.m_recmode;
00122
00123 m_kalrechi2 = orig.m_kalrechi2;
00124 m_kaldof = orig.m_kaldof;
00125 m_kaldepth =orig.m_kaldepth;
00126 m_kalbrLastLayer=orig.m_kalbrLastLayer;
00127 m_kalecLastLayer=orig.m_kalecLastLayer;
00128 }
00129
00130 return *this;
00131 }
00132
00133
00134 RecMucTrack::RecMucTrack(const RecMucTrack& source)
00135 :
00136 m_ExtTrackID (source.m_ExtTrackID),
00137 m_MdcPos (source.m_MdcPos),
00138 m_MdcMomentum (source.m_MdcMomentum),
00139 m_MucPos (source.m_MucPos),
00140 m_MucPosSigma (source.m_MucPosSigma),
00141 m_MucMomentum (source.m_MucMomentum),
00142 m_CurrentPos (source.m_CurrentPos),
00143 m_CurrentDir (source.m_CurrentDir),
00144 m_CurrentInsct (source.m_CurrentInsct),
00145 m_pHits (source.m_pHits),
00146 m_pExpectedHits(source.m_pExpectedHits),
00147 m_Intersections(source.m_Intersections),
00148 m_Directions (source.m_Directions)
00149 {
00150 m_id = source.m_id;
00151 m_status = source.m_status;
00152 m_type = source.m_type;
00153 m_numHits = source.m_numHits;
00154 m_startPart = source.m_startPart;
00155 m_endPart = source.m_endPart;
00156 m_brLastLayer = source.m_brLastLayer;
00157 m_ecLastLayer = source.m_ecLastLayer;
00158 m_numLayers = source.m_numLayers;
00159 m_maxHitsInLayer= source.m_maxHitsInLayer;
00160 m_depth = source.m_depth;
00161 m_dof = source.m_dof;
00162 m_chi2 = source.m_chi2;
00163 m_rms = source.m_rms;
00164 m_deltaPhi = source.m_deltaPhi;
00165 m_xPosSigma = source.m_xPosSigma;
00166 m_yPosSigma = source.m_yPosSigma;
00167 m_zPosSigma = source.m_zPosSigma;
00168 m_changeUnit = source.m_changeUnit;
00169 m_recmode = source.m_recmode;
00170
00171 m_kalrechi2 = source.m_kalrechi2;
00172 m_kaldof = source.m_kaldof;
00173 m_kaldepth = source.m_kaldepth;
00174 m_kalbrLastLayer= source.m_kalbrLastLayer;
00175 m_kalecLastLayer= source.m_kalecLastLayer;
00176 }
00177
00178 RecMucTrack::RecMucTrack(const DstMucTrack& dstTrack)
00179 :DstMucTrack(dstTrack)
00180 {
00181
00182 SetDefault();
00183
00184 }
00185
00186 RecMucTrack& RecMucTrack::operator=(const DstMucTrack& dstTrack)
00187 {
00188 SetDefault();
00189 DstMucTrack::operator=(dstTrack);
00190 return *this;
00191 }
00192
00193 void RecMucTrack::SetDefault()
00194 {
00195 m_brFirstLayer = -99;
00196 m_ecFirstLayer = -99;
00197 m_Good3DPart = -99;
00198 }
00199
00200
00201 RecMucTrack::~RecMucTrack()
00202 {
00203 for(int i = 0 ; i < m_pExpectedHits.size(); i++)
00204 delete m_pExpectedHits[i];
00205 }
00206
00207
00208 void
00209 RecMucTrack::setTrackId(const int trackId)
00210 {
00211 if(trackId >= 0) {
00212 m_trackId = trackId;
00213 }
00214 }
00215
00216
00217
00218 void
00219 RecMucTrack::SetMdcPos(const float x,
00220 const float y,
00221 const float z)
00222 {
00223 m_MdcPos.set(x*10, y*10, z*10);
00224 }
00225
00226 void
00227 RecMucTrack::SetMdcMomentum(const float px,
00228 const float py,
00229 const float pz)
00230 {
00231 m_MdcMomentum.set(px*1000, py*1000, pz*1000);
00232 }
00233
00234 void
00235 RecMucTrack::SetMucPos(const float x,
00236 const float y,
00237 const float z)
00238 {
00239 m_MucPos.set(x*10, y*10, z*10);
00240 m_xPos = x*10;
00241 m_yPos = y*10;
00242 m_zPos = z*10;
00243 }
00244
00245 void
00246 RecMucTrack::SetMucPosSigma(const float x,
00247 const float y,
00248 const float z)
00249 {
00250 m_MucPosSigma.set(x*10, y*10, z*10);
00251 m_xPosSigma = x*10;
00252 m_yPosSigma = y*10;
00253 m_zPosSigma = z*10;
00254 }
00255
00256 void
00257 RecMucTrack::SetMucMomentum(const float px,
00258 const float py,
00259 const float pz)
00260 {
00261 m_MucMomentum.set(px*1000, py*1000, pz*1000);
00262 m_px = px*1000;
00263 m_py = py*1000;
00264 m_pz = pz*1000;
00265 }
00266
00267 void
00268 RecMucTrack::SetExtMucPos(const float x,
00269 const float y,
00270 const float z)
00271 {
00272 m_ExtMucPos.set(x*10, y*10, z*10);
00273 }
00274
00275 void
00276 RecMucTrack::SetExtMucMomentum(const float px,
00277 const float py,
00278 const float pz)
00279 {
00280 m_ExtMucMomentum.set(px*1000, py*1000, pz*1000);
00281 }
00282
00283 void
00284 RecMucTrack::SetCurrentPos(const float x,
00285 const float y,
00286 const float z)
00287 {
00288 m_CurrentPos.set(x*10, y*10, z*10);
00289 }
00290
00291 void
00292 RecMucTrack::SetCurrentDir(const float x,
00293 const float y,
00294 const float z)
00295 {
00296 m_CurrentDir.set(x*1000, y*1000, z*1000);
00297 }
00298
00299 void
00300 RecMucTrack::SetCurrentInsct(const float x,
00301 const float y,
00302 const float z)
00303 {
00304 m_CurrentInsct.set(x, y, z);
00305 }
00306
00307
00308
00309
00310
00311
00312
00313 void
00314 RecMucTrack::SetExtTrack(RecExtTrack *extTrack)
00315 {
00316 m_ExtTrack = extTrack;
00317 if (m_ExtTrack) m_ExtTrackID = extTrack->GetTrackId();
00318 }
00319
00320
00321
00322 void
00323 RecMucTrack::setVecHits(vector<MucRecHit*>& pHits)
00324 {
00325 for(int i = 0 ; i < pHits.size(); i++)
00326 m_vecHits.push_back((pHits[i]->GetID()).get_value());
00327 }
00328
00329 void
00330 RecMucTrack::setExpHits(vector<MucRecHit*>& pHits)
00331 {
00332 for(int i = 0 ; i < pHits.size(); i++)
00333 m_expHits.push_back((pHits[i]->GetID()).get_value());
00334 }
00335
00336 void
00337 RecMucTrack::GetMdcExtTrack(Hep3Vector mdcStartPos, Hep3Vector mdcStartMomentum, int charge,
00338 Hep3Vector &mucStartPos, Hep3Vector &mucStartMomentum)
00339 {
00340
00341
00342 mdcStartPos*=10;
00343 mdcStartMomentum*=1000;
00345
00346
00347 Hep3Vector pt = mdcStartMomentum;
00348 pt.setZ(0);
00349
00350 double radius = (pt.mag() * 1e6 / kvC * 1e3) / (fabs(charge * kMagnetField)) ;
00351
00352 double deltaPhi = -1.0 * (charge / abs(charge)) * kDeltaPhi;
00353 double deltaZ = (mdcStartMomentum.z() * 1e6 / kvC * 1e3) * kDeltaPhi / (abs(charge) * kMagnetField);
00354
00355
00356 mucStartPos = mdcStartPos;
00357 mucStartMomentum = mdcStartMomentum;
00358 double phi;
00359 int iter = 0;
00360 do {
00361 phi = mucStartMomentum.getPhi() + deltaPhi;
00362 mucStartPos.setX(mucStartPos.x() + radius * kDeltaPhi * cos(phi));
00363 mucStartPos.setY(mucStartPos.y() + radius * kDeltaPhi * sin(phi));
00364 mucStartPos.setZ(mucStartPos.z() + deltaZ);
00365
00366 mucStartMomentum.setPhi(mucStartMomentum.phi() + deltaPhi);
00367 iter++;
00368
00369 }
00370 while(IsInsideSuperConductor(mucStartPos) && iter < kMdcExtIterMax);
00371
00372
00373 mucStartPos/=10;
00374 mucStartMomentum/=1000;
00376 }
00377
00378 bool
00379 RecMucTrack::IsInsideSuperConductor(Hep3Vector pos)
00380 {
00381 if(pos.mag() < kSuperConductorR && fabs(pos.z()) < kSuperConductorZ) {
00382 return true;
00383 }
00384 else {
00385 return false;
00386 }
00387 }
00388
00389
00390
00391 void
00392 RecMucTrack::AttachHit(MucRecHit* hit)
00393 {
00394
00395
00396 if (!hit) {
00397 cout << "RecMucTrack::AttachHit-E1 null hit pointer!" << endl;
00398 return ;
00399 }
00400
00401 int part = hit->Part();
00402 int gap = hit->Gap();
00403 if ( (gap < 0) || (gap >= (int)MucID::getGapNum(part)) ) {
00404
00405 cout << "Muc2DRoad::AttachHit(MucRecHit*), bad gap number = " << gap
00406 << endl;
00407 return;
00408 }
00409
00410
00411 m_pHits.push_back(hit);
00412
00413
00414
00415
00416
00417
00418 }
00419
00420
00421 void
00422 RecMucTrack::Project(const int& part, const int& gap,
00423 float& x, float& y, float& z,
00424 int& seg)
00425 {
00426 seg = -1;
00427 x = 0.0; y = 0.0; z = 0.0;
00428
00429 if ( (gap < 0) || (gap >= (int)MucID::getGapNum(part)) ) {
00430 cout << "Muc2DRoad::Project(), invalid gap = " << gap << endl;
00431 return;
00432 }
00433
00434 HepPoint3D pos = m_CurrentPos;
00435 Hep3Vector dir = m_CurrentDir;
00436
00437 vector<HepPoint3D> insctCon = MucGeoGeneral::Instance()->FindIntersections(part, gap, pos, dir);
00438 if( insctCon.size() == 0 ) {
00439 return;
00440 }
00441
00442 HepPoint3D intersection = insctCon[0];
00443
00444 x = intersection.x();
00445 y = intersection.y();
00446 z = intersection.z();
00447
00448
00449
00450 float phi;
00451 if( (x*x+y*y+z*z) < kMinor) {
00452 return;
00453 }
00454 else {
00455 if(part == 1) {
00456 phi = acos(x/sqrt(x*x+y*y));
00457 if(y < 0) phi = 2.0*kPi - phi;
00458 phi = fmod((phi + kPi/8.0), 2.0*kPi);
00459 seg = int(phi/(kPi/4.0));
00460 }
00461 else {
00462 if(x >= 0.0) {
00463 if(y >= 0.0) {
00464 seg = 0;
00465 }
00466 else {
00467 seg = 3;
00468 }
00469 }
00470 else {
00471 if(y >= 0.0) {
00472 seg = 1;
00473 }
00474 else {
00475 seg = 2;
00476 }
00477 }
00478 }
00479 }
00480 }
00481
00482 float
00483 RecMucTrack::GetHitDistance(const MucRecHit* hit)
00484 {
00485 if (!hit) {
00486 cout << "RecMucTrack:GetHitDistance-E1 null hit pointer!" << endl;
00487 return kInfinity;
00488 }
00489
00490 int part = hit->Part();
00491 int gap = hit->Gap();
00492 if ( (gap < 0) || (gap >= (int)MucID::getGapNum(part)) ) {
00493 cout << "RecMucTrack::GetHitDistance() bad gap number = " << gap << endl;
00494 return kInfinity;
00495 }
00496
00497 HepPoint3D posHit = hit->GetCenterPos();
00498 HepPoint3D posHitLocal = hit->GetGap()->TransformToGap(posHit);
00499 HepPoint3D posInsctLocal = hit->GetGap()->TransformToGap(m_CurrentInsct);
00500 int orient = hit->GetGap()->Orient();
00501
00502 float distance = -9990;
00503 if(orient == 1) distance = fabs(posInsctLocal.x() - posHitLocal.x());
00504 if(orient == 0) distance = fabs(posInsctLocal.y() - posHitLocal.y());
00505
00506 return distance;
00507 }
00508
00509
00510 float
00511 RecMucTrack::GetHitDistance2(const MucRecHit* hit)
00512 {
00513 if (!hit) {
00514 cout << "RecMucTrack:GetHitDistance-E1 null hit pointer!" << endl;
00515 return kInfinity;
00516 }
00517
00518 int part = hit->Part();
00519 int gap = hit->Gap();
00520 if ( (gap < 0) || (gap >= (int)MucID::getGapNum(part)) ) {
00521 cout << "RecMucTrack::GetHitDistance() bad gap number = " << gap << endl;
00522 return kInfinity;
00523 }
00524
00525 HepPoint3D posHit = hit->GetCenterPos();
00526 HepPoint3D posHitLocal = hit->GetGap()->TransformToGap(posHit);
00527 HepPoint3D posInsctLocal = hit->GetGap()->TransformToGap(m_CurrentInsct);
00528 int orient = hit->GetGap()->Orient();
00529
00530 float distance = -9990;
00531 if(orient == 1) distance = (posInsctLocal.x() - posHitLocal.x());
00532 if(orient == 0) distance = (posInsctLocal.y() - posHitLocal.y());
00533
00534
00535
00536 return distance;
00537 }
00538
00540 Hep3Vector
00541 RecMucTrack::CalculateInsct(const int part,
00542 const int seg,
00543 const int gap)
00544 {
00545 MucGeoGap *gapPtr = MucGeoGeneral::Instance()->GetGap(part, seg, gap);
00546 vector<int> hitSeq;
00547 for(int i = 0; i < (int)m_pHits.size(); i++) {
00548 MucRecHit *aHit = m_pHits[i];
00549 if(aHit->Part() == part &&
00550 aHit->Seg() == seg &&
00551 aHit->Gap() == gap) {
00552 hitSeq.push_back(i);
00553 }
00554 }
00555 int nHitInGap = hitSeq.size();
00556
00557
00558 HepPoint3D insctLocal = gapPtr->TransformToGap(m_CurrentInsct);
00559
00560 HepPoint3D newInsctLocal = insctLocal;
00561
00562 vector<float> x;
00563 for(int i = 0; i < nHitInGap; i++) x.push_back(0);
00564 float xInsct = 0, xNewInsct = 0;
00565
00566
00567 int orient = gapPtr->Orient();
00568 if(orient == 1) xInsct = insctLocal.x();
00569 if(orient == 0) xInsct = insctLocal.y();
00570
00571 for(int i = 0; i < nHitInGap; i++) {
00572 float xStrip, yStrip, zStrip;
00573 (m_pHits[hitSeq[i]])->GetStrip()->GetCenterPos(xStrip, yStrip, zStrip);
00574 if(orient == 1) x[i] = xStrip;
00575 if(orient == 0) x[i] = yStrip;
00576 }
00577
00578
00579
00580 if(nHitInGap > 0) {
00581 xNewInsct = xInsct * kInsctWeight;
00582
00583
00584 for(int i = 0; i < nHitInGap; i++) {
00585
00586
00587 xNewInsct += x[i] * ((1.0 - kInsctWeight) / nHitInGap);
00588 }
00589
00590 if(orient == 1) {
00591 newInsctLocal.setX(xNewInsct);
00592 newInsctLocal.setY(insctLocal.y());
00593 }
00594 if(orient == 0) {
00595 newInsctLocal.setX(insctLocal.x());
00596 newInsctLocal.setY(xNewInsct);
00597 }
00598 }
00599
00600 m_CurrentInsct = gapPtr->TransformToGlobal(newInsctLocal);
00601
00602
00603
00604 return m_CurrentInsct;
00605 }
00606
00607 void
00608 RecMucTrack::AttachInsct(Hep3Vector insct)
00609 {
00610 m_Intersections.push_back(insct);
00611 }
00612
00613 void
00614 RecMucTrack::AttachDirection(Hep3Vector dir)
00615 {
00616 m_Directions.push_back(dir);
00617 }
00618
00619 void
00620 RecMucTrack::CorrectDir()
00621 {
00622
00623 m_CurrentDir = m_CurrentInsct - m_CurrentPos;
00624 m_CurrentDir.setMag(1.0);
00625
00626 }
00627
00628 void
00629 RecMucTrack::CorrectPos()
00630 {
00631 m_CurrentPos = m_CurrentInsct;
00632 }
00633
00634 void
00635 RecMucTrack::ComputeLastGap()
00636 {
00637 int lastGap1, lastGap2;
00638 int firstGap1, firstGap2;
00639 lastGap1 = lastGap2 = -1;
00640 firstGap1 = firstGap2 = 99;
00641 vector<MucRecHit*>::const_iterator iHit;
00642 iHit = m_pHits.begin();
00643 for( ;iHit != m_pHits.end(); iHit++)
00644 {
00645 if(*iHit)
00646 {
00647 int part = (*iHit)->Part();
00648 int gap = (*iHit)->Gap();
00649
00650 if(part == 1) {
00651 if(gap > lastGap1) lastGap1 = gap;
00652 if(gap < firstGap1) firstGap1 = gap;
00653 }
00654 else {
00655 if(gap > lastGap2) { m_ecPart = part; m_endPart = part; lastGap2 = gap; }
00656 if(gap < firstGap2) firstGap2 = gap;
00657 }
00658 }
00659 }
00660
00661 m_brLastLayer = lastGap1;
00662 if(firstGap1 == 99) m_brFirstLayer = -1;
00663 else if(firstGap1>=0&&firstGap1<9) m_brFirstLayer = firstGap1;
00664
00665 m_ecLastLayer = lastGap2;
00666 if(firstGap2 == 99) m_ecFirstLayer = -1;
00667 else if(firstGap2>=0&&firstGap2<8) m_ecFirstLayer = firstGap2;
00668
00669
00670 }
00671
00672 int
00673 RecMucTrack::ComputeDepth(int method)
00674 {
00675 if( m_numLayers == 0 ) {
00676 m_depth = m_depth_3 = -99; return 0;
00677 }
00678 else if( m_numLayers == 1 && (m_brLastLayer == 0 || m_ecLastLayer == 0) ) {
00679 m_depth = m_depth_3 = 0; return 0;
00680 }
00681
00682 m_depth_3 = 0.0;
00683
00684 float brThickness = 0.0; float ecThickness = 0.0; float deltaPhi = 0.0;
00685 int betweenSeg = 0;
00686
00687 float phi = m_MucMomentum.phi();
00688 float theta = m_MucMomentum.theta();
00689
00690 vector<MucRecHit*>::const_iterator iHit;
00691 int ngap = MucID::getGapMax();
00692
00693
00694 int Seg[9]; int part, seg, gap, strip;
00695 for(int gap = 0; gap < ngap; gap++){Seg[gap] = -1;}
00696 for(iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
00697 if(*iHit) {
00698 part = (*iHit)->Part();
00699 seg = (*iHit)->Seg();
00700 gap = (*iHit)->Gap();
00701 strip = (*iHit)->Strip();
00702 if(part==1) Seg[gap] = seg;
00703 }
00704 }
00705
00706 int segTrackBr = -1;
00707 for(int gap = 0; gap <= brLastLayer(); gap++){
00708 if(Seg[gap] != -1) segTrackBr = Seg[gap];
00709 }
00710
00711 for(int gap = 0; gap <= brLastLayer(); gap++){
00712 float thickness = 0.0;
00713 if(Seg[gap] != -1 && Seg[brLastLayer()-1] != -1 && Seg[gap] != Seg[brLastLayer()-1]) betweenSeg = 1;
00714 thickness = MucGeoGeneral::Instance()->GetGap(1, 0, gap)->GetIronThickness();
00715
00716 if(sin(m_MucMomentum.theta()) != 0 ) thickness /= sin(m_MucMomentum.theta());
00717 else ;
00718 deltaPhi = m_MucMomentum.phi() - segTrackBr*kPi/4;
00719 if(Seg[gap] == -1 && betweenSeg == 1) {
00720 thickness += 40;
00721 }
00722 if(cos(deltaPhi) != 0 ) thickness /= cos(deltaPhi);
00723 else ;
00724 if(deltaPhi == 0 && Seg[brLastLayer()-1]==2 ) thickness = 0;
00725
00726
00727 if(brFirstLayer()<brLastLayer()) brThickness += thickness;
00728 else if(brFirstLayer()==brLastLayer()) {
00729 if(m_MucMomentum.mag()>1000 || brFirstLayer()<2) brThickness += thickness;
00730
00731 }
00732 else cout<<"in RecMucTrack: Wrong Gap Info"<<endl;
00733
00734
00735 }
00736
00737
00738
00739
00740
00741 for (int gap = 0; gap!=-1&&gap <= ecLastLayer(); gap++) {
00742 ecThickness += MucGeoGeneral::Instance()->GetGap(0, 0, gap)->GetIronThickness();
00743 }
00744
00745 if (cos(theta) != 0.0) ecThickness /= cos(theta);
00746 else ;
00747 ecThickness = fabs(ecThickness);
00748
00749
00750
00751 if(method == 2){
00752
00753 if((m_Good3DLine == 1) &&(m_Good3DPart == 1)){
00754 if(m_IntersectionInner[0].x()!=-9999){
00755 for(int gap = 1;gap <= brLastLayer(); gap++)
00756 {
00757 if(m_IntersectionInner[gap].x() != -9999 && m_IntersectionInner[gap-1].x() != -9999)
00758 m_depth_3 += (m_IntersectionInner[gap] - m_IntersectionOuter[gap-1]).mag();
00759 }
00760
00761 m_depth_3 += ecThickness;
00762 }
00763 else m_depth_3 = brThickness + ecThickness;
00764 }
00765 if((m_Good3DLine == 1) &&(m_Good3DPart != 1)){
00766 for(int gap = 1;gap <= ecLastLayer(); gap++)
00767 {
00768 if(m_IntersectionInner[gap].x() != -9999 && m_IntersectionInner[gap-1].x() != -9999)
00769 m_depth_3 += (m_IntersectionInner[gap] - m_IntersectionOuter[gap-1]).mag();
00770 }
00771
00772 m_depth_3 += brThickness;
00773 if (cos(theta) != 0.0) m_depth_3 += 40/cos(theta);
00774 }
00775 if(m_Good3DLine == 0) m_depth_3 = brThickness + ecThickness;
00776 if(m_depth_3>2000) m_depth_3 = brThickness + ecThickness;
00777 }
00778 else
00779 {
00780 m_depth_3 = brThickness + ecThickness;
00781
00782 }
00783
00784 double offset = 50.0;
00785
00786
00787 m_depth = m_depth_3;
00788
00789
00790 if(m_depth>2000) m_depth = -99;
00791
00792
00793 }
00794
00795 void
00796 RecMucTrack::ComputeDepth()
00797 {
00798 m_depth = -99.0;
00799
00800 float brThickness = 0.0;
00801 for (int gap = 0; gap <= brLastLayer(); gap++) {
00802 brThickness += MucGeoGeneral::Instance()->GetGap(1, 0, gap)->GetIronThickness();
00803 }
00804
00805
00806 float brThickness_2 = 0.0;
00807 for (int gap = 0; gap <= brLastLayer(); gap++) {
00808 if(HasHitInGap(1,gap)){
00809 brThickness_2 += MucGeoGeneral::Instance()->GetGap(1, 0, gap)->GetIronThickness();
00810 }
00811 }
00812
00813 float brThickness_3 = 0.0;
00814 vector<MucRecHit*>::const_iterator iHit;
00815 int ngap = MucID::getGapMax();
00816 int Seg[9]; int part, seg, gap, strip;
00817 for(int gap = 0; gap < ngap; gap++){Seg[gap] = -1;}
00818 for(iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
00819 if(*iHit) {
00820 part = (*iHit)->Part();
00821 seg = (*iHit)->Seg();
00822 gap = (*iHit)->Gap();
00823 strip = (*iHit)->Strip();
00824 if(part==1) Seg[gap] = seg;
00825 }
00826 }
00827
00828 float deltaPhi_3 = 0.0;
00829 int betweenSeg = 0;
00830
00831 int segTrackBr = -1;
00832 for(int gap = 0; gap <= brLastLayer(); gap++){
00833 if(Seg[gap] != -1) segTrackBr = Seg[gap];
00834 }
00835
00836 for(int gap = 0; gap <= brLastLayer(); gap++){
00837 float thickness = 0.0;
00838 if(Seg[gap] != -1 && Seg[brLastLayer()-1] != -1 && Seg[gap] != Seg[brLastLayer()-1]) betweenSeg = 1;
00839 thickness = MucGeoGeneral::Instance()->GetGap(1, 0, gap)->GetIronThickness();
00840 if(sin(m_MucMomentum.theta()) != 0 ) thickness /= sin(m_MucMomentum.theta());
00841 else cout<<"RecMucTrack::ComputeDepth, In Barrel,theta=0?"<<endl;
00842
00843 deltaPhi_3 = m_MucMomentum.phi() - segTrackBr*kPi/4;
00844 if(Seg[gap] == -1 && betweenSeg == 1) {
00845 cout<<"between segment"<<endl;
00846 thickness += 40;
00847 }
00848
00849 if(cos(deltaPhi_3) != 0 ) thickness /= cos(deltaPhi_3);
00850 else cout<<"RecMucTrack::ComputeDepth, In Barrel,Cos(phi)=0?"<<endl;
00851
00852 if(deltaPhi_3 == 0 && Seg[brLastLayer()-1]==2 ) thickness = 0;
00853
00854 brThickness_3 += thickness;
00855
00856 }
00857
00858
00859
00860 float phi = m_MucMomentum.phi();
00861 float deltaPhi = phi - kPi/4*(int)(phi/(kPi/4));
00862 if (deltaPhi > kPi/8) deltaPhi = kPi/4 - deltaPhi;
00863 float theta = m_MucMomentum.theta();
00864
00865
00866
00867 if (sin(theta) != 0.0) brThickness /= sin(theta);
00868 else cout << "RecMucTrack::ComputeDepth, In Barrel, Track theta = 0.0 ? " << endl;
00869
00870 brThickness /= cos(deltaPhi);
00871 brThickness = fabs(brThickness);
00872
00873
00874
00875 float ecThickness = 0.0;
00876 for (int gap = 0; gap <= ecLastLayer(); gap++) {
00877 ecThickness += MucGeoGeneral::Instance()->GetGap(0, 0, gap)->GetIronThickness();
00878 }
00879
00880
00881
00882 if (cos(theta) != 0.0) ecThickness /= cos(theta);
00883 else ;
00884 ecThickness = fabs(ecThickness);
00885
00886
00887 m_depth = brThickness + ecThickness;
00888 m_depth_3 = brThickness_3 + ecThickness;
00889 cout << "Depth " << m_depth << " Depth3 = "<<m_depth_3<<endl;
00890 m_depth = m_depth_3;
00891 double offset = 50.0;
00892 if(GetTotalHits() > 0) m_depth += offset;
00893
00894 }
00895
00896 void
00897 RecMucTrack::ComputeDistanceMatch()
00898 {
00899 bool firstHitFound = false;
00900 MucGeoGap *firstGap = 0;
00901 vector<MucRecHit*>::const_iterator iHit;
00902 vector<MucRecHit*> hitsGap0;
00903 float stripLocal[3]={0.0, 0.0, 0.0};
00904 float stripGlobal[3]={0.0, 0.0, 0.0};
00905 int nStrip = 0;
00906
00907 int part, seg, gap, strip;
00908 int barrel_gap0_exist = 0;
00909
00910 for(iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
00911 if(*iHit) {
00912 part = (*iHit)->Part();
00913 seg = (*iHit)->Seg();
00914 gap = (*iHit)->Gap();
00915 strip = (*iHit)->Strip();
00916 if(!firstHitFound && gap == 0) {
00917 firstGap = MucGeoGeneral::Instance()->GetGap(part, seg, gap);
00918 firstHitFound = true;
00919 }
00920 if(firstGap && part == firstGap->Part() && seg == firstGap->Seg() && gap == firstGap->Gap()) {
00921
00922 HepPoint3D posHit = (*iHit)->GetCenterPos();
00923 HepPoint3D posHitLocal = (*iHit)->GetGap()->TransformToGap(posHit);
00924 if(part==1&&gap==0) barrel_gap0_exist = 1;
00925
00926 stripLocal[0] += posHitLocal.x();
00927 stripLocal[1] += posHitLocal.y();
00928 stripLocal[2] += posHitLocal.z();
00929
00930 stripGlobal[0] += posHit.x();
00931 stripGlobal[1] += posHit.y();
00932 stripGlobal[2] += posHit.z();
00933
00934 nStrip++;
00935 }
00936 }
00937 }
00938
00939
00940
00941 int apart = -1, aseg = -1;
00942 int nHits = FindSegWithMaxHits(apart, aseg);
00943 MucGeoGap *fakefirstGap = 0;
00944 if(apart == -1 && aseg== -1)
00945 {m_Dist_muc_ext.set(0,0,0);}
00946 else {
00947 fakefirstGap = MucGeoGeneral::Instance()->GetGap(apart, aseg, 0);
00948 HepPoint3D fextLocal = fakefirstGap->TransformToGap(m_ExtMucPos);
00949 HepPoint3D fmucLocal = fakefirstGap->TransformToGap(m_MucPos);
00950 float dist_x = fextLocal.x() - fmucLocal.x();
00951 float dist_y = fextLocal.y() - fmucLocal.y();
00952 float dist_z = fextLocal.z() - fmucLocal.z();
00953 m_Dist_muc_ext.set(dist_x,dist_y,dist_z);
00954
00955 if (fakefirstGap->Orient() == 0) {
00956
00957 }
00958 else {
00959
00960 }
00961 }
00962
00963 float distance = -9990;
00964
00965 if (nStrip == 0 || !firstGap) {
00966 distance = -9990;
00967
00968 }
00969 else{
00970 for (int k = 0; k < 3; k++) stripLocal[k] /= nStrip;
00971 for (int k = 0; k < 3; k++) stripGlobal[k] /= nStrip;
00972
00973 m_StripPhi.set(stripGlobal[0],stripGlobal[1],stripGlobal[2]);
00974
00975 HepPoint3D extLocal = firstGap->TransformToGap(m_ExtMucPos);
00976
00977
00978
00979
00980
00981
00982 if (firstGap->Orient() == 0) {
00983 distance = extLocal.y() - stripLocal[1];
00984
00985 }
00986 else {
00987 distance = extLocal.x() - stripLocal[0];
00988
00989 }
00990 }
00991
00992 m_distance = distance;
00993
00994
00995
00996
00997 }
00998
00999 void
01000 RecMucTrack::Extend()
01001 {
01002 float x = m_ExtMucPos.x(), y = m_ExtMucPos.y(), z = m_ExtMucPos.z();
01003
01004 float step = 0.005;
01005 float myMucR = 1720.0;
01006 float myMucZ = 2110.0;
01007 int count_=0;
01008 while ( ( (fabs(x)<=myMucR) && (fabs(y)<=myMucR) && ((fabs(x-y)/sqrt(2.))<=myMucR) && ((fabs(x+y)/sqrt(2.))<=myMucR) && (fabs(z)<=myMucZ) )&&count_<1000 ) {
01009 x += step * m_MucMomentum.x();
01010 y += step * m_MucMomentum.y();
01011 z += step * m_MucMomentum.z();
01012 count_++;
01013
01014 }
01015
01016 int count = 0;
01017 while( !( (fabs(x)<=myMucR) && (fabs(y)<=myMucR) && ((fabs(x-y)/sqrt(2.))<=myMucR) && ((fabs(x+y)/sqrt(2.))<=myMucR) && (fabs(z)<=myMucZ) )&&count<1000 ) {
01018 x -= step * m_MucMomentum.x();
01019 y -= step * m_MucMomentum.y();
01020 z -= step * m_MucMomentum.z();
01021 count++;
01022
01023 }
01024
01025 if(count<1000&&count_<1000){
01026 if(fabs(x)<2600&&fabs(y)<2600&&fabs(z)<2800){
01027 m_ExtMucPos.set(x, y, z);
01028 m_MucPos.set(x,y,z);
01029 m_xPos = x;
01030 m_yPos = y;
01031 m_zPos = z;
01032 }
01033 }
01034
01035 }
01036
01037 void
01038 RecMucTrack::LineFit(int fittingMethod)
01039 {
01040 Extend();
01041
01042 int part = -1, seg = -1;
01043 int nHits = FindSegWithMaxHits(part, seg);
01044
01045
01046 if (part < 0 || seg < 0) {
01047 m_depth = 0;
01048
01049 return;
01050 }
01051
01052 float startPos[3] = {0.0, 0.0, 0.0};
01053 MucRec2DRoad road2D[2];
01054 vector<MucRecHit*>::const_iterator iHit;
01055 for (int orient = 0; orient <= 1; orient++) {
01056 road2D[orient] = MucRec2DRoad(part, seg, orient, startPos[0], startPos[1], startPos[2],fittingMethod );
01057 for(iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01058 if(*iHit) {
01059 int hitPart = (*iHit)->Part();
01060 int hitSeg = (*iHit)->Seg();
01061 int hitOrient = (*iHit)->GetGap()->Orient();
01062 if (hitPart == part && hitSeg == seg && hitOrient == orient) {
01063 road2D[orient].AttachHit(*iHit);
01064 }
01065 }
01066 }
01067
01068
01069
01070
01071
01072
01073 }
01074 MucRec3DRoad road3D(&road2D[0], &road2D[1]);
01075 float vx, vy, x0, y0;
01076 if ( road3D.GetPart() == 1 ){
01077 road3D.TransformPhiRToXY( road2D[1].GetSlope(), road2D[0].GetSlope(),
01078 road2D[1].GetIntercept(), road2D[0].GetIntercept(),
01079 vx, vy, x0, y0);
01080 }
01081 else {
01082 vx = road2D[1].GetSlope();
01083 x0 = road2D[1].GetIntercept();
01084 vy = road2D[0].GetSlope();
01085 y0 = road2D[0].GetIntercept();
01086 }
01087
01088
01089
01090
01091
01092 if (road2D[0].GetNGapsWithHits() >= 2 && road2D[1].GetNGapsWithHits() >= 2) {
01093
01094
01095 m_Good3DLine = 1;
01096 m_status = 1;
01097
01098 m_Good3DPart = road3D.GetPart();
01099 road3D.SetSimpleFitParams(vx, x0, vy, y0);
01100
01101
01102 float startx = 0.0, starty = 0.0, startz = 0.0;
01103 float startxSigma = 0.0, startySigma = 0.0, startzSigma = 0.0;
01104 float x1 = 0.0, y1 = 0.0, z1 = 0.0, x2 = 0.0, y2 = 0.0, z2 = 0.0;
01105 int gap = 0;
01106
01107 if(fittingMethod==2){
01108 for(int igap=0; igap<9;igap++){
01109 road3D.Project(igap, x1, y1, z1, x2, y2, z2);
01110 m_IntersectionInner[igap].set(x1,y1,z1);
01111 m_IntersectionOuter[igap].set(x2,y2,z2);
01112
01113 }
01114 }
01115
01116 do {
01117
01118 road3D.ProjectWithSigma(gap, startx, starty, startz, startxSigma, startySigma, startzSigma);
01119 gap++;
01120 }
01121 while ( (startx*startx + starty*starty + startz*startz) < kMinor &&
01122 gap < (int)MucID::getGapNum(part) );
01123
01124 if(fabs(startx)<2600&&fabs(starty)<2600&&fabs(startz)<2800){
01125 Hep3Vector MucPos_self;
01126 MucPos_self.set(startx, starty, startz);
01127 float dist = (MucPos_self - m_MucPos).mag();
01128
01129 if(dist < 1000 ){
01130 m_MucPos.set(startx, starty, startz);
01131 m_xPos = startx;
01132 m_yPos = starty;
01133 m_zPos = startz;
01134 }
01135 }
01136 m_MucPosSigma.set(startxSigma, startySigma, startzSigma);
01137 m_xPosSigma = startxSigma;
01138 m_yPosSigma = startySigma;
01139 m_zPosSigma = startzSigma;
01140
01141
01142 float momentum = m_MucMomentum.mag();
01143 float zDir = 1.0;
01144 if (m_MucMomentum.z() != 0.0) zDir = m_MucMomentum.z() / fabs(m_MucMomentum.z());
01145
01146 float px = road3D.GetSlopeZX()*zDir;
01147 float py = road3D.GetSlopeZY()*zDir;
01148 float pz = zDir;
01149 float segPhi = 0.25*kPi*seg;
01150
01151 if (part == 1 && px*cos(segPhi)+py*sin(segPhi) < 0.0) {
01152 px *= -1.0;
01153 py *= -1.0;
01154 pz *= -1.0;
01155 }
01156
01157
01158 m_MucMomentum.set(px, py, pz);
01159 m_MucMomentum.setMag(momentum);
01160
01161 m_px = px; m_py = py; m_pz = pz;
01162 Hep3Vector phi_mdc, phi_muc;
01163 phi_mdc.set(m_MdcMomentum.x(),m_MdcMomentum.y(),0);
01164 phi_muc.set(px,py,0);
01165 double deltaPhi = phi_mdc.angle(phi_muc);
01166
01167 m_deltaPhi = deltaPhi;
01168 m_dof = road3D.GetDegreesOfFreedom();
01169 m_chi2 = road3D.GetReducedChiSquare();
01170 if(m_chi2<0||m_chi2>1000) m_chi2 = 0;
01171 m_rms = 0.0;
01172
01173
01174 if(road2D[0].GetNGapsWithHits() >= 3 && road2D[1].GetNGapsWithHits() >= 3){
01175
01176
01177 float xInsct = 0.0, yInsct = 0.0, zInsct = 0.0, sigmax = 0.0, sigmay = 0.0, sigmaz = 0.0;
01178 float quad_x1 = 0.0, quad_y1 = 0.0, quad_z1 = 0.0, quad_x2 = 0.0, quad_y2 = 0.0, quad_z2 = 0.0;
01179 for(int ihit = 0; ihit < m_pHits.size(); ihit++){
01180 int igap = (m_pHits[ihit])->Gap();
01181 road3D.ProjectNoCurrentGap(igap, xInsct, yInsct, zInsct, sigmax, sigmay, sigmaz,
01182 quad_x1, quad_y1, quad_z1, quad_x2, quad_y2, quad_z2);
01183
01184 SetCurrentInsct(xInsct, yInsct, zInsct);
01185
01186 float dist_hit_track = GetHitDistance2(m_pHits[ihit]);
01187 m_distHits.push_back(dist_hit_track);
01188
01189
01190
01191 if(fittingMethod==2)
01192 {
01193
01194 Hep3Vector center = m_pHits[ihit]->GetCenterPos();
01195 int iPart = m_pHits[ihit]->Part();
01196 int iSeg = m_pHits[ihit]->Seg();
01197 int iGap = m_pHits[ihit]->Gap();
01198 float xHit, yHit, zHit = 0.;
01199 if (iPart == 1) {
01200 if ( iGap %2 == 1) {
01201 xHit = center.z();
01202 yHit = sqrt(center.x()*center.x()
01203 + center.y()*center.y());
01204 if(iSeg==2) yHit = center.y();
01205 }
01206 else {
01207 xHit = center.x();
01208 yHit = center.y();
01209 }
01210 }
01211 else {
01212 if ( iGap %2 == 0) {
01213 xHit = center.z();
01214 yHit = center.y();
01215 }
01216 else {
01217 xHit = center.z();
01218 yHit = center.x();
01219 }
01220 }
01221
01222 float distance1 = fabs(xHit - quad_x1)/(xHit - quad_x1) * sqrt((xHit - quad_x1)*(xHit - quad_x1) + (yHit - quad_y1)*(yHit - quad_y1));
01223 float distance2 = fabs(xHit - quad_x2)/(xHit - quad_x2) * sqrt((xHit - quad_x2)*(xHit - quad_x2) + (yHit - quad_y2)*(yHit - quad_y2));
01224
01225 float dist_quad = distance1;
01226 if(fabs(distance1) > fabs(distance2)) dist_quad = distance2;
01227
01228
01229
01230
01231
01232 if(quad_x1 == -9999) m_distHits_quad.push_back(-99);
01233 else m_distHits_quad.push_back(dist_quad);
01234
01235 }
01236 }
01237 }
01238 else{
01239 for(int ihit = 0; ihit < m_pHits.size(); ihit++){
01240 m_distHits.push_back(-99);
01241 m_distHits_quad.push_back(-99);
01242 }
01243 }
01245
01247 ComputeLastGap();
01248
01249 HepPoint3D initPos(startx-m_MucMomentum.x()/momentum, starty-m_MucMomentum.y()/momentum,startz-m_MucMomentum.z()/momentum);
01250
01251
01252 for(int igap = 0; igap < 9 ; igap++){
01253 vector<int> padID;
01254 vector<float> intersect_x;
01255 vector<float> intersect_y;
01256 vector<float> intersect_z;
01257
01258
01259 float zx, zy, x0, y0;
01260 int status = road3D.RefitNoCurrentGap(igap,zx,x0,zy,y0);
01261 Hep3Vector mom_refit;
01262 if(status == 0){
01263 float zDir = 1.0;
01264 if (m_MucMomentum.z() != 0.0) zDir = m_MucMomentum.z() / fabs(m_MucMomentum.z());
01265 float px_refit = zx * zDir;
01266 float py_refit = zy * zDir;
01267 float pz_refit = zDir;
01268 float segPhi = 0.25*kPi*seg;
01269
01270 if (part == 1 && px_refit*cos(segPhi)+py_refit*sin(segPhi) < 0.0) {
01271 px_refit *= -1.0;
01272 py_refit *= -1.0;
01273 pz_refit *= -1.0;
01274 }
01275
01276 mom_refit.setX(px_refit);
01277 mom_refit.setY(py_refit);
01278 mom_refit.setZ(pz_refit);
01279 mom_refit.setMag(momentum);
01280 }
01281 else mom_refit = m_MucMomentum;
01282
01283 HepPoint3D initPos_refit;
01284 if(status == 0){
01285 float initPosx_refit, initPosy_refit, initPosz_refit;
01286 float sigmax_refit, sigmay_refit, sigmaz_refit;
01287 road3D.Project(0, zx, x0, zy, y0, initPosx_refit, initPosy_refit, initPosz_refit);
01288 initPos_refit.setX(initPosx_refit-mom_refit.x()/momentum);
01289 initPos_refit.setY(initPosy_refit-mom_refit.y()/momentum);
01290 initPos_refit.setZ(initPosz_refit-mom_refit.z()/momentum);
01291 }
01292 else initPos_refit = initPos;
01293
01294
01295
01296
01297
01298
01299
01300 vector<Identifier> MuId = road3D.ProjectToStrip(1,igap,initPos_refit,mom_refit,padID,intersect_x,intersect_y,intersect_z);
01301
01302
01303
01304 vector<Identifier>::const_iterator mucid;
01305 int i = 0;
01306 for(mucid = MuId.begin(); mucid != MuId.end(); mucid++,i++)
01307 {
01308 MucRecHit *pHit = new MucRecHit(MucID::part(*mucid),MucID::seg(*mucid),MucID::gap(*mucid),MucID::strip(*mucid));
01309 pHit->SetPadID(padID[i]);
01310 pHit->SetIntersectX(intersect_x[i]);
01311 pHit->SetIntersectY(intersect_y[i]);
01312 pHit->SetIntersectZ(intersect_z[i]);
01313
01314 m_pExpectedHits.push_back(pHit);
01315 }
01316 }
01317
01318 if(m_endPart!=-1&&m_ecLastLayer!=-1)
01319 {
01320
01321 for(int igap = 0; igap < 8; igap++){
01322 vector<int> padID;
01323 vector<float> intersect_x;
01324 vector<float> intersect_y;
01325 vector<float> intersect_z;
01326
01327
01328 float zx, zy, x0, y0;
01329
01330 int status = road3D.RefitNoCurrentGap(igap,zx,x0,zy,y0);
01331 Hep3Vector mom_refit;
01332 if(status == 0){
01333 float zDir = 1.0;
01334 if (m_MucMomentum.z() != 0.0) zDir = m_MucMomentum.z() / fabs(m_MucMomentum.z());
01335 float px_refit = zx * zDir;
01336 float py_refit = zy * zDir;
01337 float pz_refit = zDir;
01338 float segPhi = 0.25*kPi*seg;
01339
01340 if (part == 1 && px_refit*cos(segPhi)+py_refit*sin(segPhi) < 0.0) {
01341 px_refit *= -1.0;
01342 py_refit *= -1.0;
01343 pz_refit *= -1.0;
01344 }
01345
01346 mom_refit.setX(px_refit);
01347 mom_refit.setY(py_refit);
01348 mom_refit.setZ(pz_refit);
01349 mom_refit.setMag(momentum);
01350 }
01351 else mom_refit = m_MucMomentum;
01352
01353 HepPoint3D initPos_refit;
01354 if(status == 0){
01355 float initPosx_refit, initPosy_refit, initPosz_refit;
01356 float sigmax_refit, sigmay_refit, sigmaz_refit;
01357 road3D.Project(0, zx, x0, zy, y0, initPosx_refit, initPosy_refit, initPosz_refit);
01358 initPos_refit.setX(initPosx_refit-mom_refit.x()/momentum);
01359 initPos_refit.setY(initPosy_refit-mom_refit.y()/momentum);
01360 initPos_refit.setZ(initPosz_refit-mom_refit.z()/momentum);
01361 }
01362 else initPos_refit = initPos;
01363
01364
01365
01366
01367 vector<Identifier> MuId = road3D.ProjectToStrip(m_endPart,igap,initPos_refit,mom_refit,padID,intersect_x,intersect_y,intersect_z);
01368
01369 vector<Identifier>::const_iterator mucid;
01370 int i = 0;
01371 for(mucid = MuId.begin(); mucid != MuId.end(); mucid++,i++)
01372 {
01373 MucRecHit *pHit = new MucRecHit(MucID::part(*mucid),MucID::seg(*mucid),MucID::gap(*mucid),MucID::strip(*mucid));
01374 pHit->SetPadID(padID[i]);
01375 pHit->SetIntersectX(intersect_x[i]);
01376 pHit->SetIntersectY(intersect_y[i]);
01377 pHit->SetIntersectZ(intersect_z[i]);
01378
01379 m_pExpectedHits.push_back(pHit);
01380 }
01381 }
01382 }
01383
01384
01385 for(int i=0; i< m_pExpectedHits.size(); i++)
01386 {
01387 MucRecHit *ihit = m_pExpectedHits[i];
01388
01389
01390 }
01391
01392 for(int i=0; i< m_pHits.size(); i++)
01393 {
01394 MucRecHit *ihit = m_pHits[i];
01395
01396 }
01397
01398
01399 }
01400 else {
01401
01402
01403 for(int ihit = 0; ihit < m_pHits.size(); ihit++){
01404 m_distHits.push_back(-99);
01405 m_distHits_quad.push_back(-99);
01406 }
01407 }
01408
01409 }
01410
01411 void
01412 RecMucTrack::ComputeNGapsWithHits()
01413 {
01414 int ngap = MucID::getGapMax();
01415 int gap, count = 0;
01416 vector<int> firedGap;
01417 for(gap = 0; gap < ngap; gap++) firedGap.push_back(0);
01418
01419 vector<MucRecHit*>::const_iterator iHit;
01420 for (iHit=m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01421 if (*iHit) {
01422 gap = (*iHit)->Gap();
01423 firedGap[gap] = 1;
01424 }
01425 }
01426
01427 for ( gap = 0; gap < ngap; gap++) {
01428 count += firedGap[gap];
01429 }
01430
01431 m_numHits = m_pHits.size();
01432 m_numLayers = count;
01433
01434 }
01435
01436 int
01437 RecMucTrack::GetTotalHits() const
01438 {
01439 return m_pHits.size();
01440 }
01441
01442 int
01443 RecMucTrack::GetHitInGap(const int part, const int gap) const
01444 {
01445 if ( part < 0 || part > 2 ) {
01446 cout << "RecMucTrack::GetHitInGap(), invalid part " << part << endl;
01447 return 0;
01448 }
01449 if ( (gap < 0) || (gap >= (int)MucID::getGapNum(part)) ) {
01450 cout << "RecMucTrack::GetHitInGap(), invalid gap " << gap << endl;
01451 return 0;
01452 }
01453
01454 vector<MucRecHit*>::const_iterator iHit;
01455 int hitsInGap = 0;
01456
01457 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01458
01459 if ( !(*iHit) ) {
01460 cout << "RecMucTrack::GetHitInGap() null hit pointer !" << endl;
01461 return 0;
01462 }
01463 else {
01464 if( part == (*iHit)->Part() && gap == (*iHit)->Gap() ) hitsInGap++;
01465 }
01466 }
01467
01468 return hitsInGap;
01469 }
01470
01471 int
01472 RecMucTrack::GetHitInSeg(const int part, const int seg) const
01473 {
01474 int num = GetHitInSegOrient(part, seg, 0) + GetHitInSegOrient(part, seg, 1);
01475 return num;
01476 }
01477
01478 int
01479 RecMucTrack::GetHitInSegOrient(const int part, const int seg, const int orient) const
01480 {
01481 if ( part < 0 || part > 2 ) {
01482 cout << "RecMucTrack::GetHitInSeg(), invalid part " << part << endl;
01483 return 0;
01484 }
01485 if ( (seg < 0) || (seg >= (int)MucID::getSegNum(part)) ) {
01486 cout << "RecMucTrack::GetHitInSeg(), invalid seg = " << seg << endl;
01487 return 0;
01488 }
01489
01490 vector<MucRecHit*>::const_iterator iHit;
01491 int hitsInSeg = 0;
01492
01493 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01494
01495 if ( !(*iHit) ) {
01496 cout << "RecMucTrack::GetHitInSeg(), null hit pointer !" << endl;
01497 return 0;
01498 }
01499 else {
01500 if( part == (*iHit)->Part() && seg == (*iHit)->Seg() &&
01501 orient == (*iHit)->GetGap()->Orient() ) hitsInSeg++;
01502 }
01503 }
01504
01505 return hitsInSeg;
01506 }
01507
01508 int
01509 RecMucTrack::FindSegWithMaxHits(int &findPart, int &findSeg)
01510 {
01511 int maxHits = 0;
01512 int hits = 0;
01513 for(int part = 0; part < (int)MucID::getPartNum(); part++) {
01514 for(int seg = 0; seg < (int)MucID::getSegNum(part); seg++) {
01515 hits = GetHitInSeg(part, seg);
01516 if (hits > maxHits) {
01517 maxHits = hits;
01518 findPart = part;
01519 findSeg = seg;
01520 }
01521 }
01522 }
01523
01524 return maxHits;
01525 }
01526
01527 void
01528 RecMucTrack::ComputeMaxHitsInGap()
01529 {
01530 int maxHits = 0;
01531 int hits = 0;
01532 for(int part = 0; part < (int)MucID::getPartNum(); part++) {
01533 for(int gap = 0; gap < (int)MucID::getGapNum(part); gap++) {
01534 hits = GetHitInGap(part, gap);
01535 if(hits > maxHits) maxHits = hits;
01536 }
01537 }
01538
01539 m_maxHitsInLayer = maxHits;
01540 }
01541
01542 bool
01543 RecMucTrack::HasHit(const int part, const int seg, const int gap, const int strip) const
01544 {
01545 bool found = false;
01546 vector<MucRecHit*>::const_iterator iHit;
01547
01548 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01549 if (*iHit) {
01550 if ( (*iHit)->Part() == part &&
01551 (*iHit)->Seg() == seg &&
01552 (*iHit)->Gap() == gap &&
01553 (*iHit)->Strip() == strip ) {
01554 found = true;
01555 }
01556 }
01557 }
01558
01559 return found;
01560 }
01561
01562 bool
01563 RecMucTrack::HasHitInGap(const int part, const int gap) const
01564 {
01565 if ( (gap < 0) || (gap >= (int)MucID::getGapNum(part)) ) {
01566 cout << "RecMucTrack::HasHitInGap-E2 invalid gap = " << gap << endl;
01567 return false;
01568 }
01569
01570 bool found = false;
01571 vector<MucRecHit*>::const_iterator iHit;
01572
01573 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01574 if (*iHit) {
01575 if ( (*iHit)->Part() == part &&
01576 (*iHit)->Gap() == gap ) {
01577 found = true;
01578 }
01579 }
01580 }
01581
01582 return found;
01583 }
01584
01585 int
01586 RecMucTrack::GetNSharedHits(const RecMucTrack* track2) const
01587 {
01588 if (!track2) {
01589 return 0;
01590 }
01591
01592 int count = 0;
01593 vector<MucRecHit*>::const_iterator iHit1;
01594 vector<MucRecHit*>::const_iterator iHit2;
01595 MucRecHit *hit1, *hit2;
01596
01597 for( iHit1 = m_pHits.begin(); iHit1 != m_pHits.end(); iHit1++){
01598 for( iHit2 = track2->m_pHits.begin();
01599 iHit2 != track2->m_pHits.end(); iHit2++){
01600 hit1 = (*iHit1);
01601 hit2 = (*iHit2);
01602
01603 if ( (hit1 != 0) && (hit2 != 0) ) {
01604 if (hit1->GetID() == hit2->GetID()) {
01605 count++;
01606 }
01607 }
01608 }
01609 }
01610
01611 return count;
01612 }
01613
01614 MucRecHit*
01615 RecMucTrack::GetHit(const int part, const int gap) const
01616 {
01617 if ( (gap < 0) || (gap >= (int)MucID::getGapNum(part)) ) {
01618 cout << "RecMucTrack::Hit-E1 invalid gap = " << gap << endl;
01619 return 0;
01620 }
01621
01622 vector<MucRecHit*>::const_iterator iHit;
01623
01624 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01625 if (*iHit) {
01626 if ( (*iHit)->Part() == part &&
01627 (*iHit)->Gap() == gap) {
01628 return (*iHit);
01629 }
01630 }
01631 }
01632
01633 return 0L;
01634 }
01635
01636 vector<MucRecHit*>
01637 RecMucTrack::GetHits() const
01638 {
01639 return m_pHits;
01640 }
01641
01642 vector<MucRecHit*>
01643 RecMucTrack::GetExpectedHits() const
01644 {
01645 return m_pExpectedHits;
01646 }
01647
01648 vector<long>
01649 RecMucTrack::GetHitIndices() const
01650 {
01651 vector<long> v;
01652
01653 vector<MucRecHit*>::const_iterator iHit;
01654 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01655 if (*iHit) {
01656 long index = (*iHit)->GetID();
01657 v.push_back(index);
01658
01659
01660
01661
01662
01663
01664 }
01665 }
01666 return v;
01667 }
01668
01669 void
01670 RecMucTrack::ComputeTrackInfo(int fittingmethod)
01671 {
01672 setVecHits(m_pHits);
01673 setExpHits(m_pExpectedHits);
01674 ComputeLastGap();
01675 ComputeNGapsWithHits();
01676 ComputeMaxHitsInGap();
01677
01678 ComputeDepth(fittingmethod);
01679 ComputeDistanceMatch();
01680 OutputUnitChange();
01681 }
01682
01683 void
01684 RecMucTrack::PrintHitsInfo() const
01685 {
01686 vector<MucRecHit*>::const_iterator iHit;
01687 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01688 if (*iHit) {
01689 float xl, yl, zl;
01690 (*iHit)->GetStrip()->GetCenterPos(xl, yl, zl);
01691 HepPoint3D vl(xl, yl, zl);
01692 HepPoint3D vg = (*iHit)->GetGap()->TransformToGlobal(vl);
01693
01694 cout << " part " << (*iHit)->Part()
01695 << " seg " << (*iHit)->Seg()
01696 << " gap " << (*iHit)->Gap()
01697 << " strip " << (*iHit)->Strip()
01698 << " pos (" << vg.x() << ", " << vg.y() << ", " << vg.z() << ")"
01699 << endl;
01700 }
01701 }
01702
01703 }
01704
01705 void
01706 RecMucTrack::OutputUnitChange()
01707 {
01708 if(m_changeUnit == false){
01709 m_depth/=10;
01710
01711 m_xPos /=10;
01712 m_yPos /=10;
01713 m_zPos /=10;
01714
01715 m_xPosSigma/=10;
01716 m_yPosSigma/=10;
01717 m_zPosSigma/=10;
01718
01719 m_px /= 1000;
01720 m_py /= 1000;
01721 m_pz /= 1000;
01722
01723 m_distance /= 10;
01724 }
01725
01726 m_changeUnit = true;
01727
01728 }